Vivette Girault

Béatrice Rivière

Loïc Cappanera

Loïc Cappanera

Convergence

Beatrice M Rivi Ère

AND Loic Cappanera

CONVERGENCE OF A FINITE ELEMENT METHOD FOR DEGENERATE TWO-PHASE FLOW IN POROUS MEDIA

whether they are published or not. The documents may come

1. Introduction. This work discretizes on a suitable mesh a degenerate two-phase flow system set in a polyhedral domain by a new finite element scheme that directly approximates the wetting phase pressure and saturation, similar to the formulation proposed in [?]. Mass lumping is used to compute the integrals and a suitable upwinding is used to compute the flux, guaranteeing that the discrete saturation satisfies a maximum principle. The resulting system of discrete equations is a finite element analogue of the finite volume scheme introduced and analyzed by Eymard et al in the seminal work [START_REF] Eymard | Mathematical study of a petroleum-engineering scheme[END_REF]. From the point of view of implementation, the advantage of finite elements is that they only use nodal values and a single simplicial mesh. In particular, no orthogonality property is required between the faces and the lines joining the centers of control volumes. From a theoretical point of view, owing that the finite element scheme is based on functions, some steps in its convergence analysis are simpler, but nevertheless the major difficulty in the analysis consists in proving sufficient a priori estimates in spite of the degeneracy. By following closely [START_REF] Eymard | Mathematical study of a petroleum-engineering scheme[END_REF], the degeneracy is remediated by reintroducing in the proofs discrete artificial pressures. From there, convergence of the numerical solutions is shown via a compactness argument.

Incompressible two-phase flow is a popular and important multiphase flow model in reservoirs for the oil and gas industry. Based on conservation laws at the continuum scale, the model assumes the existence of a representative elementary volume. Each wetting phase and non-wetting phase saturation satisfies a mass balance equation and each phase velocity follows the generalized Darcy law [START_REF] Peaceman | Fundamentals of numerical reservoir simulation[END_REF][START_REF] Aziz | Petroleum reservoir simulation[END_REF]. The equations of the mathematical model read

∂ t (ϕs w ) -∇ • (η w ∇p w ) = f w (s in )q -f w (s w )q, ∂ t (ϕs o ) -∇ • (η o ∇p o ) = f o (s in )q -f o (s w )q, p c (s w ) = p o -p w , s w + s o = 1, (1.1)
complemented by initial and boundary conditions. Here p w , s w , η w , f w , (respectively, p o , s o , η o , f o ), are the pressure, saturation, mobility and fractional flow of the wetting (respectively non-wetting) phase, ϕ is the porosity, s in is a given input saturation, and q, q are given flow rates. The capillary pressure, p c , is a given function that depends nonlinearly on the saturation. Because the phase mobilities are degenerate when they are evaluated at the residual phase saturations and the derivative of the capillary pressure is unbounded, this system of two coupled nonlinear partial differential equations has coefficients that vanish in parts of the domain; this degeneracy makes the numerical analysis challenging.

At the continuous level, this problem has several equivalent formulations, see [7]. They are linked to the choice of primary unknowns selected among wetting phase and non-wetting phase pressure and saturation, or capillary pressure [START_REF] Helmig | Multiphase flow and transport processes in the subsurface: a contribution to the modeling of hydrosystems[END_REF][START_REF] Bastian | A fully-coupled discontinuous Galerkin method for two-phase flow in porous media with discontinuous capillary pressure[END_REF]. A good state of the art can be found in the reference [START_REF] Arbogast | The existence of weak solutions to single porosity and simple dual-porosity models of two-phase incompressible flow[END_REF]. Up to our knowledge, the mathematical analysis of the system of equations was first done in [START_REF] Kroener | Flow of oil and water in a porous medium[END_REF][START_REF] Alt | Nonsteady flow of water and oil through inhomogeneous porous media[END_REF]. An equivalent formulation of the model, based on Chavent's global pressure that removes the degeneracy, was analyzed in [START_REF] Chen | Degenerate two-phase incompressible flow. I. Existence, uniqueness and regularity of a weak solution[END_REF][START_REF] Chen | Mathematical analysis for reservoir models[END_REF]. Since then, the global pressure formulation has been discretized and analyzed in many references [START_REF] Ohlberger | Convergence of a mixed finite element: Finite volume method for the two phase flow in porous media[END_REF][START_REF] Chen | Degenerate two-phase incompressible flow III. Sharp error estimates[END_REF][START_REF] Michel | A finite volume scheme for two-phase incompressible flow in porous media[END_REF], but unfortunately, this formulation is not used in engineering practice because the global pressure is not a physical unknown. Otherwise, with one exception, the numerical analysis of the discrete version of (1.1), has always been done under unrealistic assumptions that cannot be checked at the discrete level [START_REF] Douglas | Finite difference methods for two-phase incompressible flow in porous media[END_REF][START_REF] Epshteyn | Analysis of hp discontinuous Galerkin methods for incompressible two-phase flow[END_REF]. Related to this line of work, the discretization of a degenerate parabolic equation has been studied in the literature [START_REF] Arbogast | A nonlinear mixed finite element method for a degenerate parabolic equation arising in flow in porous media[END_REF][START_REF] Yotov | A mixed finite element discretization on non-matching multiblock grids for a degenerate parabolic equation arising in porous media flow[END_REF][START_REF] Woodward | Analysis of expanded mixed finite element methods for a nonlinear parabolic equation modeling flow into variably saturated porous media[END_REF][START_REF] Eymard | A combined finite volume-nonconforming/mixedhybrid finite element scheme for degenerate parabolic problems[END_REF]. The only paper that performs the complete numerical analysis of the discrete degenerate two-phase flow system written as above (i.e., in the form used by engineers) is the analysis on finite volumes done in reference [START_REF] Eymard | Mathematical study of a petroleum-engineering scheme[END_REF]. This motivates our extension of this work to finite elements.

The remaining part of this introduction makes precise problem (1.1). The numerical scheme is developed in Section 2. Because of the nonlinearity and degeneracy of its equations, existence of a discrete solution requires that the discrete wetting phase saturation satisfies a maximum principle. This is the first object of Section 3, the second one being basic a priori pressure estimates, after which existence is shown in Section 4. The most technical part, done in Section 5, is the derivation of an unconditional bound on an auxiliary pressure, which allows to use a compactness argument. Weak and strong convergences are proved in Section 6 and the equations satisfied by the limit are identified in Section 7, thus confirming existence of a solution of the weak formulation (1.16). Numerical results are presented in Section 8.

Model problem.

Let Ω ⊂ IR d , d = 2 or 3, be a bounded connected Lipschitz domain with boundary ∂Ω and unit exterior normal n, and let T be a final time. With the last relation in (1.1), s w is the only unknown saturation; so we set s = s w , and rewrite (1.1) almost everywhere in Ω×]0, T [ as

∂ t (ϕs) -∇ • (η w ∇p w ) = f w (s in )q -f w (s)q (1.2) -∂ t (ϕs) -∇ • (η o ∇p o ) = f o (s in )q -f o (s)q, ( 1.3)
complemented by a natural boundary condition almost everywhere on ∂Ω×]0, T [

η w ∇p w • n = 0, η o ∇p o • n = 0, (1.4)
and an initial condition almost everywhere in Ω s w (•, 0) = s 0 w , 0 ≤ s 0 w ≤ 1.

(1.5)

The fractional flows are related to the mobilities by

f w = η w η w + η o , f o = 1 -f w . (1.6)
Recall that the phase saturations sum up to 1 and the phase pressures are related by p c (s w ) = p o -p w .

(1.7)

The first part of this work, up to Section 5.7, is done under the following basic assumptions:

Assumptions

• The porosity ϕ is piecewise constant in space, independent of time, positive, bounded and uniformly bounded away from zero. • The mobility of the wetting phase η w ≥ 0 is continuous and increasing. The mobility of the non-wetting phase η o ≥ 0 is continuous and decreasing. This implies that the function f w is increasing and the function f o is decreasing. • There is a positive constant η * such that η w (s) + η o (s) ≥ η * , ∀s ∈ [0, 1].

(1.8)

• The capillary pressure p c is a continuous, strictly decreasing function in W 1,1 (0, 1).

• The flow rates at the injection and production wells, q, q ∈ L 2 (Ω×]0, T [) satisfy q ≥ 0, q ≥ 0, Ω q = Ω q.

(1.9)

• The prescribed input saturation s in satisfies almost everywhere in Ω×]0, T [ 0 ≤ s in ≤ 1.

(1.10)

Since p c , η α , f α , α = w, o are bounded above and below, it is convenient to extend them continuously by constants to IR.

Although the numerical scheme studied below does not discretize the global pressure, following [START_REF] Eymard | Mathematical study of a petroleum-engineering scheme[END_REF], its convergence proof uses a number of auxiliary functions related to the global pressure. First, we introduce the primitive g c of p c , ∀x ∈ [0, 1], g c (x) = Moreover, the derivative of g satisfies formally the identities ∀x ∈ [0, 1], η α (x)p αg (x) + g (x) = 0, α = w, o.

(1.15)

1.2. Weak variational formulation. By multiplying (1.2) and (1.3) with a smooth function v, say v ∈ C 1 (Ω × [0, T ]) that vanishes at t = T , applying Green's formula in time and space, and using the boundary and initial conditions (1.4) and (1.5), we formally derive a weak variational formulation

- T 0 Ω ϕ s ∂ t v + T 0 Ω η w ∇ p w • ∇ v = Ω ϕ s 0 v(0) + T 0 Ω f w (s in )q -f w (s)q v, T 0 Ω ϕ s ∂ t v + T 0 Ω η o ∇ p o • ∇ v = - Ω ϕ s 0 v(0) + T 0 Ω f o (s in )q -f o (s)q v.
But in general, the pressures are not sufficiently smooth to make this formulation meaningful and following [7], by using (1.15), it is rewritten in terms of the artificial pressures,

- T 0 Ω ϕ s ∂ t v + T 0 Ω η w ∇(p w + p wg (s)) + ∇ g(s) • ∇ v = Ω ϕ s 0 v(0) + T 0 Ω f w (s in )q -f w (s)q v, T 0 Ω ϕ s ∂ t v + T 0 Ω η o ∇(p o -p og (s)) -∇ g(s) • ∇ v = - Ω ϕ s 0 v(0) + T 0 Ω f o (s in )q -f o (s)q v.
(1.16)

The two formulas coincide when the pressures are slightly more regular. With the above assumptions, problem (1.16) has been analyzed in reference [START_REF] Alt | Nonsteady flow of water and oil through inhomogeneous porous media[END_REF], where it is shown that it has a solution s in L ∞ (Ω×]0, T [) with g(s) in L 2 (0, T ; H 1 (Ω)), p α , α = w, o, in L 2 (Ω×]0, T [) with both p w + p wg (s) and p o -p og (s) in L 2 (0, T ; H 1 (Ω)).

Scheme.

From now on, we assume that Ω is a polygon (d = 2) or Lipschitz polyhedron (d = 3) so it can be entirely meshed.

2.1. Meshes and discretization spaces. The mesh T h is a regular family of simplices K, with a constraint on the angle that will be used to enforce the maximum principle: each angle is not larger than π/2, see [START_REF] Casado-Diaz | Finite elements approximation of second order linear elliptic equations in divergence form with right-hand side in l 1[END_REF]. This is easily constructed in 2D. In 3D, since we only investigate convergence we can embed the domain in a triangulated box. Moreover, since the porosity ϕ is a piecewise constant, to simplify we also assume that the mesh is such that ϕ is a constant per element. The parameter h denotes the mesh size i.e., the maximum diameter of the simplices. On this mesh, we consider the standard finite element space of order one

X h = {v h ∈ C 0 ( Ω) ; ∀K ∈ T h , v h | K ∈ P 1 }. (2.1)
Thus the dimension of X h is the number of nodes, say M , of T h . Let φ i be the Lagrange basis function, that is piecewise linear, and takes the value 1 at node i and the value 0 elsewhere. As usual, the Lagrange interpolation operator I h ∈ L(C 0 ( Ω); X h ) is defined by

∀v ∈ C 0 ( Ω), I h (v) = M i=1 v i φ i , (2.2)
where v i is the value of v at the node of index i. It is easy to see that under the mesh condition, we have ∀K,

K ∇φ i • ∇φ j ≤ 0, ∀i = j. (2.3)
For a given node i, we denote by ∆ i the union of elements sharing the node i and by N (i) the set of indices of all the nodes in ∆ i . In the spirit of [START_REF] Guermond | Invariant domains and first-order continuous finite element approximations for hyperbolic systems[END_REF], we define

c ij = ∆ i ∩∆ j |∇φ i • ∇φ j |, ∀i, j. (2.4)
Recall that the trapezoidal rule on a triangle or a tetrahedron K is

K f ≈ 1 d + 1 |K| d+1 =1 f i ,
where f i is the value of the function f at the th node (vertex), with global number i , of K. For any region O, the notation |O| means the measure (volume) of O.

We define

m i = 1 d + 1 K∈∆ i |K| = 1 d + 1 |∆ i |,
and taking into account the porosity ϕ, we define more generally [START_REF] Alt | Nonsteady flow of water and oil through inhomogeneous porous media[END_REF]. It is well-known that the trapezoidal rule defines a norm on X h , • h , uniformly equivalent to the L 2 norm. Let U h ∈ X h and write

mi (ϕ) = 1 d + 1 K∈∆ i ϕ| K |K|, so that m i = mi
U h = M i=1 U i φ i .
The discrete L 2 norm associated with the trapezoidal rule is

U h h = M i=1 m i |U i | 2 1 2
.

There exist positive constants C and C, independent of h and M , such that

∀U h ∈ X h , C U h 2 L 2 (Ω) ≤ U h 2 h ≤ C U h 2 L 2 (Ω) .
(2.5)

This is also true for other piecewise polynomial functions, but with possibly different constants. The scalar product associated with this norm is denoted by (

•, •) h , ∀U h , V h ∈ X h , (U h , V h ) h = M i=1 m i U i V i . (2.6)
By analogy, we introduce the notation

∀U h , V h ∈ X h , (U h , V h ) ϕ h = M i=1 mi (ϕ)U i V i . (2.7)
The assumptions on the porosity ϕ imply that (2.7) defines a weighted scalar product associated with the weighted norm

• ϕ h , ∀U h ∈ X h , U h ϕ h = (U h , U h ) ϕ h 1 2 ,
that satisfies the analogue of (2.5), with the same constants C and C,

∀U h ∈ X h , C (min Ω ϕ) U h 2 L 2 (Ω) ≤ U h ϕ h 2 ≤ C (max Ω ϕ) U h 2 L 2 (Ω) . (2.8)
2.2. Motivation of the space discretization. While discretizing the time derivative is fairly straightforward, discretizing the space derivatives is more delicate because we need a scheme that is consistent and satisfies the maximum principle for the saturation. For the moment, we freeze the time variable and focus on consistency in space. First, we recall a standard property of functions of X h on meshes satisfying (2.3).

Proof. The first equality is obtained by using (2.3), (2.4) and the fact that

M j=1 φ j = 1,
as in [START_REF] Eymard | Finite volume methods[END_REF] (Section 12.1). For the second part, we use the symmetry of c ij and the antisymmetry of V j -V i to deduce that

- M i=1 U i j =i,j∈N (i) c ij V j -V i = 1 2 M i=1 j =i,j∈N (i) c ij U j -U i V j -V i ,
which is the desired result.

Note that c ij vanishes when j / ∈ N (i). Therefore, when there is no ambiguity it is convenient to write the above double sums on i and j with i and j running from 1 to M .

As an immediate consequence of Proposition 2.1, we have, by taking

V h = U h , ∀U h ∈ X h , ∇ U h L 2 (Ω) = 1 √ 2 M i,j=1 c ij |U j -U i | 2 1 2 .
(2.10) Now, we consider the case of the product of the gradients by a third function. Beforehand, we introduce the following notation: for indices i and j of two neighboring interior nodes, ∆ i ∩ ∆ j in two dimensions is the union of two triangles and in three dimensions the union of a number of tetrahedra bounded by a fixed constant, say L, determined by the regularity of the mesh. We shall use the following notation

c ij,K = K |∇ φ i • ∇ φ j |, w K = 1 |K] K w. (2.11)
Note that

K⊂∆ i ∩∆ j c ij,K = c ij .
(2.12)

Then we have the following proposition:

Proposition 2.2. Let (2.3) hold. With the notation (2.11), the following identity holds for all w in L 1 (Ω):

∀U h , V h ∈ X h , Ω w∇ U h • ∇ V h = - M i=1 U i M j=1 K⊂∆ i ∩∆ j c ij,K w K V j -V i , (2.13)
Proof. It is easy to prove that

Ω w∇ U h • ∇ V h = M i,j=1 d ij U i V j , (2.14)
where

d ij = ∆ i ∩∆ j w(∇ φ i • ∇ φ j ) = Ω w(∇ φ i • ∇ φ j ). (2.15)
Again, we have for any i,

M j=1 d ij = 0, and d ii = - 1≤j≤M,j =i d ij ,
and by substituting this equality into (2.14), we obtain

Ω w ∇ U h • ∇ V h = M i,j=1 U i d ij V j -V i . (2.16)
But, in view of (2.11) and (2.15), and since

∇ φ i • ∇ φ j is a constant in each element K contained in ∆ i ∩ ∆ j , d ij = - K⊂∆ i ∩∆ j c ij,K w K , (2.17)
and (2.13) follows by substituting this equation into (2.16).

Note that d ij = d ji owing to (2.17). The first consequence of Proposition 2.2 is that the right-hand side of (2.13) is a consistent approximation of (w, ∇ u • ∇ v).

Proposition 2.3. Let (2.3) hold, let u and v belong to H 2 (Ω) and w to L ∞ (Ω), and let U h = I h u, V h = I h v be defined by (2.2). Then, there exists a constant C, independent of h, M , u, v, and w, such that

Ω w∇ u • ∇ v + M i,j=1 U i K⊂∆ i ∩∆ j c ij,K w K V j -V i ≤ C h w L ∞ (Ω) u H 2 (Ω) v H 2 (Ω) .
(2.18)

Proof. In view of the identity (2.13), the left-hand side of (2.18) is bounded as follows:

Ω w ∇ u • ∇ v -∇ U h • ∇ V h ≤ w L ∞ (Ω) × ∇(u -U h ) L 2 (Ω) ∇ v L 2 (Ω) + ∇(v -V h ) L 2 (Ω) ∇ U h L 2 (Ω) .
From here, (2.18) is a consequence of standard finite element interpolation error. Now, if w is in W 1,∞ (Ω), then again, standard finite element approximation shows that there exists a constant C, independent of h, K ⊂ ∆ i ∩ ∆ j , and w, such that

w K -w L ∞ (K) ≤ C h |w| W 1,∞ (K) ≤ C h |w| W 1,∞ (Ω) . (2.19)
As a consequence, we will show that in the error formula (2.18), the average w K can be replaced by any value of w in K. Since all K in ∆ i ∩ ∆ j share the edge, say e ij , whose end points are the nodes with indices i and j, then we can pick the value of w at any point, say W i,j , of e ij . At this stage, we choose this value freely, but we prescribe that it be symmetrical with respect to i and j, i.e., W i,j = W j,i .

(2.20)

Then we have the following approximation result.

Theorem 2.4. With the assumption and notation of Proposition 2.3, there exists a constant C, independent of h and M , such that for all u, and v in H 2 (Ω) and

w in W 1,∞ (Ω), Ω w∇ u • ∇ v = - M i,j=1 U i c ij W i,j V j -V i + R, (2.21)
for any arbitrary value W i,j of w in the common edge e ij satisfying (2.20), and the remainder R satisfies

|R| ≤ C h |w| W 1,∞ (Ω) u H 2 (Ω) v H 2 (Ω) . (2.22)
Proof. We infer from (2.12) and (2.13) that

Ω w ∇ U h •∇ V h = - M i,j=1 U i V j -V i K⊂∆ i ∩∆ j c ij,K w K -W i,j - M i,j=1 U i c ij V j -V i W i,j . Let R ij = K⊂∆ i ∩∆ j c ij,K w K -W i,j ,
which is symmetric in i and j by assumption (2.20). As in Proposition 2.1, the symmetry of R ij and the anti-symmetry of V j -V i , imply

- M i,j=1 U i R ij V j -V i ≤ 1 2 M i,j=1 |R ij | U j -U i 2 1 2 M i,j=1 |R ij | V j -V i 2 1 2 . (2.23)
From the non negativity of c ij,K , (2.12), and (2.19), we infer that

|R ij | ≤ K⊂∆ i ∩∆ j c ij,K C h |w| W 1,∞ (Ω) = c ij C h |w| W 1,∞ (Ω) .
Hence, with (2.10) and standard finite element approximation, M i,j=1

U i R ij V j -V i ≤ C h |w| W 1,∞ (Ω) ∇ U h L 2 (Ω) ∇ V h L 2 (Ω) ≤ C h |w| W 1,∞ (Ω) u H 2 (Ω) v H 2 (Ω) .
The result follows by combining this inequality with (2.18).

The above considerations show that

- M i,j=1 U i c ij W i,j V j -V i is a consistent approximation of Ω w∇ u • ∇ v,
for any symmetric choice of W i,j in e ij , the common edge of ∆ i ∩ ∆ j . This will lead to the upwinded space discretization in the next subsection, see also [START_REF] Michel | A finite volume scheme for two-phase incompressible flow in porous media[END_REF]. Furthermore, for all real numbers V i and W i,j satisfying (2.20), 1 ≤ i, j ≤ M , the symmetry of c ij and anti-symmetry of

V j -V i imply M i,j=1 c ij W i,j (V j -V i ) = 0. (2.24)
2.3. Fully discrete scheme. Let τ = T N be the time step, t n = nτ , the discrete times, 0 ≤ n ≤ N . Regarding time, we shall use the standard L 2 projection ρ τ defined on ]t n-1 , t n ], for any function f in L 1 (0, T ), by

ρ τ (f ) n := ρ τ (f )| ]t n-1 ,tn] := 1 τ tn t n-1 f. (2.25)
Regarding space, we shall use a standard element-by-element L 2 projection ρ h as well as a nodal approximation operator r h defined at each node x i for any function g ∈ L 1 (Ω) by

r h (g)(x i ) = 1 |∆ i | ∆ i g, 1 ≤ i ≤ M, (2.26)
and extended to Ω by r h (g) ∈ X h . The operator ρ h is defined for any

f in L 1 (Ω) by ρ h (f )| K = ρ K (f ) where, in any element K, ρ K (f ) = 1 |K| K f. (2.27)
The initial saturation s 0 w is approximated by the operator r h ,

S 0 h = r h (s 0 w ).
(2.28)

The input saturation s in is approximated in space and time by

s in,h,τ = ρ τ (r h (s in )). (2.29)
Clearly, (1.10) implies in space and time 0 ≤ s in,h,τ ≤ 1.

In order to preserve (1.9), the functions q and q are approximated by the functions qh,τ and q h,τ defined with r h and corrected as follows:

qh,τ = ρ τ r h (q) - 1 |Ω| Ω (r h (q) -q) , q h,τ = ρ τ r h (q) - 1 |Ω| Ω (r h (q) -q) . (2.30)
Since qh,τ and q h,τ are piecewise linears in space, they are exactly integrated by the trapezoidal rule and we easily derive from (1.9) and (2.30) that we have for all n,

qn h , 1 h = q n h , 1 h . (2.31)
The set of primary unknowns is the discrete wetting phase saturation and the discrete wetting phase pressure, S n h and P n w,h , defined pointwise at time t n by:

S n h = M i=1 S n,i φ i , P n w,h = M i=1 P n,i w φ i , 1 ≤ n ≤ N.
Then the discrete non-wetting phase pressure P n o,h defined by

P n o,h = M i=1 P n,i o φ i , 1 ≤ n ≤ N,
is a secondary unknown. The upwind scheme we propose for discretizing (1.2)-(1.3) is inspired by the finite volume scheme introduced and analyzed by Eymard al in [START_REF] Eymard | Mathematical study of a petroleum-engineering scheme[END_REF]. For each time step n, 1 ≤ n ≤ N , the lines of the discrete equations are

mi (ϕ) τ (S n,i -S n-1,i ) - M j=1 c ij η w (S n,ij w )(P n,j w -P n,i w ) = m i f w (s n,i in )q n,i -f w (S n,i )q n,i , (2.32)
- mi (ϕ) τ (S n,i -S n-1,i ) - M j=1 c ij η o (S n,ij o )(P n,j o -P n,i o ) = m i f o (s n,i in )q n,i -f o (S n,i )q n,i , (2.33) P n,i o -P n,i w = p c (S n,i ), 1 ≤ i ≤ M, ( 2
τ (S n,M -S n-1,M )- M j=1 c M j η w (S n,M j w )(P n,j w -P n,M w ) -m M f w (s n,M in )q n,M -f w (S n,M )q n,M .
Then, in view of (2.24),

ÃM = M -1 i=1 A i + ÃM = M i=1 mi (ϕ) τ (S n,i -S n-1,i ) - M i=1 m i f w (s n,i in )q n,i -f w (S n,i )q n,i .
By summing in the same fashion the lines of (2.33), we obtain

M i=1 mi (ϕ) τ (S n,i -S n-1,i ) = - M i=1 m i f o (s n,i in )q n,i -f o (S n,i )q n,i .
A combination of these two equations yields

ÃM = - M i=1 m i (f w (s n,i in )+f o (s n,i in ))q n,i -(f w (S n,i )+f o (S n,i ))q n,i = - M i=1 m i (q n,i -q n,i ) = 0,
by virtue of (1.6), the definition (2.25), and (1.9).

3. In (2.32) (respectively, (2.33)), any constant can be added to P w (respectively, P o ), but in view of (2.34), the constant must be the same for both pressures. The last equation (2.35) is added to resolve this constant.

As usual, it is convenient to associate time functions S h,τ , P α,h,τ with the sequences indexed by n. These are piecewise constant in time in ]0, T [, for instance

P α,h,τ (t, x) = P n α,h (x), α = w, o, ∀(t, x) ∈ Ω×]t n-1 , t n ]. (2.38)
In view of the material of the previous subsection, we introduce the following form:

∀W h , U h , V h , Z h ∈ X h , [Z h , W h ; V h , U h ] h = M i,j=1 U i c ij W ij (V j -V i ), (2.39)
where the first argument Z h indicates that the choice of W ij depends on Z h . Such dependence, used for the upwinding, will be specified further on, but it is assumed from now on that W ij satisfies (2.20). Considering (2.24), the form satisfies the following properties,

∀Z h , W h , V h ∈ X h , [Z h , W h ; V h , 1] h = 0, (2.40) ∀Z h , W h , V h ∈ X h , [Z h , W h ; V h , V h ] h = - 1 2 M i,j=1 c ij Wij (V i -V j ) 2 . ( 2
P n o,h in X h , for 1 ≤ n ≤ N , solution of, for all θ h in X h , 1 τ (S n h -S n-1 h , θ h ) ϕ h -P n w,h , I h (η w (S n h )); P n w,h , θ h h = I h (f w (s n in,h ))q n h -I h (f w (S n h ))q n h , θ h h (2.42) - 1 τ (S n h -S n-1 h , θ h ) ϕ h -P n o,h , I h (η o (S n h )); P n o,h , θ h h = I h (f o (s n in,h ))q n h -I h (f o (S n h ))q n h , θ h h (2.43) P n o,h -P n w,h = I h (p c (S n h )), (2.44)
P n w,h , 1 h = 0, (2.45)
where the choice of η w (S n h ) in the left-hand side of (2.42) (respectively, η o (S n h ) in the lefthand side of (2.43)) is given by (2.36) (respectively (2.37)). Strictly speaking, the interpolation operator I h is introduced in (2.42) and (2.43) because the forms are defined for functions of X h , but for the sake of simplicity, since only nodal values are used, it may be dropped further on.

We shall see that under the above basic hypotheses, the discrete problem (2.42)-(2.45) has at least one solution. In the sequel, we shall use the following discrete auxiliary pressures (compare with (1.16)):

U w,h,τ = P w,h,τ + I h (p wg (S h,τ )), U o,h,τ = P o,h,τ -I h (p og (S h,τ )).
(2.46)

The following theorem is the main result of this work:

Theorem 2.6. Under the above basic hypotheses and the additional assumptions (5.42)-(5.44), the discrete solutions converge up to subsequences as follows:

lim (h,τ )→(0,0) S h,τ = s strongly in L 2 (Ω×]0, T [), lim (h,τ )→(0,0)
U w,h,τ = p w + p wg (s), lim

(h,τ )→(0,0) U o,h,τ = p o -p og (s), weakly in L 2 (0, T ; H 1 (Ω)), lim (h,τ )→(0,0) P α,h,τ = p α , weakly in L 2 (Ω×]0, T [), α = w, o,
where p w + p wg , p o -p wg , and s solve the weak formulation (1.16).

The proof of the theorem requires several steps that are covered in the remaining of this work.

3. First a priori bounds. This section is devoted to basic a priori bounds used in proving existence of a discrete solution. Existence is fairly technical and will be postponed till Section 4. The first step is a key bound on the discrete saturation. In a second step, this bound will lead to a pressure estimate and in particular to a bound on the discrete analogue of auxiliary pressures.

3.1. Maximum principle. The scheme (2.32)-(2.35) satisfies the maximum principle property. The proof given below uses a standard argument as in [16, ?].

Theorem 3.1. The following bounds hold:

0 ≤ S h,τ ≤ 1. (3.1)
Proof. As 0 ≤ s 0 w ≤ 1 almost everywhere, by construction (2.28), we immediately have

0 ≤ min Ω s 0 w ≤ S 0 h ≤ max Ω s 0 w ≤ 1.
The proof proceeds by contradiction. Assume that there is an index n ≥ 1 such that

S n-1 h ≤ 1, S n h > 1.
This means there is a node i such that

S n,i = S n h L ∞ (Ω) > 1,
and thus S n,i > S n-1,i .

Dropping the index n in the rest of the proof, (2.32) and (2.33) imply

j =i,j∈N (i) c ij η w (S ij w )(P j w -P i w ) + m i f w (s i in )q i -f w (S i )q i > 0, (3.2)
- j =i,j∈N (i) c ij η o (S ij o )(P j o -P i o ) -m i f o (s i in )q i -f o (S i )q i > 0. (3.3)
We first show that (3.2) holds true with S ij w replaced by S i . Indeed if P i w > P j w , then S ij w = S i . If P i w < P j w , then S ij w = S j , and as η w is increasing and by assumption, S j ≤ S i , η w (S ij w )(P j w -P i w ) ≤ η w (S i )(P j w -P i w ).

Finally, the term vanishes when P i w = P j w . Therefore we have in all cases j =i,j∈N (i)

c ij η w (S i )(P j w -P i w ) + m i f w (s i in )q i -f w (S i )q i > 0. ( 3.4)
A similar argument gives -

c ij η o (S i )(P j o -P i o ) -m i f o (s i in )q i -f o (S i )q i > 0. (3.5)
The substitution of (2.34) into (3.5) yields

- j =i,j∈N (i) c ij η o (S i ) (P j w -P i w ) + (p c (S j ) -p c (S i ) ) -m i f o (s i in )q i -f o (S i )q i > 0. (3.6)
Since p c is decreasing and S i ≥ S j , the second term in the above sum is negative. This implies that

- j =i,j∈N (i) c ij η o (S i )(P j w -P i w ) -m i f o (s i in )q i -f o (S i )q i > 0. (3.7)
The sum on j cancels by multiplying (3.4) by η o (S i ), (3.7) by η w (S i ), and adding the two. The sign is unchanged because either η o (S i ) or η w (S i ) is strictly positive. Hence,

m i η o (S i ) f w (s i in )q i -f w (S i )q i -m i η w (S i ) f o (s i in )q i -f o (S i )q i > 0.
By definition of f w and f o , this reduces to

η o (S i )f w (s i in ) -η w (S i )f o (s i in ) > 0. (3.8)
Now consider the function:

r(s) = η o (s)f w (s i in ) -η w (s)f o (s i in ). (3.9)
It is decreasing and r(s i in ) = 0. Then, since S i > 1 ≥ s i in , see (1.10), we have r(S i ) ≤ r(s i in ) = 0, which contradicts (3.8). The proof of the lower bound in (3.1) follows the same lines. Let us start with a lower bound that removes the degeneracy caused by the mobilities when they multiply the discrete pressures.

Lemma 3.3. Let p wg and p og be defined in (1.12). We have for all n and any i and

j η * (U n,j w -U n,i w ) 2 ≤ η w (S n,ij w )(P n,j w -P n,i w ) 2 + η o (S n,ij o )(P n,j o -P n,i o ) 2 . (3.13)
Proof. To simplify the notation, we drop the superscript n. The second mean formula for integrals gives

p wg (S j ) -p wg (S i ) = S j S i f o (s)p c (s)ds = f o (ξ)(p c (S j ) -p c (S i )), (3.14)
for some ξ between S i and S j . Using (2.34) we write

U j w -U i w = (1 -f o (ξ))(P j w -P i w ) + f o (ξ)(P j o -P i o ) = f w (ξ)(P j w -P i w ) + f o (ξ)(P j o -P i o ).
Therefore since f w + f o = 1, we have

(U j w -U i w ) 2 ≤ η w (ξ) η w (ξ) + η o (ξ) (P j w -P i w ) 2 + η o (ξ) η w (ξ) + η o (ξ) (P j o -P i o ) 2 . ( 3.15)
We now consider six cases.

1) If P i w > P j w and

P i o ≤ P j o , then η w (S ij w ) = η w (S i ) and η o (S ij o ) = η o (S j
) when P i o < P j o ; when P i o = P j o , the value of η o does not matter. From (3.10) we then have S i ≥ S j . Since η w is increasing, η w (ξ) ≤ η w (S i ) and since η o is decreasing, η o (ξ) ≤ η o (S j ). Thus we have

(U j w -U i w ) 2 ≤ η w (S ij w ) η w (ξ) + η o (ξ) (P j w -P i w ) 2 + η o (S ij o ) η w (ξ) + η o (ξ) (P j o -P i o ) 2 ,
and with (1.8)

(U j w -U i w ) 2 ≤ 1 η * η w (S ij w )(P j w -P i w ) 2 + η o (S ij o )(P j o -P i o ) 2 . (3.16)
2) If P i w > P j w and

P i o > P j o , then η w (S ij w ) = η w (S i ) and η o (S ij o ) = η o (S i ). From η o (S i )(p c (S j ) -p c (S i )) = (η o (S i ) + η w (S i )) S j S i
f o (S i )p c (s)ds, and (3.14), we derive

η o (S i )(p c (S j ) -p c (S i ))-(η o (S i ) + η w (S i ))(p wg (S j ) -p wg (S i )) = (η o (S i ) + η w (S i )) S j S i (f o (S i ) -f o (s))p c (s)ds.
As p c and f o are decreasing, the above right-hand side is negative. Hence

η o (S i )(p c (S j ) -p c (S i )) -(η o (S i ) + η w (S i ))(p wg (S j ) -p wg (S i )) ≤ 0. (3.17)
We multiply (3.17) by (P j o -P i o ) + (P j w -P i w ) < 0 and use (2.34),

η o (S i )(p c (S j )-p c (S i ))-(η o (S i )+η w (S i ))(p wg (S j )-p wg (S i )) 2(P j w -P i w ) + p c (S j ) -p c (S i ) ≥ 0.
By expanding and using the next inequality implied by (3.14)

, if f o (ξ) = 0, (p wg (S j ) -p wg (S i ))(p c (S j ) -p c (S i )) ≥ (p wg (S j ) -p wg (S i )) 2 ,
we obtain

η o (S i )(p c (S j ) -p c (S i )) 2 + 2η o (S i )(p c (S j ) -p c (S i ))(P j w -P i w ) ≥ (η o (S i ) + η w (S i ))(p wg (S j ) -p wg (S i )) 2(P j w -P i w ) + p wg (S j ) -p wg (S i ) .
When (η o (S i ) + η w (S i ))(P j w -P i w ) 2 is added to both sides, this becomes

η w (S i )(P j w -P i w ) 2 + η o (S i )(P j o -P i o ) 2 ≥ (η o (S i ) + η w (S i ))(U j w -U i w ) 2 ,
and (1.8) implies the desired result. It remains to consider the case f o (ξ) = 0, i.e., p wg (S j ) = p wg (S i ). If η o (S i ) = 0, then (3.17) yields

p c (S j ) -p c (S i ) ≤ 0 which implies P i o -P j o ≥ P i w -P j w ,
and we deduce immediately

η w (S i )(P j w -P i w ) 2 + η o (S i )(P j o -P i o ) 2 ≥ (η w (S i ) + η o (S i ))(P j w -P i w ) 2 ≥ η * (P j w -P i w ) 2 .
When η o (S i ) = 0, we have trivially

η w (S i )(P j w -P i w ) 2 + η o (S i )(P j o -P i o ) 2 = η w (S i )(P j w -P i w ) 2 ≥ η * (P j w -P i w ) 2 .
3) If P i w ≤ P j w and

P i o > P j o , then η w (S ij w ) = η w (S j ) and η o (S ij o ) = η o (S i
) in the case of a strict inequality; also S i ≤ S j . Then (3.15) and the monotonic properties of η w and η o yield (3.13). If P i w = P j w , then according to (3.11), S i ≤ S j and the same conclusion holds.

4) If P i w ≤ P j w and P i o = P j o , then from (3.12), we have S i ≤ S j and with (3.15)

(U j w -U i w ) 2 ≤ η w (ξ) η w (ξ) + η o (ξ) (P j w -P i w ) 2 ≤ η w (S ij w ) η w (ξ) + η o (ξ) (P j w -P i w ) 2 ,
which is the desired result.

5) Similarly, if P i w = P j w and P i o < P j o , then from (3.11), we have S j ≤ S i and with (3.15)

(U j w -U i w ) 2 ≤ η o (ξ) η w (ξ) + η o (ξ) (P j o -P i o ) 2 ≤ η o (S ij o ) η w (ξ) + η o (ξ) (P j o -P i o ) 2 .
6) If P i w < P j w and P i o < P j o , (3.13) follows from the second case by switching i and j. The pressure bound in the next theorem is the one that arises naturally from the left-hand side of (2.42) and (2.43).

Theorem 3.4. There exists a constant C, independent of h and τ , such that

τ N n=1 M i,j=1 c ij η w (S n,ij w )(P n,i w -P n,j w ) 2 + η o (S n,ij o )(P n,i o -P n,j o ) 2 ≤ C. (3.18)
Proof. We test (2.42) by P n w,h , (2.43) by P n o,h , add the two equations, multiply by τ and sum over n from 1 to N . By using (2.44) and (2.41), we obtain

- N n=1 S n h -S n-1 h , I h p c (S n h ) ϕ h + 1 2 N n=1 τ α=w,o M i,j=1 c ij η α (S n,ij α )(P n,i α -P n,j α ) 2 = N n=1 τ α=w,o f α (s n in,h )q n h -f α (S n h )q n h , P n α,h h . (3.19)
Following [START_REF] Eymard | Mathematical study of a petroleum-engineering scheme[END_REF], the first term in (3. [START_REF] Helmig | Multiphase flow and transport processes in the subsurface: a contribution to the modeling of hydrosystems[END_REF]) is treated with the primitive g c of p c , see (1.11). Indeed, by the mean-value theorem, there exists ξ between S n,i and S n-1,i such that

g c (S n,i ) -g c (S n-1,i ) = -(S n,i -S n-1,i )p c (ξ).
As the function p c is decreasing, then p c (ξ) ≥ p c (S n,i ) when S n,i ≥ S n-1,i and p c (ξ) ≤ p c (S n,i ) when S n,i ≤ S n-1,i . In both cases, we have

g c (S n,i ) -g c (S n-1,i ) ≤ -(S n,i -S n-1,i )p c (S n,i ),
and owing that ϕ is positive and constant in time, (3.19) can be replaced by the inequality

g c (S N h ) -g c (S 0 h ), 1 ϕ h + 1 2 N n=1 τ α=w,o M i,j=1 c ij η α (S n,ij α )(P n,i α -P n,j α ) 2 ≤ N n=1 τ α=w,o f α (s n in,h )q n h -f α (S n h )q n h , P n α,h h . (3.20)
As the first term in the above left-hand side is bounded, owing to the continuity of g c and boundedness of S h,τ , it suffices to handle the right-hand side. Let us drop the superscript n and treat one term in the time sum. Following again [START_REF] Eymard | Mathematical study of a petroleum-engineering scheme[END_REF], in view of Lemma 3.3 we use the auxiliary pressures p wg and p wo , defined in (1.12). Clearly, (1.14) and (2.34) imply

P i w + p wg (S i ) + p og (S i ) + p c (0) = P i o , ∀i. (3.21)
Using this, a generic term, say Y , in the right-hand side of (3.20) can be expressed as

Y = qh -q h , U w,h h + f o (s in,h )q h -f o (S h )q h , p c (0) h + f o (s in,h )q h -f o (S h )q h , p og (S h ) h -f w (s in,h )q h -f w (S h )q h , p wg (S h ) h = T 1 + • • • + T 4 .
We now bound each term T i . For T 1 , (2.31) implies that any constant β can be added to U w,h , in particular β can be chosen so that the sum has zero mean value in Ω. Hence, considering the generalized Poincaré inequality

∀v ∈ H 1 (Ω), v L 2 (Ω) ≤ C Ω v + ∇ v L 2 (Ω) , (3.22)
with a constant C, depending only on the domain Ω, we have

U w,h + β h ≤ C U w,h + β L 2 (Ω) ≤ C ∇ U w,h L 2 (Ω) ,
with another constant C. Then Young's inequality yields

|T 1 | ≤ C 2 2η * qh -q h 2 h + η * 4 ∇ U w,h 2
L 2 (Ω) ,
and with Lemma 3.3, this becomes

|T 1 | ≤ C 2 2η * qh -q h 2 h + 1 4 M i,j=1 c ij η w (S ij )(P j w -P i w ) 2 + η o (S ij )(P j o -P i o ) 2 .
The term T 2 is easily bounded since p c (0) is a number, and so are the terms T 3 and T 4 , in view of the boundedness of the saturation and the continuity of p og and p wg . We thus have

|T 2 + T 3 + T 4 | ≤ C( qh L 1 (Ω) + q h L 1 (Ω) ).
Then substituting these bounds for each n into (3.20), we obtain

1 4 τ N n=1 M i,j=1 c ij η w (S n,ij w )(P n,i w -P n,j w ) 2 + η o (S n,ij o )(P n,i o -P n,j o ) 2 ≤ C qh,τ -q h,τ 2 L 2 (Ω×]0,T [) + qh,τ L 1 (Ω×]0,T [) + q h,τ L 1 (Ω×]0,T [) ,
thus proving (3.18).

By combining Theorem 3.4 with Lemma 3.3, we immediately derive a bound on the discrete auxiliary pressures. The bound (3.23) with α = o follows from the same with α = w, (1.14), and (2.34).

Theorem 3.5. We have for α = w, o

η * ∇ U α,h,τ 2
L 2 (Ω×]0,T [) ≤ C, (3.23)
with the constant C of (3.18).

Existence of numerical solution.

We fix n ≥ 1 and assume there exists a solution

(S n-1 h , P n-1 w,h ) at time t n-1 with 0 ≤ S n-1 h ≤ 1.
We want to show existence of a solution (S n h , P n w,h ) by means of the topological degree [START_REF] Deimling | Nonlinear Functional Analysis[END_REF][START_REF] Dinca | Brouwer degree and applications[END_REF]. Let θ be a constant parameter in [0, 1]. For any continuous function f : [0, 1] → IR, we define the transformed function f :

[0, 1] → IR by ∀s ∈ [0, 1], f (s) = f (ts + (1 -t)θ).
Since θ is fixed, when t = 0, f (s) = f (θ), a constant independent of s. Now, (2.45) implies that any solution P w,h,τ of (2.42)-(2.45) belongs to the following subspace X 0,h of X h ,

X 0,h = {Λ h ∈ X h ; Ω Λ h = 0}. ( 4.1)
This suggests to define the mapping

F : [0, 1] × X h × X 0,h → X h × X 0,h by F(t, ζ, Λ) = (A h , A h + B h ),
where

A h , respectively B h , solves for all Θ h ∈ X h , (A h , Θ h ) = 1 τ (ζ h -S n-1 h , Θ h ) ϕ h -Λ h , I h ( η w (ζ h )); Λ h , Θ h h -I h ( f w (s n in,h ))tq n h -I h ( f w (ζ h ))tq n h , Θ h h , (4.2)
(B h , Θ h ) = - 1 τ (ζ h -S n-1 h , Θ h ) ϕ h -P o,h , I h ( η o (ζ h )); P o,h , Θ h h -I h ( f o (s n in,h ))tq n h -I h ( f o (ζ h ))tq n h , Θ h h , (4.3)
and P o,h is defined by 4.4) determine uniquely A h and B h , and it is easy to check that A h + B h belongs to X 0,h . The mapping t → F(t, ζ h , Λ h ) is continuous. Indeed, since the space has finite dimension, we only need to check continuity of the upwinding. By splitting x into its positive and negative part, x = x + + x -, the upwind term, say η w (ζ ij w )(P j w -P i w ) reads

P o,h = Λ h -I h ( p c (ζ h )). ( 4
η w (ζ ij w )(P j w -P i w ) = η w (tζ i + (1 -t)θ) (P j w -P i w ) -+ η w (tζ j + (1 -t)θ) (P j w -P i w ) + ,
which is continuous with respect to t.

We remark that

F(1, ζ h , Λ h ) = 0 implies that (ζ h , Λ h ) solves (2.42)-(2.45). Conversely, if (ζ h , Λ h ) solves (2.42)-(2.45) then F(1, ζ h , Λ h ) = 0. Thus showing existence of a solution to the problem (2.42)-(2.45) is equivalent to showing existence of a zero of F(1, ζ h , Λ h ).
Before proving existence of a zero, we use the estimates established in the previous section to determine an a priori bound of any zero

(ζ h , Λ h ) of F(1, ζ h , Λ h ).

A priori bounds on (ζ h , Λ h

). In the following we consider t ∈ [0, 1] and

(ζ h , Λ h ) ∈ X h × X 0,h that satisfy F(t, ζ h , Λ h ) = 0. (4.5)
We first show that ζ h satisfies a maximum principle.

Proposition 4.1. The following bounds hold for all (t, ζ h , Λ h ) satisfying (4.5):

0 ≤ ζ h ≤ 1. (4.6)
Proof. Either t ∈]0, 1] or t = 0. The proof for t ∈]0, 1] follows closely the argument used in proving Theorem 3.1 and is left to the reader. For t = 0 we proceed again by contradiction. Assume first that ζ h L ∞ (Ω) > 1, i.e., there is a node i such that

ζ i = ζ h L ∞ (Ω) > 1 ≥ S n-1,i .
As t = 0, (4.5) reduces to

j =i c ij η w (θ)(Λ i -Λ j ) > 0, - j =i c ij η o (θ)(Λ i -Λ j ) > 0, ∀1 ≤ i ≤ M.
Since η o and η w are non-negative functions satisfying (1.8), the inequalities above yield a contradiction. A similar argument is used to show that ζ h ≥ 0.

Next we show the following bound on Λ h .

Proposition 4.2.

There is a constant C such that for all t ∈ [0, 1] we have

η * M i,j=1 c ij Λ j -Λ i + p wg (tζ j + (1 -t)θ) -p wg (tζ i + (1 -t)θ) 2 ≤ C. (4.7)
Proof. The proof follows closely that of Theorem 3.5. First we show there exists a constant

C 1 independent of t such that M i,j=1 c ij η w (tζ ij w + (1 -t)θ)(Λ j -Λ i ) 2 + η o (tζ ij o + (1 -t)θ)(P j o,h -P i o,h ) 2 ≤ C 1 ,
with P o,h defined in (4.4). This bound is obtained via arguments similar to those used in proving Theorem 3.4. The main difference is that the formula is neither summed over n nor multiplied by the time step τ . As a consequence, the constant C 1 includes a term of the form τ -1 g c L ∞ (Ω) arising from the bound of the discrete time derivative. To finish the proof we must show that

η * Λ j -Λ i + p wg (tζ j + (1 -t)θ) -p wg (tζ i + (1 -t)θ) 2 ≤ η w (tζ ij + (1 -t)θ)(Λ j -Λ i ) 2 +η o (tζ ij o + (1 -t)θ)(P j o -P i o ) 2 .
By (1.8), this is trivially satisfied when t = 0. When t ∈]0, 1], the argument is the same as in the proof of Lemma 3.3.

Propositions 4.1 and 4.2 are combined to obtain a bound on

ζ h h + Λ h h . Proposition 4.3. There exists a constant R 1 > 0, independent of t ∈ [0, 1], such that any solution (ζ h , Λ h ) of (4.5) satisfies ζ h h + Λ h h ≤ R 1 . (4.8)
Proof. According to Proposition 4.1, there exists a constant C 1 independent of t such that

ζ h h ≤ C 1 .
To establish a bound on Λ h h , we infer from (1.12) that the function |p wg | is bounded by p c (0) -p c (1) because f o is bounded by one and p c is a decreasing function. Thus (4.7) implies that there exists a constant C 2 independent of t that satisfies

M i,j=1 c ij Λ j -Λ i 2 ≤ C 2 , i.e., ∇Λ h L 2 (Ω) ≤ √ C 2 √ 2 , (4.9)
owing to (2.10). As Λ h ∈ X 0,h , the generalized Poincaré inequality (3.22) shows there exists a constant C 3 independent of t such that

Λ h L 2 (Ω) ≤ C 3 .
Then the equivalence of norm (2.5) yields

Λ h h ≤ C 4 ,
and (4.8) follows by setting R 1 = C 1 + C 4 , a constant independent of t.

Proof of existence.

For any R > 0, let B R denote the ball

B R = {(ζ h , Λ h ) ∈ X h × X 0,h : ζ h h + Λ h h ≤ R}, (4.10)
and let R 0 = R 1 + 1, where R 1 is the constant of (4.8). Since all solutions (ζ h , Λ h ) of (4.5) are in the ball B R 1 , this function has no zero on the boundary ∂B R 0 . Existence of a solution of (2.42)-(2.45) follows from the following result:

Theorem 4.4. The equation F(1, ζ h , Λ h ) = 0 has at least one solution (ζ h , Λ h ) ∈ B R 0 .
Proof. The proof proceeds in two steps. First, we show that the system with t = 0 has a solution:

F(0, ζ h , Λ h ) = 0.
This is a square linear system in finite dimension, so existence is equivalent to uniqueness.

Thus we assume that it has two solutions, and for convenience, we still denote by (ζ h , Λ h ) the difference between the two solutions. The system reads mi

τ ζ i h - j =i,j∈N (i) c ij η w (θ)(Λ j -Λ i ) = 0, 1 ≤ i ≤ M, (4.11)
- mi τ ζ i h - j =i,j∈N (i) c ij η o (θ)(Λ j -Λ i ) = 0, 1 ≤ i ≤ M, (4.12)
i

m i Λ i = 0. (4.13)
We add the first two equations, multiply by Λ i , and sum over i. Then (2.10) and (2.41) imply that Λ h is a constant and finally (4.13) shows that this constant is zero. This yields

ζ h = 0.
Next, we argue on the topological degree. Since the topological degree of a linear map is the sign of its determinant, we have

d(F(0, ζ h , Λ h ), B R 0 , 0) = 0.
We also know that

d(F(t, ζ h , Λ h ), B R 0 , 0) is independent of t since the mapping t → F(t, ζ h , Λ h ) is continuous and for every t ∈ [0, 1], if F(t, ζ h , Λ h ) = 0, then (ζ h , Λ h ) does not belong to ∂B R 0 . Therefore we have d(F(1, ζ h , Λ h ), B R 0 , 0) = d(F(0, ζ h , Λ h ), B R 0 , 0) = 0. This implies that F(1, ζ h , Λ h ) has a zero (ζ h , Λ h ) ∈ B R 0 .
5. Additional pressure estimates. The pressure estimates (3.18) and (3.23) are not sufficient to pass to the limit in the scheme (2.42)-(2.45). These are nonlinear equations and we need strong convergences that do not stem directly from (3.18) and (3.23).

Following [START_REF] Eymard | Mathematical study of a petroleum-engineering scheme[END_REF], we propose to derive a bound for the gradient of g, see (1.13), at s = S h,τ . Then, under suitable assumptions on the behavior of η w , η o , and p c , we shall prove the strong convergence of g(S h,τ ) in L 2 (Ω×]0, T [) and in turn the strong convergence of S h,τ in L 2 (Ω×]0, T [).

Estimating the gradient of g(S h,τ ) is a long and intricate process; it is based on the fact that |g(S n,j ) -g(S n,i )| 2 ≤ C f w (S n,i ) -f w (S n,j ) g(S n,i ) -g(S n,j ) , see (5.59). Therefore, we must derive a bound for the product of the gradients of g and f w . This is split into several steps.

5.1. A preliminary inequality. Our starting step is the following inequality:

Proposition 5.1. There exists a constant C 1 independent of h and τ such that

- N n=1 τ α=o,w P n α,h , η α (S n α,h ); f α (S n h ), P n α,h h = R 1 , (5.1)
where the remainder

R 1 satisfies |R 1 | ≤ C 1 .
Proof. By testing (2.42) with I h f w (S n h ) and (2.43) with I h f o (S n h ), adding the resulting equalities, and multiplying by τ , we obtain

N n=1 S n h -S n-1 h , f w (S n h ) -f o (S n h ) ϕ h - N n=1 τ α=o,w P n α,h , η α (S n α,h ); f α (S n h ), P n α,h h = T 0 qh,τ , α=o,w f α (s in,h,τ )f α (S h,τ ) h -q h,τ , α=o,w (f α (S h,τ )) 2 h ≤ 4 q L 1 (Ω×]0,T [) , (5.2)
in view of (1.6) and (1.9). To control the time difference of S h,τ , we introduce the global flux defined by

∀x ∈ [0, 1], G(x) = x 0 f w (s) -f o (s) ds, (5.3)
and we write

(S n h -S n-1 h ) f w (S n h ) -f o (S n h ) = (S n h -S n-1 h )G (S n h ). But by (1.6), G (x) = 2f w (x) -1 is increasing. Hence, considering that G(S n h ) -G(S n-1 h ) = (S n h -S n-1 h )G (c),
for some c between S n h -S n-1 h , we easily check that

G(S n h ) -G(S n-1 h ) ≤ (S n h -S n-1 h )G (S n h ).
Thus, the properties of ϕ imply

N n=1 S n h -S n-1 h , f w (S n h ) -f o (S n h ) ϕ h ≥ (G(S N h ), 1) ϕ h -(G(S 0 h ), 1) ϕ h .
But the boundedness of S h,τ , the continuity of f α , and the properties of ϕ imply

(G(S N h ), 1) ϕ h -(G(S 0 h ), 1) ϕ h ≤ C ,
with a constant C independent of h and τ . By substituting these inequalities into (5.2) we derive (5.1) with

C 1 = 4 q L 1 (Ω×]0,T [) + C .
5.2. Some discrete total flux inequalities. In this section, it is convenient to work directly on the scheme (2.32)-(2.33). For each index i, the sum of the equations (2.32) and (2.33) give, for 1

≤ i ≤ M and 1 ≤ n ≤ N , - j =i,j∈N (i) c ij η w (S n,ij w )(P n,j w -P n,i w ) + η o (S n,ij o )(P n,j o -P n,i o ) = m i (q n,i -q n,i ).
Following [START_REF] Eymard | Mathematical study of a petroleum-engineering scheme[END_REF], this suggests to define a discrete anti-symmetric upwinded total flux,

F n,ij = -η w (S n,ij w )(P n,j w -P n,i w ) -η o (S n,ij o )(P n,j o -P n,i o ); (5.4) it satisfies j =i,j∈N (i) c ij F n,ij = m i (q n,i -q n,i ). (5.5)
This identity yields a first bound for the discrete total flux.

Proposition 5.2. The discrete total flux F n,ij satisfies the following bounds for α = w, o:

N n=1 τ M i,j=1 f 2 α (S n,i )c ij F n,ij ≤ 2 q L 1 (Ω×]0,T [) . (5.6)
Proof. The statement follows by multiplying (5.5) with τ f 2 α (S n,i ), and summing

N n=1 τ M i,j=1 f 2 α (S n,i )c ij F n,ij = N n=1 τ M i=1 m i f 2 α (S n,i )(q n,i -q n,i ) ≤ 2 q L 1 (Ω×]0,T [) .
To simplify some of the calculations below, it is convenient to drop the time superscript n, when there is no ambiguity, and restore it when needed.

By using the relation (1.7), F i,j can also be written as

F ij = -η w (S ij w ) + η o (S ij o ) (P j w -P i w ) -η o (S ij o )(p c (S j ) -p c (S i )) = -η w (S ij w ) + η o (S ij o ) (P j o -P i o ) + η w (S ij w )(p c (S j ) -p c (S i )).
(5.7)

In order to insert it into (5.1), we bring forward F ij in the expressions for η α (S ij α )(P j α -P i α ), α = w, o. Starting from the identity

η w (S ij w )(P j w -P i w ) = f w (S ij w ) η w (S ij w ) + η o (S ij o ) (P j w -P i w ) + η o (S ij o )(p c (S j ) -p c (S i )) -η o (S ij o )(p c (S j ) -p c (S i )) + η o (S ij w ) -η o (S ij o ) (P j w -P i w ) ,
the expression (5.7) leads to

η w (S ij w )(P j w -P i w ) = f w (S ij w ) -F ij -η o (S ij o )(p c (S j )-p c (S i ))+ η o (S ij w )-η o (S ij o ) (P j w -P i w ) . (5.8) Similarly, η o (S ij o )(P j o -P i o ) = f o (S ij o ) -F ij +η w (S ij w )(p c (S j )-p c (S i ))+ η w (S ij o )-η w (S ij w ) (P j o -P i o ) . (5.9)
We also introduce the anti-symmetric quantities that collect the terms other than F ij in (5.8) and (5.9),

C ij w = η o (S ij o ) p c (S j ) -p c (S i ) -η o (S ij w ) -η o (S ij o ) (P j w -P i w ), (5.10)
C ij o = -η w (S ij w ) p c (S j ) -p c (S i ) -η w (S ij o ) -η w (S ij w ) (P j o -P i o ). (5.11)
With this notation, we have

η α (S ij α )(P j α -P i α ) = f α (S ij α ) -F ij -C ij α , α = w, o.
Thus, the term that is summed over i in (5.1) has the expression

- α=w,o f α (S i ) j =i,j∈N (i) c ij η α (S ij α )(P j α -P i α ) = α=w,o f α (S i ) j =i,j∈N (i) c ij f α (S ij α ) F ij + C ij α .
(5.12)

Now, we reintroduce the superscript n and to simplify, we set

A 1,i,n = α=w,o f α (S n,i ) M j=1 c ij f α (S n,ij α )F n,ij , (5.13)
A α,i,n = f α (S n,i ) M j=1 c ij f α (S n,ij α )C n,ij α . (5.14)
With this notation, our next proposition is derived by substituting (5.12)-(5.14) into (5.1).

Proposition 5.3. We have, with the remainder R 1 of (5.1),

N n=1 τ M i=1 A 1,i,n + N n=1 τ M i=1 α=w,o A α,i,n = R 1 .
(5.15)

We must transform suitably each term in this sum to bring forward g. Let us start with the first term of (5.15), i.e., the combination of the discrete total flux.

Combination of the discrete total flux.

To simplify, let A 1 denote the first term,

A 1 = N n=1 τ M i,j=1 α=w,o f α (S n,i )c ij f α (S n,ij α )F n,ij .
Inspired by (5.6), we introduce the difference

A 2 = A 1 - N n=1 τ M i,j=1 f 2 w (S n,i ) + f 2 o (S n,i ) c ij F n,ij .
Clearly, A 2 collects the discrepancies arising from the upwinding,

A 2 = N n=1 τ M i,j=1 α=w,o f α (S n,i )c ij f α (S n,ij α ) -f α (S n,i ) F n,ij .
(5.16)

As (5.6) yields

A 1 = A 2 + R 2 , with |R 2 | ≤ 4 q L 1 (Ω×]0,T [) , (5.17)
a bound for A 1 stems from a bound for A 2 . To this end, in view of (5.16), it is useful to consider the four subsets of indices j ∈ N (i), j = i, union and intersection:

N w (i) = {j ∈ N (i) ; P n,j w > P n,i w }, N o (i) = {j ∈ N (i) ; P n,j o > P n,i o } N w,S (i) = {j ∈ N (i), j = i ; P n,j w = P n,i w , S n,j ≥ S n,i }, N o,S (i) = {j ∈ N (i), j = i ; P n,j o = P n,i o , S n,j ≤ S n,i }, UN (i) = N w (i) ∪ N o (i) ∪ N w,S (i) ∪ N o,S (i), N F (i) = {j ∈ N (i) ; P n,i
w > P n,j w and P n,i o > P n,j o }.

(5.18)

Strictly speaking, these subsets should we written with the superscript n, but we omit it for the sake of simplicity. Then we have the following bound for A 2 :

Proposition 5.4. There exists a constant C 2 , independent of h and τ , such that

A 2 = - 1 2 N n=1 τ M i=1 j∈U N (i) c ij f w (S n,j ) -f w (S n,i ) 2 F n,ij + R 3 , (5.19)
where the remainder R 3 satisfies

|R 3 | ≤ C 2 = 2 q L 1 (Ω×]0,T [) .
Proof. Let us drop the superscript n. By definition, f w (S ij w ) -f w (S i ) = 0 when P i w > P j w and when P i w = P j w and S i > S j . Similarly, f o (S ij o ) -f o (S i ) = 0 when P i o > P j o and when P i o = P j o and S n,i < S j . Therefore, the nth term in A 2 , say a 2 , reduces to

a 2 = M i=1 α=w,o f α (S i ) j∈Nα(i)∪N α,S (i) c ij f α (S j ) -f α (S i ) F ij .
By expanding the products, this can be written

a 2 = - 1 2 M i=1 α=w,o j∈Nα(i)∪N α,S (i) c ij f 2 α (S i ) -f 2 α (S j ) + (f α (S i ) -f α (S j )) 2 F ij .
(5.20)

Since c ij vanishes when j is not a neighbor of i, we have, by interchanging i and j and using the anti-symmetry of F ij and the symmetry of c ij ,

- M i=1 j∈Nw(i) c ij f 2 w (S j )F ij = M i=1,j=1,P j w <P i w c ij f 2 w (S i )F ij . (5.21)
Similarly,

- M i=1 j∈N w,S (i) c ij f 2 w (S j )F ij = M i=1,j=1,P i w =P j w ,S i ≥S j c ij f 2 w (S i )F ij . (5.22) Hence - 1 2 M i=1 j∈Nw(i) c ij f 2 w (S i ) -f 2 w (S j ))F ij = - 1 2 M i=1,j=1,P i w =P j w c ij f 2 w (S i )F ij ,
and

- 1 2 M i=1 j∈N w,S (i) c ij f 2 w (S i ) -f 2 w (S j ) F ij = - 1 2 M i=1,j=1,P i w =P j w c ij f 2 w (S i )F ij ,
because there is no contribution from the indices i, j such that P i w = P j w , S i = S j since in this case the factor F ij = 0. The same is true for the non-wetting phase. Thus

- 1 2 α=w,o M i=1 j∈Nα(i)∪N α,S (i) c ij f 2 α (S i ) -f 2 α (S j ) F ij = - 1 2 α=w,o M i=1,j=1 c ij f 2 α (S i )F ij .
By comparing with (5.6), we see that

1 2 N n=1 τ M i=1 α=w,o j∈Nα(i)∪N α,S (i) c ij f 2 α (S n,i ) -f 2 α (S n,j ) F n,ij ≤ 2 q L 1 (Ω×]0,T [) .
( 5.23) This and the equality

f o (S n,j ) -f o (S n,i ) 2 = f w (S n,j ) -f w (S n,i ) 2 ,
readily imply (5.19).

Now, we set

A ij = c ij f w (S j ) -f w (S i ) 2 F ij , a 3 = - 1 2 M i=1 j∈U N (i) A ij .
(5.24)

The next proposition simplifies the expression for a 3 .

Proposition 5.5. We have

a 3 = M i=1 j∈N F (i) c ij f w (S j ) -f w (S i ) 2 F ij .
(5.25)

Proof. By expanding the indices in the set UN (i), interchanging the indices i and j, and using the anti-symmetry of A ij , we derive

a 3 = 1 2 α=w,o P i α >P j α A ij + P i w =P j w ,S j ≤S i A ij + P i o =P j o ,S i ≤S j A ij .
Now, we split the first two sums above as follows:

α=w,o P i α >P j α A ij = 2 P i w >P j w ,P i o >P j o A ij + P i w >P j w ,P i o ≤P j o A ij + P i o >P j o ,P i w ≤P j w A ij .
This leads to

a 3 = j∈N F (i) A ij + 1 2 P i w >P j w ,P i o ≤P j o A ij + P i o >P j o ,P i w ≤P j w A ij + P i w =P j w ,S j ≤S i A ij + P i o =P j o ,S i ≤S j A ij .
The anti-symmetry of A ij gives

P i w >P j w ,P i o ≤P j o A ij = - P j w >P i w ,P j o <P i o A ij - P j w >P i w ,P j o =P i o A ij .
By substituting and applying twice again the anti-symmetry of A ij , we derive

a 3 = j∈N F (i) A ij + 1 2 P i o >P j o ,P i w =P j w A ij + P i w >P j w ,P i o =P j o A ij - P i w =P j w ,S i ≤S j A ij - P i o =P j o ,S j ≤S i A ij .
(5.26) Note that

P i o >P j o ,P i w =P j w A ij = P i o ≥P j o ,P i w =P j w A ij ,
since the additional term is zero. Therefore, in view of first (3.11) and next (3.12),

P i o >P j o ,P i w =P j w A ij = P i w =P j w ,S i ≤S j A ij , P i w >P j w ,P i o =P j o A ij = P i o =P j o ,S i ≥S j A ij .
Thus all terms multiplying 1 2 in (5.26) are cancelled and we recover (5.25). By applying (5.17) and Propositions 5.4 and 5.5, A 1 has the following expression: Proposition 5.6. We have

N n=1 τ M i,j=1 α=w,o f α (S n,i )c ij f α (S n,ij α )F n,ij = N n=1 τ M i=1 j∈N F (i) c ij f w (S n,j ) -f w (S n,i ) 2 F n,ij + R 4 , (5.27)
where

|R 4 | ≤ 6 q L 1 (Ω×]0,T [) .
(5.28)

This settles the contribution of the first term of (5.15); the second terms are handled in the next subsection.

5.4. Terms involving the capillary pressure and mobility. These are the terms A α,i,n defined in (5.14). By virtue of the anti-symmetry of C ij , we can write for α = w, o

M i,j=1 f α (S i )c ij f α (S ij α ) C ij α = - 1 2 M i,j=1 f α (S j ) -f α (S i ) c ij f α (S ij α ) C ij α .
(5.29)

Owing to (1.6), the term with α = o in the right-hand side is

1 2 M i,j=1 f w (S j )-f w (S i ) c ij f o (S ij o ) C ij o . Therefore, α=w,o M i=1 A α,i,n = 1 2 M i,j=1 c ij f w (S n,j ) -f w (S n,i ) -f w (S n,ij w ) C n,ij w + f o (S n,ij o ) C n,ij o . (5.30)
Let K ij denote the symmetric term

K ij := c ij f w (S n,j ) -f w (S n,i ) -f w (S n,ij w ) C n,ij w + f o (S n,ij o ) C n,ij o ;
by virtue of this symmetry, we have

α=w,o M i=1 A α,i,n = P n,i w >P n,j w K n,ij + 1 2 P n,i w =P n,j w K n,ij .
(5.31) 5.5. Combining all terms. By substituting (5.27) and (5.31) into (5.15), we obtain the next lemma.

Lemma 5.7. We have

- N n=1 τ α=o,w P n α,h , η α (S n α,h ); f α (S n h ), P n α,h h = N n=1 τ M i=1 j∈N F (i) c ij f w (S n,j ) -f w (S n,i ) 2 F n,ij - P n,i w >P n,j w c ij f w (S n,j ) -f w (S n,i ) f w (S n,ij w ) C n,ij w -f o (S n,ij o ) C n,ij o - 1 2 P n,i w =P n,j w c ij f w (S n,j ) -f w (S n,i ) f w (S n,ij w ) C n,ij w -f o (S n,ij o ) C n,ij o + R 4 .
(5.32)

with R 4 bounded by (5.28).

Thus, to bring forward g, we must suitably combine the terms of the above sum over i, and this is done by examining all pairs of indices (i, j) involved in (5.32), i.e., the pairs of indices in the following sets: (i) P i w > P j w and P i o > P j o , (ii) P i w > P j w and P i o < P j o , (iii) P i w > P j w and P i o = P j o , (iv) P i w = P j w and P i o > P j o , (v) P i w = P j w and P i o < P j o . Note that the sixth case that would be P i w = P j w and P i o = P j o brings no information because it implies that S i = S j .

For the argument below, we shall use the following intermediate result.

Proposition 5.8. For each indices i and j, there exist (non unique) points α and α between S i and S j such that

g(S j ) -g(S i ) = -η o (α)f w (α) p c (S j ) -p c (S i ) = -η w (α )f o (α ) p c (S j ) -p c (S i ) . (5.33)
Proof. By the definition (1.13),

g(S j ) -g(S i ) = - S j S i η o (x)f w (x)p c (x) dx = - S j S i η w (x)f o (x)p c (x) dx.
(5.34)

Since the functions η o f w and η w f o are continuous and do not change sign between S i and S j , (5.33) follows from the second mean formula for integrals.

To simplify, the superscript n is dropped.

5.5.1. The case P i w > P j w and P i o > P j o . The following holds: Proposition 5.9. Let P i w > P j w and P i o > P j o ; then the factor of τ in (5.32) satisfies

c ij f w (S j ) -f w (S i ) f w (S j ) -f w (S i ) F ij -f w (S i )C ij w -f o (S i )C ij o ≥ c ij f w (S j ) -f w (S i ) g(S j ) -g(S i
) .

(5.35)

Proof. Let E ij denote the left-hand side of (5.35). In the case P i w > P j w and P i o > P j o , an expansion of F ij , and

C ij α yields E ij = = -c ij f w (S j ) -f w (S i ) 2 η w (S i ) + η o (S i ) P j w -P i w -c ij f w (S j ) -f w (S i ) p c (S j ) -p c (S i ) η o (S i ) f w (S j ) -f w (S i ) + 2 η w (S i )η o (S i ) η w (S i ) + η o (S i )
.

As P i w > P j w , the first line in the above right-hand side is nonnegative and hence

E ij ≥ -c ij f w (S j ) -f w (S i ) p c (S j ) -p c (S i ) η o (S i ) f w (S j ) + f w (S i ) .
Now, either S i ≤ S j or S i > S j . If S i ≤ S j , then f w (S j ) ≥ f w (S i ) and p c (S j ) ≤ p c (S i ) because f w is increasing and p c is decreasing. This implies in particular that

E ij ≥ -c ij f w (S j ) -f w (S i ) p c (S j ) -p c (S i ) η o (S i )f w (S j ) ≥ -c ij f w (S j ) -f w (S i ) p c (S j ) -p c (S i ) η o (x)f w (x),
for any x between S i and S j . Then (5.35) follows from the first part of (5.33).

If S i > S j , then f w (S j ) -f w (S i ) ≤ 0, p c (S j ) -p c (S i ) ≥ 0, and we infer from (1.7) that E ij reads

E ij =c ij f w (S j ) -f w (S i ) η w (S i ) + η o (S i ) f w (S j ) -f w (S i ) P i o -P j o + f w (S j ) -f w (S i ) η w (S i ) + η o (S i ) p c (S j ) -p c (S i ) + p c (S i ) -p c (S j ) η o (S i ) f w (S j ) -f w (S i ) + 2 η w (S i )η o (S i ) η w (S i ) + η o (S i )
.

Since P i o -P j o > 0, the first line in the above right-hand side is nonnegative, and thus

E ij ≥c ij f w (S i ) -f w (S j ) p c (S j ) -p c (S i ) f w (S i ) -f w (S j ) η w (S i ) + η o (S i ) + f w (S j ) -f w (S i ) η o (S i ) + 2 η w (S i )η o (S i ) η w (S i ) + η o (S i ) ,
which reduces to

E ij ≥ c ij f w (S i ) -f w (S j ) p c (S j ) -p c (S i ) η w (S i ) f o (S j ) + f o (S i ) .
This leads for instance to

E ij ≥ -c ij f w (S j ) -f w (S i ) p c (S j ) -p c (S i ) η w (S i )f o (S j ) ≥ -c ij f w (S j ) -f w (S i ) p c (S j ) -p c (S i ) η w (x)f o (x),
for any x between S i and S j . Then (5.35) follows from the second part of (5.33).

5.5.2. The case P i w > P j w and P i o < P j o . We have Proposition 5.10. Let P i w > P j w and P i o < P j o ; then the factor of τ in (5.32) satisfies

-c ij f w (S j ) -f w (S i ) f w (S i )C ij w -f o (S j )C ij o ≥ c ij f w (S j ) -f w (S i ) g(S j ) -g(S i
) .

(5.36)

Proof. Let E ij denote the left-hand side of (5.36). We have P j w -P i w < 0 and P j o -P i o > 0. Then S ij w = S i and S ij o = S j ; also S j ≤ S i which implies that η w (S j ) ≤ η w (S i ) and η o (S j ) ≥ η o (S i ). The expression for C ij α becomes (see (5.10) and (5.11))

C ij w = η o (S j ) p c (S j ) -p c (S i ) -η o (S i ) -η o (S j ) (P j w -P i w ), C ij o = -η w (S i ) p c (S j ) -p c (S i ) -η w (S j ) -η w (S i ) (P j o -P i o )
.

Hence

E ij = -c ij f w (S j ) -f w (S i ) f w (S i ) η o (S j ) p c (S j ) -p c (S i ) + η o (S j ) -η o (S i ) (P j w -P i w ) + f o (S j ) η w (S i ) p c (S j ) -p c (S i ) + η w (S j ) -η w (S i ) (P j o -P i o )) = -c ij f w (S j ) -f w (S i ) (P j o -P i o ) f w (S i ) η o (S j ) -η o (S i ) + f o (S j ) η w (S j ) -η w (S i ) + p c (S j ) -p c (S i ) f w (S i )η o (S i ) + f o (S j )η w (S i ) = -c ij f w (S j ) -f w (S i ) (P j o -P i o ) f w (S i ) η o (S j ) -η o (S i ) + f o (S j ) η w (S j ) -η w (S i ) + p c (S j ) -p c (S i ) η w (S i ) f o (S i ) + f o (S j ) .
It follows from the above considerations that

-c ij f w (S j ) -f w (S i ) p c (S j ) -p c (S i ) η w (S i ) f o (S i ) + f o (S j ) ≥ -c ij f w (S j ) -f w (S i ) p c (S j ) -p c (S i ) η w (S i )f o (S j ) ≥ -c ij f w (S j ) -f w (S i ) p c (S j ) -p c (S i ) η w (x)f o (x) ,
for any x between S i and S j . Now, the sign of the factor

f w (S i ) η o (S j ) -η o (S i ) + f o (S j ) η w (S j ) -η w (S i )
is not clear. If it is nonnegative, then the whole term

-c ij f w (S j ) -f w (S i ) (P j o -P i o ) f w (S i ) η o (S j ) -η o (S i ) + f o (S j ) η w (S j ) -η w (S i )
is also nonnegative,

E ij ≥ -c ij f w (S j ) -f w (S i ) p c (S j ) -p c (S i ) η w (x)f o (x) .
and (5.36) follows from (5.33). If

f w (S i ) η o (S j ) -η o (S i ) + f o (S j ) η w (S j ) -η w (S i ) < 0,
then we rewrite E ij in terms of P w ,

E ij = -c ij f w (S j ) -f w (S i ) (P j w -P i w ) f w (S i ) η o (S j ) -η o (S i ) + f o (S j ) η w (S j ) -η w (S i ) + p c (S j ) -p c (S i ) f w (S j )η o (S j ) + η w (S j )f o (S j ) .
Since the first line is now nonnegative, we infer

E ij ≥ -c ij f w (S j ) -f w (S i ) p c (S j ) -p c (S i ) η o (S j ) f w (S j ) + f w (S i ) ≥ -c ij f w (S j ) -f w (S i ) p c (S j ) -p c (S i ) η o (x)f w (x),
again for any x between S i and S j , and the result follows from (5.33).

5.5.3. The case P i w > P j w and P i o = P j o . In this case, p c (S j ) -p c (S i ) > 0, S j ≤ S i , and we have the following result: Proposition 5.11. Let P i w > P j w and P i o = P j o ; then the factor of τ in (5.32) satisfies

-c ij f w (S j ) -f w (S i ) f w (S i )C ij w -f o (S j )C ij o ≥ c ij f w (S j ) -f w (S i ) g(S j ) -g(S i
) .

(5.37)

Proof. Let E ij denote the left-hand side of (5.37). We have

C ij w = η o (S i ) p c (S j ) -p c (S i ) and C ij o = -η w (S i ) p c (S j ) -p c (S i ) .
Hence

E ij = -c ij f w (S j ) -f w (S i ) f w (S i )η o (S i ) + f o (S j )η w (S i ) p c (S j ) -p c (S i ) = -c ij f w (S j ) -f w (S i ) p c (S j ) -p c (S i ) η o (S i )f w (S i ) + η w (S i )f o (S j ) .
From here, (5.37) is derived as in the end of the proof of Proposition 5.10.

5.5.4. The case P i w = P j w and P i o > P j o . In this case, p c (S j ) ≤ p c (S i ) and S j ≥ S i . We have the following result: Proposition 5.12. Let P i w = P j w and P i o > P j o ; then the factor of τ in (5.32) satisfies

- 1 2 c ij f w (S j ) -f w (S i ) f w (S j )C ij w -f o (S i )C ij o ≥ 1 2 c ij f w (S j ) -f w (S i ) g(S j ) -g(S i
) .

(5.38)

Proof. In this case, C ij w = η o (S i ) p c (S j ) -p c (S i ) and C ij o = -η w (S i ) p c (S j ) -p c (S i ) . Then the left-hand side E ij of (5.38) is

E ij = - 1 2 c ij f w (S j ) -f w (S i ) p c (S j ) -p c (S i ) f w (S j )η o (S i ) + f o (S i )η w (S i ) ≥ - 1 2 c ij f w (S j ) -f w (S i ) p c (S j ) -p c (S i ) f w (S j )η o (S i ),
and the proof of (5.38) proceeds as above.

5.5.5. The case P i w = P j w and P i o < P j o . In this case, p c (S j ) ≥ p c (S i ) and hence S i ≥ S j . We have the following result: Proposition 5.13. Let P i w = P j w and P i o < P j o ; then the factor of τ in (5.32) satisfies

- 1 2 c ij f w (S j ) -f w (S i ) f w (S i )C ij w -f o (S j )C ij o ≥ 1 2 c ij f w (S j ) -f w (S i ) g(S j ) -g(S i ) .
(5.39)

Proof. In this case, C ij w = η o (S j ) p c (S j ) -p c (S i ) and C ij o = -η w (S j ) p c (S j ) -p c (S i ) . Then the left-hand side E ij of (5.38) is

E ij = - 1 2 c ij f w (S j ) -f w (S i ) p c (S j ) -p c (S i ) f w (S i )η o (S j ) + f o (S j )η w (S j ) ≥ - 1 2 c ij f w (S j ) -f w (S i ) p c (S j ) -p c (S i ) η o (S j )f w (S i ),
and the proof of (5.39) ends as above.

5.6. Auxiliary bound for the gradient of g. The following theorem is the first outcome of this section.

Theorem 5.14. There exists a constant C, independent of h and τ , such that

T 0 Ω ∇(I h (f α (S h,τ ))) • ∇(I h (g(S h,τ ))) ≤ C, α = w, o.
(5.40)

Proof. Owing to (1.6), it suffices to prove (5.40) when α = w. By applying Propositions 5.9-5.13 to Lemma 5.7 and combining with Proposition 5.1, we readily derive that

N n=1 τ M i=1 j∈N (i),P n,i w ≥P n,j w c ij f w (S n,j ) -f w (S n,i ) g(S n,j ) -g(S n,i ) ≤ C, (5.41)
with a constant C independent of h and τ . Therefore, (5.40) will follow if we bound the summand for all j such that P n,i w < P n,j w . But the symmetry of the summand implies that

M i=1 j∈N (i),P n,i w <P n,j w c ij f w (S n,j ) -f w (S n,i ) g(S n,j ) -g(S n,i ) = M i=1 j∈N (i),P n,i w >P n,j w c ij f w (S n,j ) -f w (S n,i ) g(S n,j ) -g(S n,i ) .
Hence

Ω ∇(I h (f w (S n h )))•∇(I h (g(S n h ))) = 2 M i=1 j∈N (i),P n,i w >P n,j w c ij f w (S n,j ) -f w (S n,i ) g(S n,j ) -g(S n,i ) + M i=1 j∈N (i),P n,i w =P n,j w c ij f w (S n,j ) -f w (S n,i ) g(S n,j ) -g(S n,i ) ,
and (5.40), with another constant C, follows by substituting this equality into (5.41).

5.7.

Bound for the gradient of g. In order to deduce from (5.40) a direct bound for the gradient of g, we need to sharpen the assumptions on the mobility.

Assumptions on the mobility.

Here we assume that η w (0) = η o (1) = 0, and η α ∈ W 1,∞ (0, 1), α = w, o. Furthermore, we assume that, for all x ∈]0, 1[,

α w x θw-1 ≤ η w (x) ≤ 1 α w x θw-1 , θ w ≥ 1, 0 < α w ≤ 1, (5.42) α o (1 -x) θo-1 ≤ -η o (x) ≤ 1 α o (1 -x) θo-1 , θ o ≥ 1, 0 < α o ≤ 1, (5.43) 1 α 3 x β 3 -1 (1 -x) β 4 -1 ≥ -p c (x) ≥ α 3 x β 3 -1 (1 -x) β 4 -1 , β 3 , β 4 > 0, 0 < α 3 ≤ 1.
(5.44) From (5.42) and (5.43), we deduce respectively, for all x ∈]0, 1[,

α w θ w x θw ≤ η w (x) ≤ 1 α w θ w
x θw , (5.45)

α o θ o (1 -x) θo ≤ η o (x) ≤ 1 α o θ o (1 -x) θo . (5.46)
The sum of these two inequalities reads for all x ∈]0, 1[,

α w θ w x θw + α o θ o (1 -x) θo ≤ η w (x) + η o (x) ≤ 1 α w θ w x θw + 1 α o θ o (1 -x) θo .
Let denote the lower bound in this inequality. It is easy to check that is a nonnegative continuous function of x on [0, 1], hence uniformly continuous. Therefore, it is bounded and as it does not vanish in this interval, it is bounded away from zero. Thus there exists a positive constant C min such that

∀x ∈ [0, 1], C min ≤ (x) := α w θ w x θw + α o θ o (1 -x) θo ≤ C max , (5.47)
where

C max = max x∈[0,1] 1 α w θ w x θw + 1 α o θ o (1 -x) θo .
(5.48) 5.7.2. Properties of the derivatives of f w and g. By definition, we have

g (x) = - η w (x)η o (x) η w (x) + η o (x) p c (x),
which is positive in ]0, 1[. Considering (1.8), (5.45), and (5.46), we infer

g (x) ≤ 1 η * α 3 1 α w θ w 1 α o θ o x θw-1+β 3 (1 -x) θo-1+β 4 , (5.49)
thus implying that g is a bounded function, i.e., g is Lipschitz continuous. Note that the Lipschitz constant L of g is bounded by

L ≤ 1 αα 3 1 α w θ w 1 α o θ o max x∈[0,1] x θw-1+β 3 (1 -x) θo-1+β 4 .
(5.50)

On the other hand, (5.44)-(5.47) yield for all x ∈]0, 1[,

g (x) ≥ α 3 C max α w θ w α o θ o x θw-1+β 3 (1 -x) θo-1+β 4 > 0.
(5.51)

Thus g ∈ W 1,∞ (0, 1) is a strictly monotonic increasing function on [0, 1] with range [0, β] for some β > 0, hence invertible with inverse g -1 ∈ W 1,∞ (0, β). Now, we turn to f w . By definition, we have

f w (x) = 1 (η w (x) + η o (x)) 2 η o (x)η w (x) -η w (x)η o (x) .
(5.52)

The inequalities (5.42)-(5.47) imply that

f w (x) ≥ 1 C 2 max α o α w 1 θ o x θw-1 (1 -x) θo + 1 θ w x θw (1 -x) θo-1 . Thus, ∀x ∈ [0, 3 4 ], f w (x) ≥ α o α w C 2 max θ o 1 4 θo x θw-1 , (5.53) and ∀x ∈ [ 1 4 , 1], f w (x) ≥ α o α w C 2 max θ w 1 4 θw (1 -x) θo-1 .
(5.54)

Let us use these results to compare g and f w . It follows from (5.49) that

∀x ∈ [0, 3 4 ], g (x) ≤ 1 αα 3 1 α w θ w 1 α o θ o θ o C 2 max α o α w α o α w C 2 max θ o x θw-1 ,
and by setting

C 1 = 1 αα 3 1 α w θ w 1 α o θ o 4 θo θ o C 2 max α o α w
and comparing with (5.53), we obtain

∀x ∈ [0, 3 4 ], g (x) ≤ C 1 f w (x). (5.55) Similarly, ∀x ∈ [ 1 4 , 1], g (x) ≤ 1 αα 3 1 α w θ w 1 α o θ o θ w C 2 max α o α w α o α w C 2 max θ w (1 -x) θo-1 ,
so that, by setting

C 2 = 1 αα 3 1 α w θ w 1 α o θ o 4 θw θ w C 2 max
α o α w and comparing with (5.54), we deduce

∀x ∈ [ 1 4 , 1], g (x) ≤ C 2 f w (x).
(5.56)

This leads to the desired relation between the derivative of f w and g:

∀x ∈ [0, 1], g (x) ≤ Cf w (x), (5.57)
where

C = max(C 1 , C 2 ).
The main result of this section follows by combining (5.57) with (5.40).

Theorem 5.15. Under the assumptions (5.42)-(5.44) on the derivatives of the mobilities and capillary pressure, there exists a constant C, independent of h and τ , such that

∇(I h (g(S h,τ ))) L 2 (Ω×]0,T [) ≤ C.
(5.58)

Proof. Let (i, j) be any pair of indices. If S n,i ≤ S n,j , then by (5.57),

f w (S n,j ) -f w (S n,i ) = S n,j S n,i f w (x) dx ≥ C S n,j
S n,i g (x) dx = C(g(S n,j ) -g(S n,i )).

As g is increasing, we have g(S n,j ) -g(S n,i ) ≥ 0. Therefore

f w (S n,j ) -f w (S n,i ) g(S n,j ) -g(S n,i ) ≥ C|g(S n,j ) -g(S n,i )| 2 .
(5.59) By changing both signs, the same result holds when S n,j < S n,i . Then (5.58) follows from (5.40).

6.

Convergence. The interpolants of p αg (S h,τ ), g(S h,τ ), and p c (S h,τ ) play an important part in this work, see Theorems 3.5 and 5.15, and (2.44). Therefore, we begin by studying convergence properties first of I h (g(S h,τ )) and I h (p αg (S h,τ )), α = w, o, and next I h (p c (S h,τ )). Some results will stem from an interesting relation between differences in values of S h,τ and g(S h,τ ).

Properties of I h (g(S h,τ

)) and I h (p αg (S h,τ )), α = w, o .

Convergence properties of I h (g(S h,τ )). Let K be an element of T h with vertices a

i , 1 ≤ i ≤ d + 1 (local numbers); then K |g(S h,τ )(t n )| 2 = K |g( d+1 i=1 S n,i φ i )| 2 .
As 0 ≤ S n,i , φ i ≤ 1 and g is increasing, we have

0 ≤ d+1 i=1 S n,i φ i ≤ d+1 i=1 S n,i , g( d+1 i=1 S n,i φ i ) ≤ d+1 i=1 g(S n,i ). Hence K |g(S h,τ )(t n )| 2 ≤ (d + 1)|K| d+1 i=1 (g(S n,i )) 2 .
As a consequence, there exist constants C, D, E, independent of n, h and τ , such that

g(S h,τ )(t n ) L 2 (Ω) ≤ C I h (g(S h,τ )(t n )) h ≤ D I h (g(S h,τ )(t n )) ϕ h ≤ E I h (g(S h,τ )(t n )) L 2 (Ω) , (6.1
) owing to (2.5) and (2.8). These inequalities carry over to the norm in L 2 (Ω×]0, T [). Now, let us prove the following convergence property of I h (g(S h,τ )). Lemma 6.1. Under the assumptions of Theorem 5.15, we have

lim (h,τ )→(0,0) g(S h,τ ) -I h (g(S h,τ )) L 2 (Ω×]0,T [) = 0. (6.2)
Proof. For any x in any element K of T h , we have

I h (g(S h,τ ))(x, t n ) -g(S h,τ )(x, t n ) = d+1 i=1 g(S n,i )φ i (x) -g d+1 i=1 S n,i φ i (x) .
As S h,τ is a polynomial of degree one in K, it attains its maximum and its minimum in space at vertices of K, say g(S n, ) and g(S n,r ) are its maximum and minimum respectively. Thus, recalling that g is a nonnegative monotonically increasing function,

d+1 i=1 g(S n,i )φ i (x) ≤ g(S n, ), g d+1 i=1 S n,i φ i (x) ≥ g(S n,r ).
Hence

I h (g(S h,τ )) -g(S h,τ ) 2 L 2 (Ω×]0,T [) ≤ N n=1 τ K∈T h |K||g(S n, ) -g(S n,r )| 2 . (6.3)
For any node i, let

κ i = Max |K|,
where the maximum is taken over all elements K in ∆ i . Then we can readily check that

N n=1 τ K∈T h |K||g(S n, ) -g(S n,r )| 2 ≤ C N n=1 τ M i=1 κ i j∈N (i) g(S n,j ) -g(S n,i ) 2 ,
where C is a bound for the maximum number of elements that share a common edge, bound independent of h and τ by virtue of the regularity of the mesh. Now, recall the classical formula in each d-simplex K,

K |∇ φ i • ∇ φ j | = 1 d 2 |F i ||F j | |K| |n i • n j |, (6.4)
where F i is the face opposite to the vertex a i and n i is the exterior (to K) unit normal to the face F i . The regularity of the mesh implies that there exists a constant c 0 , independent of h and τ , such that

|n i • n j | ≥ c 0 .
Hence, using again the regularity of the mesh, we obtain

K |∇ φ i • ∇ φ j | ≥ C h d-2 K ,
and denoting by ρ ij the minimum of h K for all K in ∆ i ∩ ∆ j , we deduce

c ij ≥ Cρ d-2 ij , (6.5)
with another constant C independent of h and τ . By collecting these results, we derive

I h (g(S h,τ )) -g(S h,τ ) 2 L 2 (Ω×]0,T [) ≤ C N n=1 τ M i=1 κ i j∈N (i) 1 ρ d-2 ij c ij g(S n,j ) -g(S n,i ) 2 . (6.6)
With another application of the regularity of the mesh, this becomes

I h (g(S h,τ )) -g(S h,τ ) 2 L 2 (Ω×]0,T [) ≤ Ch 2 ∇(I h (g(S h,τ ))) 2 L 2 (Ω×]0,T [) , (6.7)
(note that the power of h is independent of the dimension) and the limit (6.2) follows from Theorem 5.15.

6.1.2. Relation between g(S n,j ) -g(S n,i ) and S n,j -S n,i . Here, we derive an upper bound for S n,j -S n,i in terms of g(S n,j ) -g(S n,i ). Lemma 6.2. Under the assumptions (5.42)-(5.44) on the derivatives of the mobilities and capillary pressure, there exists a constant C, independent of h and τ , such that for all i, j, and n

|S n,j -S n,i | ≤ C |g(S n,j ) -g(S n,i )| 1 r , (6.8)
where r = max(θ o + β 4 , θ w + β 3 ) > 1.

Proof. To simplify, we set c = S n,i , d = S n,j and assume c < d. From (5.51), it follows that

g(d) -g(c) ≥ α 3 C max α w θ w α o θ o d c x θw+β 3 -1 (1 -x) θo+β 4 -1 . (6.9)
For the sake of brevity, we do not specify the constant factor in (6.9) and write

g(d) -g(c) ≥ C 1 d c x θw+β 3 -1 (1 -x) θo+β 4 -1 .
Now, we argue according to the positions of c and d. There are four cases.

1) If 1 8 ≤ c ≤ 7 8 , then (6.9) gives

g(d) -g(c) ≥ C 1 1 8 θw+β 3 -1 d c (1 -x) θo+β 4 -1 = C 1 1 8 θw+β 3 -1 1 θ o + β 4 (1 -c) θo+β 4 -(1 -d) θo+β 4 . But (1 -c) θo+β 4 -(1 -d) θo+β 4 = (1 -c) θo+β 4 -1 (d -c) + (1 -d) (1 -c) θo+β 4 -1 -(1 -d) θo+β 4 -1 ≥ (d -c) 1 8 θo+β 4 -1 . Hence g(d) -g(c) ≥ C 1 θ o + β 4 1 8 θw+β 3 +θo+β 4 -2 (d -c). (6.10)
2) If c > 7 8 , then d > 7 8 and (6.9) gives

g(d) -g(c) ≥ C 1 7 8 θw+β 3 -1 1 θ o + β 4 (1 -c) θo+β 4 -(1 -d) θo+β 4 . Let us set a = 1 -d, b = d -c γ = θ o + β 4 -1 > 0.
We can also write

(1 -c) θo+β 4 -(1 -d) θo+β 4 = a γ+1 1 + b a γ+1 -1 .
It is easy to check that the function

x → (1 + x) γ+1 -1 -x γ+1
vanishes at x = 0 and is strictly monotonic increasing, hence is strictly positive for x > 0. Hence

a γ+1 1 + b a γ+1 -1 > a γ+1 b a γ+1 = b γ+1 . Thus (1 -c) θo+β 4 -(1 -d) θo+β 4 ≥ (d -c) θo+β 4 ,
and g(d) -g(c) ≥ C 1 7 8

θw+β 3 -1 1 θ o + β 4 (d -c) θo+β 4 . (6.11)
3) If c < 1 8 and d < 7 8 , then the integrand 1 -x ≥ 1 -d > 1 8 and by the above argument, d -c)

g(d) -g(c) ≥ C 1 1 8 θo+β 4 -1 1 θ w + β 3 d θw+β 3 -c θw+β 3 ≥ C 1 1 8 θo+β 4 -1 1 θ w + β 3 (d -c) θw+β 3 . (6.12) 4) If c < 1 8 and d > 7 8 , then c < 1 6 (d -c) < 1 2 (d -c) < d. Therefore, we can write g(d) -g(c) ≥ C 1 1 2 (d-c) 1 6 (d-c) x θw+β 3 -1 (1 -x) θo+β 4 -1 ≥ C 1 1 2 θo+β 4 -1 1 θ w + β 3 1 2 (d -c) θw+β 3 - 1 6 (
θw+β 3 ≥ C 1 1 2 θo+θw+β 3 +β 4 -1 1 θ w + β 3 1 - 1 3 θw+β 3 (d -c) θw+β 3 . (6.13) Since d -c ≤ 1, θ o + β 4 > 1, and θ w + β 3 > 1, we have in all cases g(d) -g(c) ≥ C 2 (d -c) max(θo+β 4 ,θw+β 3 ) ,
where C 2 is the minimum of the constant factors in (6.10)-(6.13).

The convergence to zero of the differences I h (p αg (S h,τ )) -p αg (S h,τ ), α = w, o, follows from this lemma and Theorem 5.15.

Lemma 6.3. Under the assumptions (5.42)-(5.44) on the derivatives of the mobilities and capillary pressure, there exists a constant C, independent of h and τ , such that

I h (p αg (S h,τ )) -p αg (S h,τ ) L 2 (Ω×]0,T [) ≤ C h γα , α = w, o, (6.14)
where γ w = β 3 r , γ o = β 4 r and in both cases, r is the exponent of Lemma 6.2. Proof. Let us start with α = w. Arguing as in the proof of Lemma 6.1, with -p wg (monotonic increasing) instead of g, the analogue of (6.6) holds for -p wg (S h,τ ), with the same notation

I h (p wg (S h,τ )) -p wg (S h,τ ) 2 L 2 (Ω×]0,T [) ≤ C N n=1 τ M i=1 j∈N (i) κ i c ij c ij p wg (S n,j ) -p wg (S n,i ) 2 ,
(6.15) and the result will stem from an adequate upper bound for p wg (S n,j ) -p wg (S n,i ), for all neighbors j of i. To this end, we proceed as in Lemma 6.2. Let c = S n,i , d = S n,j and suppose again that c < d; then by (1.12), (5.44), (5.46), and (5.47),

|p wg (S n,j ) -p wg (S n,i )| ≤ 1 C min α 3 α o θ o d c x β 3 -1 (1 -x) θo+β 4 -1 , (6.16)
that we write as

|p wg (S n,j ) -p wg (S n,i )| ≤ C 1 d c x β 3 -1 (1 -x) θo+β 4 -1 .
Here, the discussion reduces to three cases.

1) If 1 8 ≤ c ≤ 7 8 , since θ o + β 4 -1 > 0, d c x β 3 -1 (1 -x) θo+β 4 -1 ≤ 8 1-β 3 d c (1 -x) θo+β 4 -1 ≤ 8 1-β 3 (d -c). (6.17) 2) Likewise, if c > 7 8 , d c x β 3 -1 (1 -x) θo+β 4 -1 ≤ 8 7 1-β 3 (d -c). (6.18) 3) If c < 1 8 , d c x β 3 -1 (1 -x) θo+β 4 -1 ≤ d c x β 3 -1 = 1 β 3 d β 3 -c β 3 ≤ 1 β 3 d -c β 3 . (6.19)
Indeed, by Jensen's inequality, valid for 0 < β 3 ≤ 1,

d = c + (d -c) ≤ c β 3 + (d -c) β 3 1 β 3 , i.e., d β 3 ≤ c β 3 + (d -c) β 3 .
Consequently, in all cases, |p wg (S n,j ) -p wg (S n,i )| ≤ C 2 S n,j -S n,i β 3 . (6.20)

Thus, by substituting into (6.15), applying Lemma 6.2, and setting γ w = β 3 r , we infer

I h (p wg (S h,τ )) -p wg (S h,τ ) 2 L 2 (Ω×]0,T [) ≤ C N n=1 M i=1 j∈N (i) κ i c ij τ c ij S n,j -S n,i | 2β 3 ≤ C N n=1 M i=1 j∈N (i) κ i c ij τ c ij A 2γw ij ,
where A ij = |g(S n,j ) -g(S n,i )|. Note that r > β 3 , hence γ w < 1. Then

I h (p wg (S h,τ ))-p wg (S h,τ ) 2 L 2 (Ω×]0,T [) ≤ C N n=1 M i=1 j∈N (i) κ i c ij (τ c ij ) 1-γw τ c ij γw A 2γw ij ≤ C N n=1 M i=1 j∈N (i) τ c ij A 2 ij γw N n=1 M i=1 j∈N (i) κ i c ij 1 1-γw τ c ij 1-γw . But N n=1 M i=1 j∈N (i) κ i c ij 1 1-γw τ c ij 1-γw ≤ C T |Ω| 1-γw sup i,j κ i c ij γw ,
and (6.14) with α = w follows from (6.5), the regularity of the mesh, and Theorem 5.15.

When α = o, the proof is based on fact that -p og is nonnegative, monotonically increasing, and satisfies the inequality

-p og (x) ≤ 1 C min 1 α 3 α w θ w x 0 x θw+β 3 -1 (1 -x) β 4 -1 .
By comparing with (6.16), we see that the above argument carries over to p og with β 3 replaced by β 4 .

Finally, with the notation of Lemma 6.3, the following bound regarding p c (S h,τ ) follows from (6.14) and (1.14), and the fact that p c (0) is a constant:

I h (p c (S h,τ )) -p c (S h,τ ) L 2 (Ω×]0,T [) ≤ C h γ , (6.21)
where γ = 1 r min(β 3 , β 4 ). Proof. The convergence (6.22) means that for all ψ ∈ L 1 (Ω×]0, T [),

Ω×]0,T [ S h,τ ψ → Ω×]0,T [ sψ and Ω×]0,T [ (1 -S h,τ )ψ → Ω×]0,T [ (1 -s)ψ.
We argue by contradiction. Suppose that s > 1 on a set of positive measure, say D, and take ψ = (s -1) + , the positive part of s -1. Then

0 ≤ Ω×]0,T [ 1 -S h,τ ψ → Ω×]0,T [ (1 -s)(s -1) + = D (1 -s)(s -1) + ,
thus contradicting the fact that (1 -s) < 0 on D. This proves that s ≤ 1. The proof that s ≥ 0 is similar.

Regarding the pressure, the bound (3.23) yields weak convergence, up to a subsequence, of the gradient of U α,h,τ . We can deduce weak convergence of the sequences themselves by applying (3.22). Indeed,

Ω U w,h,τ = U w,h,τ , 1 h = I h (p wg (S h,τ )), 1 h ,
owing to (2.45). Then the properties of p wg and the boundedness of S h,τ imply that

I h (p wg (S h,τ )), 1 h ≤ C.
Similarly,

Ω U o,h,τ = I h (p wg (S h,τ )) + p c (0), 1 h ,
a bounded quantity. Then we infer from (3.22) that

U α,h,τ L 2 (Ω×]0,T [) ≤ C, α = w, o. (6.24)
With this, (3.23) implies that there exist functions Wα ∈ L 2 (0, T ; H 1 (Ω)), α = w, o, and a subsequence of h and τ (not indicated) such that, lim

(h,τ )→(0,0) U α,h,τ = Wα , weakly in L 2 (0, T ; H 1 (Ω)). (6.25)
Likewise, the function I h (g(S h,τ )) is bounded in L 2 (Ω×]0, T [) and it follows from this and (5.58) that there exists a function K ∈ L 2 (0, T, H 1 (Ω)) such that, up to a subsequence, lim (h,τ )→(0,0)

I h (g(S h,τ )) = K weakly in L 2 (0, T, H 1 (Ω)). (6.26)
This implies in particular that for almost every time t, I h (g(S h,τ )) converges strongly in L 2 (Ω). But as is well-known, these convergences are not sufficient to pass to the limit in the nonlinear terms: they must be supplemented by a bound for a fractional derivative in time of S h,τ that yields compactness in time. This will stem via a bound for a fractional derivative in time of g(S h,τ ).

Compactness in time.

Following the argument introduced by Kazhikhov, see [START_REF] Lions | On some questions in boundary value problems of mathematical physics[END_REF], and recalling that • ϕ h is equivalent to the L 2 norm in finite dimension, we want to derive first a fractional estimate in time for I h (g(S h,τ )) and next for g(S h,τ ). The following lemma is a preliminary bound written in terms of sums of the pointwise values in time. Lemma 6.5. Under the assumptions of Theorem 5.15, there exist constants C, independent of h and τ , such that for all integers 1 ≤ ≤ N -1,

N - m=1 τ g(S m+ h ) -g(S m h ) ϕ h 2 ≤ C( τ ), N - m=1 τ g(S m+ h ) -g(S m h ) 2 L 2 (Ω) ≤ C( τ ). (6.27)
Proof. The starting point is the inequality

N - m=1 τ g(S m+ h ) -g(S m h ) ϕ h 2 ≤ L N - m=1 τ g(S m+ h ) -g(S m h ), S m+ h -S m h ϕ h , (6.28)
owing that g is Lipschitz continuous and increasing. Thus, by writing

S m+ h -S m h = k=1 S m+k h -S m+k-1 h ,
testing each line of (2.42) taken at level m + k with I h g(S m+ h ) -g(S m h ) , and applying (6.28), we obtain

N - m=1 τ g(S m+ h ) -g(S m h ) ϕ h 2 ≤ L N - m=1 τ k=1 τ f w (s m+k in,h )q m+k h -f w (S m+k h )q m+k h , g(S m+ h ) -g(S m h ) h + P m+k w,h , I h (η w (S m+k h )); P m+k w,h , I h g(S m+ h ) -g(S m h ) h . (6.29)
It is easy to check that, on one hand, with r = or r = 0,

P m+k w,h , I h (η w (S m+k h )); P m+k w,h , I h g(S m+r h ) h = 1 2 M i,j=1
g(S m+r,j ) -g(S m+r,i ) c ij η w (S m+k,ij w

) P m+k,j w -P m+k,i w ≤ 1 4 M i,j=1 c ij η w (S m+k,ij w ) g(S m+r,j ) -g(S m+r,i ) 2 + P m+k,j w -P m+k,i w 2 ≤ 1 4 M i,j=1 c ij η w (1) g(S m+r,j ) -g(S m+r,i ) 2 + η w (S m+k,ij w ) P m+k,j w -P m+k,i w 2 ,
since η w is increasing and S h,τ is bounded by one. On the other hand,

f w (s m+k in,h )q m+k h -f w (S m+k h )q m+k h , g(S m+ h ) -g(S m h ) h ≤ C qm+k L 1 (Ω) + q m+k L 1 (Ω) ,
where here and below, C denotes constants that are independent of , h, and τ . Therefore, in view of (2.10)

N - m=1 τ g(S m+ h ) -g(S m h ) ϕ h 2 ≤ L N - m=1 τ 1 8 η w (1)( τ ) r= ,0 ∇ I h (g(S m+r h )) 2 L 2 (Ω) + 1 2 k=1 τ M i,j=1 c ij η w (S m+k,ij w ) P m+k,j w -P m+k,i w 2 + C k=1 τ qm+k L 1 (Ω) + q m+k L 1 (Ω) ≤ 1 8 η w (1)L( τ ) N m=1+ τ ∇ I h (g(S m h )) 2 L 2 (Ω) + N - m=1 τ ∇ I h (g(S m h )) 2 L 2 (Ω) + 1 2 L N - m=1 τ k=1 τ M i,j=1 c ij η w (S m+k,ij w ) P m+k,j w -P m+k,i w 2 + C qm+k L 1 (Ω) + q m+k L 1 (Ω)
.

(6.30) By (5.58), il suffices to bound the terms in the last line above. This is achieved by interchanging the sums over m and k. Let n = m + k; n runs from 2 to N and m runs from max(1, n -) to min(n -1, N -). Thus

N - m=1 τ k=1 τ M i,j=1 c ij η w (S m+k,ij w ) P m+k,j w -P m+k,i w 2 = N n=2 τ min(n-1,N -) m=max(1,n-) τ M i,j=1 c ij η w (S n,ij w ) P n,j w -P n,i w 2 . But min(n -1, N -) -max(1, n -) ≤ -1. Hence N - m=1 τ k=1 τ M i,j=1 c ij η w (S m+k,ij w ) P m+k,j w -P m+k,i w 2 ≤ ( τ ) N n=2 τ M i,j=1
c ij η w (S n,ij w ) P n,j w -P n,i w 2 , (6.31) and we know from (3.18) that this last sum over n is bounded. In the same fashion,

N - m=1 τ k=1 τ qm+k L 1 (Ω) + q m+k L 1 (Ω) ≤ ( τ ) q L 1 (Ω×]0,T [) + q L 1 (Ω×]0,T [) . (6.32)
Then, under the assumptions of Theorem 5.15, (6.27) follows by substituting (5.58), (6.31), and (6.32) into (6.30). The second inequality stems from the first and (6.1).

The next theorem transforms (6.27) into integrals.

Theorem 6.6. Under the assumptions of Theorem 5.15, there exists a constant C, independent of h, and τ , such that for all real numbers δ, 0 < δ < T ,

T -δ 0 g(S h,τ (t + δ)) -g(S h,τ (t) 2 L 2 (Ω) dt ≤ Cδ. (6.33) Similarly, T -δ 0 I h g(S h,τ (t + δ)) -g(S h,τ (t) 2 L 2 (Ω) dt ≤ Cδ, (6.34)
with another constant C, independent of h, and τ .

Proof. The argument is not new, see for instance [START_REF] Walkington | Compactness properties of the DG and CG time stepping schemes for parabolic equations[END_REF], but we recall it for the reader's convenience. The discussion depends on the value of δ; there are three cases:

(i) 0 < δ < τ , (ii) δ = τ , 1 ≤ ≤ N -1, (iii) δ = τ ( + η), 1 ≤ ≤ N -1, 0 < η < 1. (i) If 0 < δ < τ , we have for all f in L 1 (0, T ) T -δ 0 f (t) dt = N -2 m=0 t m+1 -δ tm f (t) dt + t m+1 t m+1 -δ f (t) dt + t N -δ t N -1 f (t) dt.
On the interval (t m , t m+1 -δ), by convention, see (2.38)

applied to S h,τ , S h,τ (t + δ) = S m+1 h = S h,τ (t) and on the interval (t m+1 -δ, t m+1 ), S h,τ (t + δ) = S m+2 h . Therefore T -δ 0 g(S h,τ (t + δ)) -g(S h,τ (t)) 2 L 2 (Ω) dt = N -2 m=0 δ g(S m+2 h ) -g(S m+1 h ) 2 L 2 (Ω) = N -1 m=1 δ g(S m+1 h ) -g(S m h ) 2 L 2 (Ω) ,
and (6.27) with = 1 yields (6.33) for this value of δ.

(ii) Let δ = τ , for instance consider = 2. We have

t m+1 tm g(S h,τ (t + 2τ )) -g(S h,τ (t)) 2 L 2 (Ω) dt = τ g(S m+3 h ) -g(S m+1 h ) 2 L 2 (Ω) . Thus T -2τ 0 g(S h,τ (t + δ)) -g(S h,τ (t)) 2 L 2 (Ω) dt = τ N -2 m=1 g(S m+2 h ) -g(S m h ) 2 L 2 (Ω) ,
and (6.27) with = 2 yields (6.33). This argument applies to any with 1 ≤ ≤ N -1.

(iii) Let δ = τ ( + η) for some 0 < η < 1, and for instance consider again = 2. Then for any integer m, 0 ≤ m ≤ N -3,

t m+1 -τ η tm g(S h,τ (t + τ (2 + η)) -g(S h,τ (t)) 2 L 2 (Ω) = τ (1 -η) g(S m+3 h ) -g(S m+1 h ) 2 L 2 (Ω) ,
and

t m+1 t m+1 -τ η g(S h,τ (t + τ (2 + η)) -g(S h,τ (t)) 2 L 2 (Ω) = τ η g(S m+4 h ) -g(S m+1 h ) 2 L 2 (Ω) . Hence T -τ (2+η) 0 g(S h,τ (t + τ (2 + η))) -g(S h,τ (t)) 2 L 2 (Ω) dt = N -2 m=1 τ (1 -η) g(S m+2 h ) -g(S m h ) 2 L 2 (Ω) + N -3 m=1 τ η g(S m+3 h ) -g(S m h ) 2 L 2 (Ω) .
Then (6.27) with = 2 and = 3 implies that

T -τ (2+η) 0 g(S h,τ (t + τ (2 + η)) -g(S h,τ (t)) 2 L 2 (Ω) dt ≤ C 2τ (1 -η) + 3τ η = Cτ (2 + η).
More generally, we have

T -τ ( +η) 0 g(S h,τ (t + τ ( + η)) -g(S h,τ (t)) 2 L 2 (Ω) dt = τ (1 -η) N - m=1 g(S m+ h ) -g(S m h ) 2 L 2 (Ω) + τ η N -( +1) m=1 g(S m+ +1 h ) -g(S m h ) 2 L 2 (Ω) ,
and we apply (6.27) with to the first sum and + 1 to the second sum. This proves the first part of (6.33). The proof of (6.34) follows likewise from (6.27) with ϕ = 1.

6.4. Strong convergence. With Theorem 6.6, it follows from Kolmogorov's theorem that the sequence I h (g(S h,τ )) is compact in L 2 (Ω×]0, T [), see [START_REF] Lions | On some questions in boundary value problems of mathematical physics[END_REF]. Thus, again up to a subsequence, I h (g(S h,τ )) converges strongly in L 2 (Ω×]0, T [). Since it converges weakly to K in L 2 (0, T ; H 1 (Ω)) ( K belongs also to L ∞ (Ω×]0, T [)), uniqueness of the limit implies lim (h,τ )→(0,0)

I h (g(S h,τ )) = K strongly in L 2 (Ω×]0, T [). (6.35)
By Lemma 6.1, this also implies lim

(h,τ )→(0,0) g(S h,τ ) = K strongly in L 2 (Ω×]0, T [). (6.36)
From here, let us prove the strong convergence of S h,τ . Recall that g is invertible with range ]0, β[ and inverse g

-1 ∈ W 1,∞ (]0, β[). Let F h,τ = g(S h,τ ); then S h,τ = g -1 (F h,τ ).
The strong convergence of F h,τ and the continuity of g -1 imply the strong convergence of S h,τ to g -1 ( K) in L 2 (Ω×]0, T [), and since S h,τ converges weakly to s, uniqueness of the limit implies that s = g -1 ( K), i.e., lim (h,τ )→(0,0) S h,τ = s = g -1 ( K) strongly in L 2 (Ω×]0, T [). (6.37) This strong convergence and the continuity of g, p αg , α = w, o, and p c , also imply that lim

(h,τ )→(0,0) g(S h,τ ) = g(s), lim (h,τ )→(0,0) p αg (S h,τ ) = p αg (s), α = w, o, lim (h,τ )→(0,0) p c (S h,τ ) = p c (s), (6.38
) all strongly in L 2 (Ω×]0, T [). Furthermore Lemma 6.3 and (6.21) yield lim (h,τ )→(0,0)

I h (p αg (S h,τ )) = p αg (s), lim (h,τ )→(0,0) I h (p c (S h,τ )) = p c (s), strongly in L 2 (Ω×]0, T [).
(6.39) In view of (6.25), this convergence implies that P α,h,τ converges weakly in L 2 (Ω×]0, T [) to some function pα ∈ L 2 (Ω×]0, T [), α = w, o. Furthermore, uniqueness of the limit implies that Wα , the limit function of U α,h,τ has the form Ww = pw + p wg (s), Wo = po -p og (s).

(6.40) 7. Identification of the limit. Let us pass to the limit in the equations of the scheme. This is done in several steps because we do not have convergence of the pressure gradient.

7.1. The upwind diffusions. Since the discrete auxiliary pressures U α,h,τ converge weakly to Wα in L 2 (0, T ; H 1 (Ω)), instead of treating directly the upwind diffusion terms P α,h,τ , I h (η α (S h,τ )); P α,h,τ , θ h h , we begin with P α,h,τ , I h (η α (S h,τ )); U α,h,τ , θ h h . 7.1.1. Discrete auxiliary pressure. Let us start with the wetting phase, the treatment of the non-wetting phase being much the same.

Let v be a smooth function, say v ∈ C 1 ( Ω × [0, T ]) and let V h,τ = ρ τ (I h (v)). Assume for the moment that s, the limit of S h,τ , is sufficiently smooth, say s ∈ W 1,∞ (Ω×]0, T [) and let sτ = s(t n ) in ]t n-1 , t n ]. Then assumption (5.42) implies

1 τ tn t n-1 η w (s) dt -η w (s n τ ) L ∞ (Ω) ≤ Cτ η w L ∞ (0,1) ∂ t s L ∞ (Ω×]0,T [) . (7.1)
Proceeding as in Section 2.2, we treat the upwinding in several steps and consider first

T 0 Ω η w (s)∇ U w,h,τ • ∇ V h,τ = T 0 Ω ∇ U w,h,τ • ∇ V h,τ ρ τ (η w (s)) -η w (s τ ) + η w (s τ ) . (7.2)
But in view of (7.1),

T 0 Ω ∇ U w,h,τ • ∇ V h,τ ρ τ (η w (s)) -η w (s τ )) ≤ Cτ η w L ∞ (0,1) ∂ t s L ∞ (Ω×]0,T [) × U w,h,τ L 2 (0,T ;H 1 (Ω) V h,τ L 2 (0,T ;H 1 (Ω)) ,
and the boundedness of all factors of τ , owing to (3.23) and the regularity of v, implies lim (h,τ )→(0,0)

T 0 Ω ∇ U w,h,τ • ∇ V h,τ ρ τ (η w (s)) -η w (s τ )) = 0. (7.3)
Next the weak convergence of U w,h,τ to Ww in L 2 (0, T ; H 1 (Ω)), the strong convergence of V h,τ to v in L ∞ (0, T ; W 1,∞ (Ω)), the continuity of η w , the regularity of s, and (7.3) imply lim (h,τ )→(0,0)

T 0 Ω η w (s τ )∇ U w,h,τ • ∇ V h,τ = T 0 Ω η w (s)∇ Ww • ∇ v.
Let us expand the expression in the above left-hand side. With the notation (2.11), in view of Proposition 2.2 we have

T 0 Ω η w (s τ )∇ U w,h,τ • ∇ V h,τ = - N n=1 τ M i,j=1 U n,i w K⊂∆ i ∩∆ j c ij,K (η w (s n τ )) K V n,j -V n,i .
By symmetry, this becomes

T 0 Ω η w (s n τ )∇ U w,h,τ •∇ V h,τ = 1 2 N n=1 τ M i,j=1 K⊂∆ i ∩∆ j c ij,K (η w (s n τ )) K U n,j w -U n,i w V n,j -V n,i . (7.4) Hence lim (h,τ )→(0,0) 1 2 N n=1 τ M i,j=1 K⊂∆ i ∩∆ j c ij,K (η w (s n τ )) K U n,j w -U n,i w V n,j -V n,i = T 0 Ω η w (s)∇ Ww •∇ v.
(7.5) According to (5.42) and the regularity of s, η w (s) belongs to L ∞ (0, T ; W

R ij = K⊂∆ i ∩∆ j c ij,K (η w (s n τ )) K -η w (s n,ij w,τ ) .
By proceeding as in Theorem 2.4 and applying (3.23), the regularity of v, and the approximation properties of I h , we obtain

1 2 N n=1 τ M i,j=1 R ij U n,j w -U n,i w V n,j -V n,i ≤ 1 2 N n=1 τ M i,j=1 |R ij | U n,j w -U n,i w 2 1 2 M i,j=1 |R ij | V n,j -V n,i 2 1 2 ≤ C h η w L ∞ (0,1) ∇ s L ∞ (Ω×]0,T [) ∇ U w,h,τ L 2 (Ω×]0,T [) ∇ V h,τ L 2 (Ω×]0,T [) ≤ C h η w L ∞ (0,1) ∇ s L ∞ (Ω×]0,T [) |v| H 1 (0,T ;H 2 (Ω)) .
With (7.5), this implies lim (h,τ )→(0,0)

1 2 N n=1 τ M i,j=1 c ij η w (s n,ij w,τ ) U n,j w -U n,i w V n,j -V n,i = T 0 Ω η w (s)∇ Ww •∇ v. (7.7)
To recover T 0 P w,h,τ , I h (η w (S h,τ )); U w,h,τ , V h,τ h , we write η w (s n,ij w,τ ) = η w (s n,ij w,τ ) -η w (S n,ij w ) + η w (S n,ij w ), and we must examine the convergence of

X := 1 2 N n=1 τ M i,j=1 c ij η w (s n,ij w,τ ) -η w (S n,ij w ) U n,j w -U n,i w V n,j -V n,i .
On the one hand, owing to the smoothness of v, we have

|V n,j -V n,i | ≤ Ch i ∇ v L ∞ (Ω×]0,T [) , (7.8)
where h i is the length of the edge whose endpoints are the vertices i and j. On the other hand,

|η w (s n,ij w,τ ) -η w (S n,ij w )| ≤ C η w L ∞ (0,1) |s n,ij w,τ -S n,ij w |. Hence |X| ≤ C ∇ v L ∞ (Ω×]0,T [) ∇ U w,h,τ L 2 (Ω×]0,T [) N n=1 τ M i,j=1 c ij h 2 i |s n,ij w,τ -S n,ij w | 2 1 2 .
It is easy to check that M i,j=1

c ij h 2 i |s n,ij w,τ -S n,ij w | 2 ≤ C M i=1 m i |s n,i τ -S n,i | 2 . Therefore |X| ≤ C ∇ v L ∞ (Ω×]0,T [) ∇ U w,h,τ L 2 (Ω×]0,T [) N n=1 τ I h (s n τ ) -S n h,τ 2 L 2 (Ω) 1 2 = C ∇ v L ∞ (Ω×]0,T [) ∇ U w,h,τ L 2 (Ω×]0,T [) I h (s τ ) -S h,τ L 2 (Ω×]0,T [) ,
where we have used the equivalence (2.5). Then, we write

I h (s τ ) -S h,τ L 2 (Ω×]0,T [) ≤ I h (s τ ) -sτ L 2 (Ω×]0,T [) + sτ -s L 2 (Ω×]0,T [) + s -S h,τ L 2 (Ω×]0,T [) ,
and the approximation properties of I h , the strong convergence of sτ to s and of S h,τ to s, all in L 2 (Ω×]0, T [) imply that lim (h,τ )→(0,0)

N n=1 τ M i,j=1 c ij η w (s n,ij w,τ ) -η w (S n,ij w ) U n,j w -U n,i w V n,j -V n,i = 0. (7.9)
A combination of (7.9) and (7.7) yields the intermediate convergence result when the limit function s is smooth, lim (h,τ )→(0,0)

- N n=1 τ P w,h,τ , I h (η w (S h,τ )); U w,h,τ , V h,τ h = T 0 Ω η w (s)∇ Ww • ∇ v. (7.10)
It remains to lift the regularity restriction on s. Let (S m ) m≥1 be a sequence of smooth functions that tend to s in L 2 (Ω×]0, T [). Then for each ε > 0, there exists an integer M 0 such that

S M 0 -s L 2 (Ω×]0,T [) ≤ ε. (7.11)
From (7.11), the projection properties, and the fact that M 0 is fixed, we infer

ρ τ (η w (s)) -η w (s) L 2 (Ω×]0,T [) ≤ ρ τ (η w (s) -η w (S M 0 )) L 2 (Ω×]0,T [) + ρ τ (η w (S M 0 )) -η w (S M 0 ) L 2 (Ω×]0,T [) + η w (S M 0 ) -η w (s) L 2 (Ω×]0,T [) ≤ (2 ε + C τ ) η w L ∞ (0,1) . (7.12)
Now, we replace (7.2) by

T 0 Ω η w (s)∇ U w,h,τ • ∇ V h,τ = T 0 Ω ρ τ η w (s) -η w (S M 0 ) ∇ U w,h,τ • ∇ V h,τ + T 0 Ω ρ τ η w (S M 0 ) ∇ U w,h,τ • ∇ V h,τ = T 0 Ω ρ τ η w (s) -η w (S M 0 ) ∇ U w,h,τ • ∇ V h,τ + T 0 Ω η w (S M 0 )∇ U w,h,τ • ∇ V h,τ . (7.13)
For the first term, owing to (7.11), the projection properties, and (5.42), we have

T 0 Ω ρ τ η w (s)-η w (S M 0 ) ∇ U w,h,τ • ∇ V h,τ ≤ ∇ V h,τ L ∞ (Ω×]0,T [) × ∇ U w,h,τ L 2 (Ω×]0,T [) η w (s) -η w (S M 0 ) L 2 (Ω×]0,T [) ≤ ε η w L ∞ (0,1) ∇ V h,τ L ∞ (Ω×]0,T [) ∇ U w,h,τ L 2 (Ω×]0,T [) .
Then the uniform boundedness of U w,h,τ and V h,τ yield

T 0 Ω ρ τ η w (S) -η w (S M 0 ) ∇ U w,h,τ • ∇ V h,τ ≤ C ε, (7.14)
with a constant C independent of h and τ . Thus, we must examine the limit of the second term. Since M 0 is fixed and S M 0 is smooth, by reproducing the previous steps, we derive the analogue of (7.7) for the function S M 0 , lim (h,τ )→(0,0)

1 2 N n=1 τ M i,j=1 c ij η w ((S M 0 ) n,ij w,τ ) U n,j w -U n,i w V n,j -V n,i = T 0 Ω η w (S M 0 )∇ Ww • ∇ v = T 0 Ω η w (s)∇ Ww • ∇ v + R, (7.15)
where

|R| = T 0 Ω η w (S M 0 ) -η w (s) ∇ Ww • ∇ v ≤ η w L ∞ (0,1) S M 0 -s L 2 (Ω×]0,T [) ∇ Ww L 2 (Ω×]0,T [) ∇ v L ∞ (Ω×]0,T [) ≤ C ε. (7.16)
To relate the left-hand side of (7.15) to P w,h,τ , I h (η w (S h,τ )); U w,h,τ , V h h , we split

η w ((S M 0 ) n,ij w,τ ) = η w (S n,ij w ) + η w ((S M 0 ) n,ij w,τ ) -η w (S n,ij w ),
and examine the convergence of

Y := 1 2 N n=1 τ M i,j=1 c ij η w ((S M 0 ) n,ij w,τ ) -η w (S n,ij w ) U n,j w -U n,i w V n,j -V n,i .
By arguing as above and using the interpolant I h , we derive

|Y | ≤ C η w L ∞ (0,1) ∇ v L ∞ (Ω×]0,T [) ∇ U w,h,τ L 2 (Ω×]0,T [) I h ((S M 0 ) τ ) -S h,τ L 2 (Ω×]0,T [) .
Finally, we write

I h ((S M 0 ) τ ) -S h,τ L 2 (Ω×]0,T [) ≤ I h ((S M 0 ) τ ) -(S M 0 ) τ L 2 (Ω×]0,T [) + (S M 0 ) τ -S M 0 L 2 (Ω×]0,T [) + S M 0 -s L 2 (Ω×]0,T [) + s -S h,τ L 2 (Ω×]0,T [) ≤ C h S M 0 L ∞ (0,T ;H 2 (Ω)) + Cτ S M 0 H 1 (0,T ;L 2 (Ω)) + ε + s -S h,τ L 2 (Ω×]0,T [) , so that |Y | ≤ C(h + τ + ε) + s -S h,τ L 2 (Ω×]0,T [) .
(7.17)

In the next theorem, the limit (7.10) when s is only in L 2 (Ω×]0, T [) follows by combining (7.13)-(7.17). The same argument holds when w is replaced by o.

Theorem 7.1. Let v ∈ C 1 ( Ω × [0, T ]) be a smooth function and let V h,τ = I h (v)(t n ) in ]t n-1 , t n ].
Under the assumptions and notation on the mobility (5.42)-(5.46), lim (h,τ )→(0,0)

- T 0 P α,h,τ , I h (η α (S h,τ )); U α,h,τ , V h,τ h = T 0 Ω η α (s)∇ Wα • ∇ v (7.18)
where s is the strong limit of S h,τ and Wα the weak limit of U α,h,τ , α = w, o.

7.1.2. The additional term . This paragraph is dedicated to the limit of

T 0 P α,h,τ , I h (η α (S h,τ )); I h (p αg (S h,τ )), V h,τ h , α = w, o.
It shall be split below, as suggested by the following observation, derived from (1.12) and (1.13):

η w (S ij w )p wg (S j ) + g(S j ) = S j 0 f o (x) η w (S ij w ) -η w (x) p c (x) dx, η o (S ij o )p og (S j ) + g(S j ) = S j 0 f w (x) η o (S ij o ) -η o (x) p c (x) dx.
Thus, we add and subtract g and write by applying (2.9),

T 0 P α,h,τ , I h (η α (S h,τ )); I h (p αg (S h,τ )), V h,τ h = N n=1 τ M i,j=1
V n,i c ij η α (S n,ij α ) p αg (S n,j ) -p αg (S n,i ) + g(S n,j ) -g(S n,i )

+ T 0 Ω ∇ g(S h,τ ) • ∇ V h,τ = T 1 + T 2 .
Since lim (h,τ )→(0,0)

T 2 = T 0 Ω ∇ g(s) • ∇ v, (7.19)
we must prove that the first term tends to zero. When α = w, it has the form

T 1 = - 1 2 N n=1 τ M i,j=1 c ij S n,j S n,i f o (x) η w (S n,ij w ) -η w (x) p c (x) dx V n,j -V n,i , (7.20)
with an analogous expression in the non-wetting phase. Then (7.8) yields,

|T 1 | ≤ C 2 ∇ v L ∞ (Ω×]0,T [) N n=1 τ M i,j=1 h i c ij S n,j S n,i f o (x) η w (S n,ij w ) -η w (x) p c (x) dx . (7.21)
Showing that T 1 is small requires a technical argument that we split into several steps.

Proposition 7.2. For the wetting phase, we have

| S j S i f o (x) η w (S ij w ) -η w (x) p c (x) dx ≤ -η w (S j ) -η w (S i ) p wg (S j ) -p wg (S i ) . (7.22)
For the non-wetting phase, the corresponding expression is bounded by

| S j S i f w (x) η o (S ij o ) -η o (x) p c (x) dx ≤ η o (S j ) -η o (S i ) p og (S j ) -p og (S i ) . (7.23)
Proof. Let us prove (7.22), the proof of (7.23) being similar. The discussion depends on the respective values of S j and S i . There are two cases: S i < S j or S i > S j . Of course S i = S j brings nothing.

1) If S i < S j and S ij w = S i , then η w (S ij w )-η w (x) = η w (S i )-η w (x), and, as p wg is decreasing,

0 ≤ S j S i f o (x)(-p c (x)) η w (x) -η w (S i ) dx ≤ -η w (S j ) -η w (S i ) p wg (S j ) -p wg (S i ) . If S ij w = S j , then η w (S ij w ) -η w (x) = η w (S j ) -η w (x), and 0 ≤ S j S i f o (x)(-p c (x)) η w (S j ) -η w (x) dx ≤ -η w (S j ) -η w (S i ) p wg (S j ) -p wg (S i ) . 2) If S i > S j and S ij w = S i , then 0 ≤ S i S j f o (x)(-p c (x)) η w (S i ) -η w (x) dx ≤ -η w (S i ) -η w (S j ) p wg (S i ) -p wg (S j ) .
Finally, suppose that S ij w = S j . Then 0

≤ S i S j f o (x) p c (x) η w (S j ) -η w (x) dx ≤ -η w (S i ) -η w (S j ) p wg (S i ) -p wg (S j ) .
This proves (7.22).

By substituting (7.22) into (7.21), we arrive at

|T 1 | ≤ C 2 ∇ v L ∞ (Ω×]0,T [) N n=1 τ M i,j=1
h i c ij -η w (S n,j ) -η w (S n,i ) p wg (S n,j ) -p wg (S n,i ) , (7.24) with an analogous bound in the non-wetting phase. Up to the factor h i , they behave like

T 0 Ω ∇(I h (η α (S h,τ ))) • ∇(I h (p αg (S h,τ ))), α = w, o.
Thus T 1 tends to zero if this quantity is bounded or is of the order of a negative power of h that is larger than -1. We have no direct bound for it, but as we do have a bound for (5.40), we can gain some insight by relating the two integrands. Again, we examine the wetting phase, the treatment of the non-wetting phase being the same. The proposition below will be applied to x 1 = S n,i and x 2 = S n,j . The condition x 1 < x 2 is not a restriction because if it does not hold, the indices i and j can be interchanged without changing the value of the two integrands.

T 0 Ω ∇(I h (f α (S h,τ ))) • ∇(I h (g(S h,τ ))), see
Proposition 7.3. Under the assumptions and notation on the mobility (5.42)-(5.46), we have for all pairs x 1 , x 2 with 0 ≤

x 1 < x 2 ≤ 3 4 , η w (x 2 ) -η w (x 1 ) p wg (x 1 ) -p wg (x 2 ) ≤ C(x θw 2 -x θw 1 )(x β 3 2 -x β 3 1 ), (7.25) f w (x 2 ) -f w (x 1 ) g(x 2 ) -g(x 1 ) ≥ C(x θw 2 -x θw 1 )(x θw+β 3 2 -x θw+β 3 1
).

(7.26)

Similarly, we have for all pairs x 1 , x 2 with 1 4 ≤ x 1 < x 2 ≤ 1,

η w (x 2 ) -η w (x 1 ) p wg (x 1 ) -p wg (x 2 ) ≤ C(x 2 -x 1 ) (1 -x 1 ) θo+β 4 -(1 -x 2 ) θo+β 4 , (7.27) f w (x 2 ) -f w (x 1 ) g(x 2 ) -g(x 1 ) ≥ C (1 -x 1 ) θo -(1 -x 2 ) θo (1 -x 1 ) θo+β 4 -(1 -x 2 ) θo+β 4 .
(7.28) Finally, we have for all pairs x 1 , x 2 with 0 ≤ x 1 ≤ 1 4 and 3 4 ≤ x 2 ≤ 1,

η w (x 2 ) -η w (x 1 ) p wg (x 1 ) -p wg (x 2 ) ≤ C f w (x 2 ) -f w (x 1 ) g(x 2 ) -g(x 1
) . (7.29)

All constants C above are independent of x 1 and x 2 .

Proof. According to (5.42),

η w (x 2 ) -η w (x 1 ) ≤ 1 α w θ w (x θw 2 -x θw 1 ).
Next, recalling that p wg (x) = f o (x)p c (x), we have, owing to (5.44), (5.46), and (1.8),

p wg (x 1 ) -p wg (x 2 ) = x 2 x 1 f o (x)(-p c (x)) dx ≤ 1 η * 1 α 3 1 α o θ o x 2 x 1 x β 3 -1 (1 -x) θo+β 4 -1 dx ≤ 1 η * 1 α 3 1 α o θ o x 2 x 1 x β 3 -1 dx ≤ 1 η * 1 α 3 β 3 1 α o θ o (x β 3 2 -x β 3 1
), (7.30) and (7.25), valid on [0, 1], follows from these two inequalities.

For (7.26), we use (5.53) that gives

f w (x 2 ) -f w (x 1 ) ≥ α o α w C 2 max 1 θ o θ w 1 4 θo (x θw 2 -x θw 1 ), (7.31)
and we use (5.51) that gives

g(x 2 ) -g(x 1 ) ≥ α 3 C max α w θ w α o θ o 1 θ w + β 3 1 4 θo+β 4 -1 (x θw+β 3 2 -x θw+β 3 1
).

The product of the two leads to (7.26).

Regarding (7.27), (7.25), albeit valid for all x ∈ [0, 1], is not adequate for the comparison we have in mind, and instead we use that

η w (x) ≤ 1 α w , which implies that η w (x 2 ) -η w (x 1 ) ≤ 1 α w (x 2 -x 1 ).
Similarly, we use

-p wg (x) ≤ 1 η * 1 α 3 1 α o θ o 4 1-β 3 (1 -x) θo+β 4 -1 , so that p wg (x 1 ) -p wg (x 2 ) ≤ 1 η * 1 α 3 1 α o θ o 1 θ o + β 4 4 1-β 3 (1 -x 1 ) θo+β 4 -(1 -x 2 ) θo+β 4 ),
thus proving (7.27). Next, by applying (5.54), we have

f w (x 2 ) -f w (x 1 ) ≥ 1 C 2 max α o α w θ w θ o 1 4 θw (1 -x 1 ) θo -(1 -x 2 ) θo .
Likewise, by applying (5.51), we obtain

g(x 2 ) -g(x 1 ) ≥ α 3 C max α w θ w α o θ o 1 θ o + β 4 1 4 θw-1+β 3 (1 -x 1 ) θo+β 4 -(1 -x 2 ) θo+β 4 .
The product of the two yields (7.28).

Finally, when 0 ≤ x 1 ≤ 1 4 and 3 4 ≤ x 2 ≤ 1, since both η w and -p wg are both increasing, they satisfy

η w (x 2 ) -η w (x 1 ) p wg (x 1 ) -p wg (x 2 ) ≤ η w (1) -p wg (1) > 0.
Likewise, as both f w and g are increasing, they satisfy

f w (x 2 ) -f w (x 1 ) g(x 2 ) -g(x 1 ) ≥ f w ( 3 4 ) -f w ( 1 4 ) g( 3 4
) -g(

) =: D > 0.

Hence

η w (x 2 ) -η w (x 1 ) p wg (x 1 ) -p wg (x 2 ) ≤ - 1 D η w p wg (1) f w (x 2 ) -f w (x 1 ) g(x 2 ) -g(x 1
) , whence (7.29). Clearly all constants involved are independent of x 1 and x 2 .

It stems from (7.25) and (7.26), that the two left-hand sides cannot be compared when x 1 and x 2 are too small. The same observation applies to (7.27) and (7.28) when 1 -x 1 and 1 -x 2 are too small. But in this case, there is no need for comparison because the expression we want to bound is sufficiently small, as is shown in the next proposition where again, x 1 = S n,i and x 2 = S n,j .

Proposition 7.4. Suppose that x 1 < x 2 ≤ h γ 1 i for some exponent γ 1 > 0 such that

γ 1 > 1 θ w + β 3 . (7.32) Then h i η w (x 2 ) -η w (x 1 ) p wg (x 1 ) -p wg (x 2 ) ≤ C h 2 i h δ 1 i , (7.33)
with the constant C of (7.26), whereas

η w (x 2 ) -η w (x 1 ) p wg (x 1 ) -p wg (x 2 ) ≤ Cx θw+β 3 2 ,
with the constant C of (7.25). Hence

h i η w (x 2 ) -η w (x 1 ) p wg (x 1 ) -p wg (x 2 ) ≤ C h i x θw 2 f w (x 2 ) -f w (x 1 ) g(x 2 ) -g(x 1 ) ,
with another constant C independent of x 1 , x 2 , and h i . Now, we use the assumption that x 2 > h γ 1 i . Then, owing to (7.40),

h i x θw 2 ≤ h 1-γ 1 θw i ≤ h δ 1 i ,
and we recover (7.39).

2) When x 1 > 1 2 x 2 , we infer from the next to last inequality in (7.30) that

p wg (x 1 ) -p wg (x 2 ) ≤ 1 η * 1 α 3 1 α o θ o x 2 -x 1 x β 3 -1 1 ≤ 1 η * 1 α 3 1 α o θ o 2 1-β 3 1 x 1-β 3 2 x 2 -x 1 .
Thus, on the one hand,

η w (x 2 ) -η w (x 1 ) p wg (x 1 ) -p wg (x 2 ) ≤ C 1 x 1-β 3 2 x 2 -x 1 x θw 2 -x θw 1 , (7.41)
where C is the above constant divided by α w . On the other hand, we use the lower bound (7.31) for the difference in f w and we need a lower bound for the difference in g. It is derived from (5.51),

g(x 2 ) -g(x 1 ) ≥ α 3 C max α w θ w α o θ o 1 4 θo+β 4 -1 x θw+β 3 -1 1 (x 2 -x 1 ) ≥ α 3 C max α w θ w α o θ o 1 4 θo+β 4 -1 1 2 θw+β 3 -1 x θw+β 3 -1 2 (x 2 -x 1 ).
(7.42) Hence (7.31) and (7.42) yield

f w (x 2 ) -f w (x 1 ) g(x 2 ) -g(x 1 ) ≥ Cx θw+β 3 -1 2 x θw 2 -x θw 1 (x 2 -x 1 ), (7.43)
with the product of the constants of (7.31) and (7.42). Then by combining (7.41) and (7.43), we deduce that

η w (x 2 ) -η w (x 1 ) p wg (x 1 ) -p wg (x 2 ) ≤ C h γ 1 θw i f w (x 2 ) -f w (x 1 ) g(x 2 ) -g(x 1 ) ;
which is (7.39) when δ 1 satisfies (7.40).

The case when 1 -x 1 is not too small is handled by the next proposition. and thus

η w (x 2 ) -η w (x 1 ) p wg (x 1 ) -p wg (x 2 ) ≤ C 1 (1 -x 1 ) θo-1 f w (x 2 ) -f w (x 1 ) g(x 2 ) -g(x 1 ) ≤ C 1 h γ 2 (θo-1) i f w (x 2 ) -f w (x 1 ) g(x 2 ) -g(x 1
) , whence (7.45).

In view of (7.32), (7.34), (7.38), and (7.40), let us choose

δ 1 = δ 1 = β 3 2θ w + β 3 , γ 1 = 2 2θ w + β 3 . ( 7
δ 2 = δ 2 = 2 + β 4 2θ o + β 4 , γ 2 = 2 2θ o + β 4 , (7.48)
satisfies (7.35), (7.37), (7.44), (7.46). Then the desired limit follows by collecting these results.

Lemma 7.7. Under the assumptions and notation on the mobility (5.42)-(5.46), the term T 1 defined in (7.20) tends to zero, with a similar limit in the non-wetting phase, i.e., lim (h,τ )→(0,0)

N n=1 τ M i,j=1 c ij S n,j S n,i f o (x) η w (S n,ij w ) -η w (x) p c (x) dx V n,j -V n,i = 0, lim (h,τ )→(0,0) N n=1 τ M i,j=1 c ij S n,j S n,i f w (x) η o (S n,ij o ) -η o (x) p c (x) dx V n,j -V n,i = 0.
(7.49)

Proof. We prove the first limit. Here the parameters of Propositions 7.4 and 7.5 are chosen by (7.47) and (7.48). It stems from the above considerations that, for each index n, the set of all indices (i, j) from 1 to M can be partitioned into three subsets,

O 1 = {(i, j) ; 0 ≤ S n,i < S n,j ≤ 3 4 }, O 2 = {(i, j) ; 1 4 ≤ S n,i < S n,j ≤ 1}, O 3 = {(i, j) ; 0 ≤ S n,i ≤ 1 4 and 3 4 ≤ S n,j ≤ 1}.
In turn, O 1 and O 2 can each be partitioned into two subsets

O 1,1 = {(i, j) ∈ O 1 ; S n,j ≤ h γ 1 T }, O 1,2 = {(i, j) ∈ O 1 ; S n,j > h γ 1 T }, O 2,1 = {(i, j) ∈ O 2 ; 1 -S n,i ≤ h γ 2 T }, O 2,2 = {(i, j) ∈ O 2 ; 1 -S n,i > h γ 2 T }.
To simplify, let

A i,j = c ij S n,j S n,i f o (x) η w (S n,ij w ) -η w (x) p c (x) dx V n,j -V n,i .
In view of (7.33) and (7.36), for all pairs (i, j) in O ,1 , = 1, 2, A i,j satisfies

|A i,j | ≤ C ∇ v L ∞ (Ω×]0,T [) h 2+δ i c ij .
Owing to (7.39) and (7.45), for all pairs (i, j) in O ,2 , = 1, 2, we have

|A i,j | ≤ C ∇ v L ∞ (Ω×]0,T [) h δ i c ij f w (S n,j ) -f w (S n,i ) g(S n,j ) -g(S n,i ) .
Finally, for all pairs (i, j) in O 3 ,

|A i,j | ≤ C ∇ v L ∞ (Ω×]0,T [) h i c ij f w (S n,j ) -f w (S n,i ) g(S n,j ) -g(S n,i ) .
According to (5.40), the sum of the terms over all (i, j) in O ,2 and O 3 tends to zero. For the remaining terms, observe that by definition,

h 2 i c ij ≤ C|∆ i ∩ ∆ j |,
so that the sum over all (i, j) in O ,1 is bounded by Ch δ i that also tends to zero, whence the first part of the limit (7.49). The same limit to zero holds for the non-wetting phase. With (7. [START_REF] Helmig | Multiphase flow and transport processes in the subsurface: a contribution to the modeling of hydrosystems[END_REF]), this lemma leads to the desired limit of the term with the auxiliary pressures. (7.51) 7.2. Convergence of the right-hand sides. In order to pass to the limit in the right-hand sides of (2.42)-(2.43) it is convenient to replace the quadrature formulas by integrals. Since the quadrature formulas are exact for polynomials of degree one, this is achieved by approximating some functions with the operator ρ h , see (2.27). As s in belongs to L ∞ (Ω×]0, T [), standard approximation properties of ρ τ and r h and a density argument imply lim Similarly, since q belongs to L 2 (Ω×]0, T [), lim (h,τ )→(0,0) ρ τ (ρ K (q)) = q in L 2 (K×]0, T [).

Also the (constant in space) correction added to ρ τ (r h (q)) satisfies lim (h,τ )→(0,0)

ρ τ 1
|Ω| Ω (r h (q) -q) = 0 in L 2 (Ω×]0, T [).

Therefore lim (h,τ )→(0,0)

qh,τ = q in L 2 (Ω×]0, T [).

(7.54)

With the same function V h,τ , consider the first term in the right-hand sides of (2.42)-(2.43)

X := N n=1 τ I h (f α (s n in,h ))q n h , V n h h = T 0 I h (f α (s in,h,τ ))q h,τ , V h,τ h .
By definition of the quadrature formula, X has the following expression:

X = N n=1 τ K∈ Ω |K| d + 1 d+1 =1 f α (s n, i in,h,τ )q n, i h,τ V n, i h,τ .
By inserting f α (ρ τ (ρ K (s in ))) and ρ τ (ρ K (q)), this becomes

X = N n=1 τ K∈ Ω |K| d + 1 d+1 =1 f α (s n, i in,h,τ ) -f α (ρ τ (ρ K (s in ))) qn, i h,τ V n, i h,τ + N n=1 τ K∈ Ω |K| d + 1 d+1 =1 f α (ρ τ (ρ K (s in ))) qn, i h,τ -ρ τ (ρ K (q)) V n, i h,τ + T 0 Ω f α (ρ τ (ρ K (s in )))ρ τ (ρ K (q))V h,τ = X 1 + X 2 + X 3 ,
since the last summand is a polynomial of degree one. We have lim (h,τ )→(0,0)

X 3 = T 0 Ω f o (s in ) q v.
It remains to show that X 1 and X 2 tend to zero. For X 1 , since f o and f w have the same derivative (up to the sign), we deduce from (5.52), (5.42), (5.43), (5.45), (5.46), and (5.47) that f α is bounded in [0, 1]; hence

|f α (s n, i in,h,τ ) -f α (ρ τ (ρ K (s in )))| ≤ C|s n, i in,h,τ -ρ τ (ρ K (s in ))|.
Thus, the summand is bounded by polynomials and the equivalence of norms yields

|X 1 | ≤ C v L ∞ (Ω×]0,T [) s in,h,τ -ρ τ (ρ K (s in )) L 2 (Ω×]0,T [) qh,τ L 2 (Ω×]0,T [) ,
that tends to zero with (h, τ ). It is easy to check that the same holds for X 2 . Hence lim (h,τ )→(0,0)

T 0 I h (f α (s in,h,τ ))q h,τ , V h,τ h = T 0 Ω f α (s in ) q v.
(7.55)

The argument for the second term in the right-hand side of (2.42) is much the same; we insert ρ τ (ρ K (s)) and we use the fact that lim I h (f α (s in,h,τ ))q h,τ -I h (f α (S h,τ ))q hτ , V h,τ h = T 0 Ω f α (s in ) q -f α (s) q v.

(7.57) 7.3. The full scheme. It remains to pass to the limit in the time derivative, say in (2.42), summed over n, and tested with the same V h,τ as previously, except that here we take v(T ) = 0. After summation by parts, this term reads (V n+1,i -V n,i )S n,i .

N n=1 (S n h -S n-1 h , V n h ) ϕ h = - N -1 n=1 (V n+1 h -V n h , S n h ) ϕ h -(V
By inserting ρ K (V n+1,i -V n,i ) in each element, this becomes

(V n+1 h -V n h , S n h ) ϕ h = (V n+1 h -V n h -ρ h (V n+1 -V n ), S n h ) ϕ h + Ω ϕρ h (V n+1 -V n )S n h .
The first term has the bound

(V n+1 h -V n h -ρ h (V n+1 -V n ), S n h ) ϕ h ≤ ϕ L ∞ (Ω) V n+1 h -V n h -ρ h (V n+1 -V n ) h S n h h .
Since the functions are piecewise polynomials, the equivalence of norms yields

N -1 n=1 (V n+1 h -V n h -ρ h (V n+1 -V n ), S n h ) ϕ h ≤ C ϕ L ∞ (Ω) × N -1 n=1 τ 1 τ I h (v n+1 -v n ) -ρ h (V n+1 -V n ) 2 L 2 (Ω) 1 2 N -1 n=1 τ S n h 2 L 2 (Ω) 1 2 .
Then the regularity of v, the approximation properties of I h and ρ h and the boundedness of S h,τ imply that lim (h,τ )→(0,0)

N -1 n=1 (V n+1 h -V n h -ρ h (V n+1 -V n ), S n h ) ϕ h = 0.
Similarly, it is easy to check from the convergence of S h,τ that -lim (h,τ )→(0,0)

N -1 n=1 Ω ϕρ h (V n+1 -V n )S n h = - T 0 Ω ϕ(∂ t v)s.
The treatment of the initial term is the same. Hence lim (h,τ )→(0,0)

N n=1 (S n h -S n-1 h , V n h ) ϕ h = - T 0 Ω ϕ(∂ t v)s - Ω ϕs 0 v. (7.59)
By combining (7.59), with (7.9) and (7.57), we readily see that the limit functions s, pα and p αg (s) satisfy the weak formulation (1.16). This proves Theorem 2.6.

8. Numerical validation. This section proposes a numerical validation of our algorithm with a two dimensional finite difference code. Details on the algorithm implemented are given. A problem with manufactured solutions is then considered to study the convergence properties of our algorithm.

We note that to facilitate the implementation of this algorithm in a two dimensional finite difference code, the source terms of the equations (2.32)-(2.33) have been replaced by functions denoted by f 1 and f 2 .

8.2. Numerical test with a manufactured solution. The numerical validation of the algorithm is done by approximating the analytical solutions defined by P w (t, x, y) = 2 + x 2 y -y 2 + x 2 sin(t + y), ( The mobilities η w and η o , introduced in section 1.1, are defined as follows:

η w = 4S 2 , η o = 0.4(1 -S) 2 . (8.7)
The capillary pressure is based on the Brooks-Corey model, it reads:

P c = AS -0.5
if S > 0.05, A(1.5 -10S) × 0.05 -0.5 otherwise. (8.8) where A is a constant set to 50. The term sources f 1 and f 2 are computed accordingly. The convergence tests are performed on a set of six uniform grids with respective mesh size h ∈ {0.2, 0.1, 0.05, 0.025, 0.0125, 0.00625}. The convergence properties are evaluated by using a time step τ set to the mesh size h with a final time T = 1. As the time derivatives and the saturations S n+1,ij w , S n+1,ij o are computed with first order time approximation, we expect the convergence rate in the L 2 norm to be of order one. The results of the convergence tests are presented in Table 8.1. The theoretical order of convergence, equal to one, is recovered for both unknowns which confirms the correct behavior of the algorithm. 9. Conclusions. This paper formulates a P 1 finite element method to solve the immiscible two-phase flow problem in porous media. The unknowns are the phase pressure and saturation, which are the preferred unknowns in industrial reservoir simulators. The numerical method employs mass lumping for integration and an upwind flux technique. As a consequence, the saturation is shown to be bounded between zero and one. The discrete approximations of pressure and saturation converge to the weak solution as the time step and mesh sizes tend to zero.

L 2 -

1 x

 1 p c (s)ds. (1.11) Since p c is a continuous function on [0, 1], the function g c belongs to C 1 ([0, 1]). Next, we introduce the auxiliary pressures p wg , p wo , and g, ∀x ∈ [0, 1], p wg (x) = x 0 f o (s)p c (s)ds, p og (x) = x 0 f w (s)p c (s)ds, (1.12) ∀x ∈ [0, 1], g(x) = -x 0 η w (s)η o (s) η w (s) + η o (s) p c (s)ds. (1.13) Owing to (1.6), ∀x ∈ [0, 1], p wg (x) + p og (x) = x 0 p c (s)ds = p c (x) -p c (0). (1.14)

6. 2 .. 22 ) 6 . 4 .

 22264 Weak convergence. All constants below are independent of h and τ .The bound (3.1) on the discrete saturation S h,τ implies that there exists a function s ∈ L ∞ (Ω×]0, T [) and a subsequence of (h, τ ) not indicated, such that lim (h,τ )→(0,0) S h,τ = s weakly* in L ∞ (Ω×]0, T [). (6Proposition The limit function s satisfies ∀(x, t) a.e. in Ω×]0, T [, 0 ≤ s(x, t) ≤ 1.(6.23)

Theorem 7. 8 . 0 P

 80 Let v ∈ C 1 ( Ω × [0, T ]) be a smooth function and let V h,τ (t) = I h (v)(t n ) in ]t n-1 , t n ].Under the assumptions and notation on the mobility (5.42)-(5.46), lim (h,τ )→(0,0)T 0 P α,h,τ , I h (η α (S h,τ )); I h (p αg (S h,τ )), V h,τ h = T 0 Ω ∇ g(s) • ∇ v, α = w, o,(7.50) where s is the limit of S h,τ .Finally, Theorems 7.1 and 7.8, together with (1.15) and (6.40), give the desired convergence of the upwind diffusion terms.Theorem 7.9. With the notation and assumptions of Theorem 7.1, we have for allfunctions v ∈ C 1 ( Ω × [0, T ]), lim (h,τ )→(0,0) -T α,h,τ , I h (η α (S h,τ )); P α,h,τ , V h,τ h = T 0 Ω η w (s)∇(p w + p wg (s)) + ∇ g(s) • ∇ v if α = w, = T 0 Ω η o (s)∇(p o -p og (s)) -∇ g(s) • ∇ v if α = o.

(

  h,τ )→(0,0)ρ τ (ρ K (s in )) = s in in L ∞ (K×]0, T [).(7.52)Then the continuity of f α , for α = w, o, yields lim(h,τ )→(0,0) f α (ρ τ (ρ K (s in ))) = f α (s in ) in L ∞ (K×]0, T [).(7.53)

( 0 (

 0 h,τ )→(0,0) S h,τ -ρ τ (ρ K (s)) L 2 (Ω×]0,T [) = 0.Then the argument used for the first term readily gives lim (h,τ )→(0,0)T I h (f α (S h,τ ))q hτ , V h,τ h = T 0 Ω f α (s) q v. (7.56)By combining (7.55) and (7.56), we obtain convergence of the right-hand sides, lim (h,τ )→(0,0) T 0

8. 1 .)m i f n+1,i 1 , 1 ≤)c) = m i f n+1,i 2 , 1 ≤

 11121 Implementation of the model. To avoid dealing with nonlinear terms, we implement a modified version of the algorithm proposed in section 2.3. The main difference consists of approximating the terms S n+1,ij w , S n+1,ij o and P n+1 o with first time order extrapolation. For each node 1 ≤ i ≤ M , the unknowns (S n+1,i , P n+1,i w ) are computed as the solution of the following problem: mi ∆t (S n+1,i -S n,i ) -j =i,j∈N (i) c ij η w (S * ,n+1,ij w i ≤ M, ,i -S n,i )j =i,j∈N (i) c ij η o (S * ,n+1,ij o ij η o (S * ,n+1,ij o )(p * ,n+1,j c -p * ,n+1,i c i ≤ M,(8.2) where the pressure P o has been substituted with P w + p c with respect to (2.34). The solution P w is enforced to satisfy (2.35) a posteriori by subtracting its integral M i=1 m i P i w after solving the above problem. The terms S * ,n+1,ij w and S * ,n+1,ij o are approximated at time iteration n by setting them to S n,ij w and S n,ij o . Eventually, the capillary pressure p * ,n+1 cis approximated with a first order Taylor expansion with respect to the saturation S, it reads: -S n ).

  ) S(t, x, y) = 0.2(2 + 2xy + cos(t + x)),(8.5) on the computational domain Ω = [0, 1] 2 . Dirichlet boundary conditions are applied on ∂Ω on both unknowns P w and S. The initial conditions of the problem satisfy (8.4)-(8.5). The porosity of the domain is set to: φ(t, x, y) = 0.2(1 + xy). (8.6)

  .34)

		S n,ij o	=	 		S n,i if P n,i o > P n,j o S n,j if P n,i o < P n,j o	(2.37)
					min(S n,i , S n,j ) if P n,i o = P n,j o
	We observe that
				S n,ij w	= S n,ji w , S n,ij o	= S n,ji o ,
	so that, if we interpret in (2.32) (respectively, (2.33)) η w (S n,ij w ) (respectively, η o (S n,ij o )) as W i,j , then (2.20) and hence (2.24) hold.
	Remark 2.5. Before setting (2.32)-(2.35) in variational form, observe that:
	1. The scheme (2.32)-(2.35) forms a square system in the primary unknowns, S n h and P n w .
	2. Formula (2.32) is also valid for i = M . Indeed, we pass to the left-hand side the right-
	hand side of (2.32) and set A i the resulting line of index i. Let ÃM denote what should be
	the line of index M , i.e.,
	ÃM =	mM (ϕ)
						M
						m i P n,i w = 0.	(2.35)
						i=1
	Here i runs from 1 to M -1 in (2.32) and from 1 to M in (2.33); the upwind values S n,ij w , S n,ij o
	are defined by
		S n,ij w	=	 		S n,i if P n,i w > P n,j w S n,j if P n,i w < P n,j w	(2.36)
					max(S n,i , S n,j ) if P n,i w = P n,j w

  3.2. First pressure bounds.The following properties will be used frequently: Lemma 3.2. The fact that p c is strictly decreasing and (2.34) yield the following:

	P i w > P j w , and P i o ≤ P j o	implies S i ≥ S j ,	(3.10)
	If P i w = P j w , then P i o ≥ P j o if and only if S i ≤ S j .	(3.11)
	If P i o = P j o , then P i w ≤ P j w , if and only if S i ≤ S j .	(3.12)

  norm of error Water pressure P w Water saturation S

	h	n df	Error	Rate	Error	Rate
	0.2	25	8.50E-3	-	4.21E-3	-
	0.1	100	4.15E-3	1.03	2.30E-3	0.87
	0.05	400	2.08E-3	1.00	1.14E-4	1.01
	0.025	1600	1.04-3	1.00	5.57E-4	1.03
	0.0125	6400	5.23E-4	0.99	2.75E-4	1.02

Table 8 . 1

 81 Results of convergence tests where the mesh size is denoted by h and the number of degrees of freedom per unknown by n df . The time step τ is set to h and errors are computed at final time T = 1.

j =i,j∈N (i)

ported in part by NSF-DMS 1 Proposition 2.1. Under condition

where 0 < δ 1 ≤ γ 1 (θ w + β 3 ) -1.

(7.34)

Similarly, suppose that 1 -x 2 < 1 -x 1 ≤ h γ 2 i for some exponent γ 2 > 0 such that

.35)

Then h i η w (x 2 ) -η w (x 1 ) p wg (x 1 ) -p wg (x 2 ) ≤ C h 2 i h δ 2 i , (7.36)

In both cases, the constants C are independent of x 1 , x 2 , and h i .

Proof. In the first case, according to (7.25), the choice (7.34) and (7.32) on γ 1 , we have

with the constant C of (7.25), which gives (7.33). In the second case, the same argument leads to

, with the constant C of (7.27), thus implying (7.36) with the choice (7.37) for δ 2 and the condition (7.35) on γ 2 . Now, we turn to the case when x 2 is not too small.

Proposition 7.5. In addition to (7.32), suppose that the exponent γ 1 of Proposition 7.4 satisfies

) , (7.39)

Again, the constant C is independent of x 1 , x 2 , and h i .

x 2 , and we examine each case. 1) When x 1 ≤ 1 2 x 2 , formula (7.26) leads to

)

Proposition 7.6. In addition to (7.35), suppose that the exponent γ 2 of Proposition 7.4 satisfies

Then

) , (7.45)

Again, the constant C is independent of x 1 , x 2 , and h i .

Proof. The proof is analogous to that of Proposition 7.5, but we sketch the steps for the reader's convenience. We skip the constants' details, but stress that they are independent of x 1 , x 2 , and h i . Again, there are two possibilities, either 1 -

, and we examine each case. 1) In the first case,

Hence

With (7.44) and (7.46), this implies (7.45).

2) In the second case, we have on the one hand,

so that η w (x 2 ) -η w (x 1 ) p wg (x 1 ) -p wg (x 2 ) ≤ C(x 2 -x 1 ) 2 (1 -x 1 ) θo+β 4 -1 .

On the other hand, f w (x 2 ) -f w (x 1 ) ≥ C(x 2 -x 1 )(1 -x 1 ) θo-1 , and g(x 2 ) -g(x 1 ) ≥ C(x 2 -x 1 )(1 -x 1 ) θo+β 4 -1 ,