Upstream mobility Finite Volumes for the Richards equation in heterogenous domains

This paper is concerned with the Richards equation in a heterogeneous domain, each subdomain of which is homogeneous and represents a rocktype. Our first contribution is to rigorously prove convergence toward a weak solution of cell-centered finite-volume schemes with upstream mobility and without Kirchhoff's transform. Our second contribution is to numerically demonstrate the relevance of locally refining the grid at the interface between subregions, where discontinuities occur, in order to preserve an acceptable accuracy for the results computed with the schemes under consideration.

The conditions to be satisfied by φ, λ, η and S will be elaborated on later. In a homogeneous medium, these physical properties are uniform over Ω, i.e., φpxq " φ 0 , λpxq " λ 0 , ηps, xq " η 0 psq, Spp, xq " S 0 ppq for all x P Ω. In a heterogeneous medium, the dependence of φ, λ, η and S on x must naturally be taken into account. The quantity s, called saturation, measures the relative volumic presence of water in the medium. The quantity p is the pressure. Let T ą 0 be a finite time horizon. We designate by Q T " p0, T qˆΩ the space-time domain of interest. Our task is to find the saturation field s : Q T Ñ r0, 1s and the pressure field p : Q T Ñ R so as to satisfy • the interior equations φpxq B t s`div F " 0 in Q T , (1.1a) F`λpxq ηps, xq ∇pp´ g¨xq " 0 in Q T , (1.1b) s´Spp, xq " 0 in Q T ; (1.1c) • the boundary conditions F¨npxq " 0 on p0, T qˆΓ N , (1.1d) ppt, xq " p D pxq on p0, T qˆΓ D ; (1.1e) • the initial data sp0, xq " s 0 pxq in Ω. (1.1f) The partial differential equation (1.1a) expresses the water volume balance. The flux F involved in this balance is given by the Darcy-Muskat law (1.1b), in which g is the gravity vector and is the known constant density of water, assumed to be incompressible. It is convenient to introduce ψ "´ g¨x, ϑ " p`ψ, (1.2) referred to respectively as gravity potential and hydraulic head. In this way, the Darcy-Muskat law (1.1b) can be rewritten as F`λpxq ηps, xq ∇pp`ψq " F`λpxq ηps, xq ∇ϑ " 0.
Equation (1.1c) connecting the saturation s and the pressure p is the capillary pressure relation. The boundary BΩ is split into two non-overlapping parts, viz., where Γ N is open and Γ D is closed, the latter having a positive pd´1q-dimensional Hausdorff measure ν d´1 pΓ D q ą 0. The no-flux Neumann condition (1.1d) is prescribed on p0, T qˆΓ N , where npxq is the outward normal unit vector at x P Γ N . The Dirichlet condition (1.1e) with a known Lipschitz function p D P W 1,8 pΩq is imposed on p0, T qˆΓ D . Note that, in our theoretical development, the function p D is assumed to be defined over the whole domain Ω, which is stronger than a data p D P L 8 pΓ D q given only on the boundary. The assumption that p D does not depend on time can be removed by following the lines of [14], but we prefer here not to deal with timedependent boundary data in order to keep the presentation as simple as possible. Finally, the initial data s 0 P L 8 pΩ; r0, 1sq in (1.1f) is also a given data.
In this work, we restrict ourselves to a specific type of heterogeneous media, defined as follows. We assume that the domain Ω can be partitioned into several connected polyhedral subdomains Ω i , 1 ď i ď I. Technically, this means that if Γ i,j denotes the interface between Ω i and Ω j (which can be empty for some particular choices of ti, ju), then (1.4) with Γ " Ť i‰j Γ i,j . Each of these subdomains corresponds to a distinctive rocktype. Inside each Ω i , the physical properties are homogeneous. In other words, λpxq " λ i , ηps, xq " η i psq, Spp, xq " S i ppq for all x P Ω i . Therefore, system (1.1) is associated with φpxq " where 1 Ωi stands for the characteristic function of Ω i . For all i P t1, . . . , Iu, we assume that φ i P p0, 1s and λ i ą 0. Furthermore, we require that η i is increasing on r0, 1s, η i p0q " 0, where µ ą 0 is the (known) viscosity of water. In addition to the assumption that Sp¨, xq, defined in (1.5b), is absolutley continuous and nondecreasing, the functions S i are also subject to some generic requirements commonly verified the models available in the literature: for each i P t1, . . . , Iu, there exists p i ď 0 such that S i is increasing on p´8, p i s, lim pÑ´8 S i ppq " 0, S i " 1 on rp i ,`8q.
(1.6b) This allows us to define an inverse S´1 i : p0, 1s Ñ p´8, p i s such that S i˝S´1 i psq " s for all s P p0, 1s. We further assume that for all i P t1, . . . , Iu the function S i is bounded in L 1 pR´q, or equivalently, that S´1 i P L 1 p0, 1q. It thus makes sense to consider the capillary energy density functions e i : RˆΩ i Ñ R`defined by e i ps, xq " For all x P Ω i , the function e i p¨, xq is nonnegative, convex since S´1 i is monotone, and bounded on r0, 1s as a consequence of the integrability of S i . For technical reasons that will appear clearly later on, we further assume that a η i˝Si P L 1 pR´q, @i P t1, . . . , Iu. At the interface Γ i,j between Ω i and Ω j , i ‰ j, any solution of (1.1a)-(1.1c) satisfies the matching conditions F i¨ni`Fj¨nj " 0 on p0, T qˆΓ i,j , (1.10a) In the continuity of the normal fluxes (1.10a), which is enforced by the conservation of water volume, n i denotes the outward normal to BΩ i and F i¨ni stands for the trace of the normal component of F | Q i,T on p0, T qˆBΩ i . In the continuity of pressure (1.10b), which also results from (1.1a)-(1.1c), p i denotes the trace on p0, T qˆBΩ i of the pressure p | Q i,T in the i-th domain.

Stability features and notion of weak solutions
We wish to give a proper sense to the notion of weak solution for problem (1.1). To achieve this purpose, we need a few mathematical transformations the definition of which crucially relies on a fundamental energy estimate at the continuous level. The calculations below are aimed at highlighting this energy estimate and will be carried out in a formal way, in constrast to those in the fully discrete setting. Multiplying (1.9a) by p´p D , invoking (1.7), integrating over Ω i and summing over i, we end up with d dt We now integrate by parts the second term. Thanks to the matching conditions (1.10) and the regularity of p D , we obtain A :" It follows from the flux value (1.9b) that Young's inequality, combined with the boundedness of ∇p D , ∇ψ, λ and η, yields for some C ě 0 depending only on λ, η, ψ, µ, Ω and p D . Let us define the energy E : r0, T s Ñ R`by Eptq " Integrating (1.11) w.r.t. time results in EpT q`1 2 Estimate (1.12) is the core of our analysis. However, it is difficult to use in its present form since η i psq " η i pS i ppqq vanishes as p tends to´8, so that the control of ∇p degenerates. To circumvent this difficulty, we resort to the nonlinear functions (customarily referred to as the Kirchhoff transforms) Θ i : R Ñ R, Φ i : R Ñ R, and Υ : RˆΩ Ñ R respectively defined by the notion of Υ being due to [23]. Bearing in mind that EpT q ě 0, estimate (1.12) implies that (1.14) As Φ i˝Θ´1 i is Lipschitz continuous, this also gives rise to a L 2 pQ i,T q-estimate on ∇Φ i ppq. The functions ř i Θ i ppq1 Ωi and ř i Φ i ppq1 Ωi are in general discontinuous across the interfaces Γ i,j , unlike Υppq. Since the functions Υ˝Θ´1 i are Lipschitz continuous, we can readily infer from (1.14) that for some C depending on T , Ω, ∇p D 8 , the S i L 1 pR´q 's and the last equality being due to (1.6a). Moreover, Υppq´Υpp D q vanishes on p0, T qˆΓ D . Poincaré's inequality provides a L 2 pQ T qestimate on Υppq since Γ D has positive measure and since Υpp D q is bounded in Ω. In view of assumption (1.8), the functions Θ i and Υ are bounded on R´. Besides, for p ě 0, η i˝Si ppq " 1{µ, so that Θ i ppq " p a λ i {µ and Υppq " min 1ďiďI p a λ i {µ. It finally comes that from which we infer a L 2 pQ i,T q-estimate on Θ i ppq. Putting V " u P H 1 pΩq | u | Γ D " 0 ( , the above estimates suggest the following notion of weak solution for our problem. The expression (1.17d) is a reformulation of the original one (1.9b) in a quasilinear form which is suitable for analysis, even though the physical meaning of the Kirchhoff transform Φ i ppq is unclear. While the formulation (1.17c) should be thought of as a weak form of (1.9a), (1.10a), (1.1f), and (1.1d), the condition Υppq´Υpp D q P L 2 pp0, T q; V q contains (1.10b) and (1.1e).

Goal and positioning of the paper
We are now in a position to clearly state the two objectives of this paper.
The first objective is to put forward a rigorous proof that, for problem (1.1) with heterogeneous data (1.5), cell-centered finite-volume schemes with upstream mobility such as described in §2.2, do converge towards a weak solution (in the sense of Definition 1.1) as the discretization parameters tend to 0. Such mathematically assessed convergence results are often dedicated to homogeneous cases: see for instance [3,30,40] for schemes involving the Kirchhoff transforms for Richards' equation, [1] for a upstream mobility CVFE approximation of Richards' equation in anisotropic domains, [17][18][19] for schemes for two-phase flows involving the Kirchhoff transform, and [31,34] for upstream mobility schemes for two-phase porous media flows. For flows in highly heterogeneous porous media, rigorous mathematical results have been obtained for schemes involving the introduction of additional interface unknowns and Kirchhoff's transforms (see for instance [6,11,12,23]), or under the non-physical assumption that the mobilities are strictly positive [28,29]. It was established very recently in [7] that cell-centered finite-volumes with (hybrid) upwinding also converge for two-phase flows in heterogeneous domains, but with a specific treatment of the interfaces located at the heterogeneities. Here, the novelty lies in the fact that we do not consider any specific treatment of the interface in the design of the scheme.
The second objective is of more practical nature. Even though our analysis still holds without any specific treatment of the interface, it is well-known that cell-centered upstream mobility finitevolumes can be inaccurate in the presence of heterogeneities. This observation motivated several contributions (see for instance [23,24,29,35]) where skeletal (i.e., edge or vertex) unknowns where introduced in order to enforce the continuity of the pressures at the interfaces Γ i,j . By means of extensive numerical simulations in §6, we will show that without local refinement of the grid at the interface, the method still converges, but with a degraded order. Our ultimate motivation is to propose an approach which consists in adding very thin cells on both sides of the interface before using the cell centered scheme under study. Then the scheme appears to behave better, with first-order accuracy. Moreover, one can still make use of the parametrized cut-Newton method proposed in [4] to compute the solution to the nonlinear system corresponding to the scheme. This method appears to be very efficient, while it avoids the possibly difficult construction of compatible parametrizations at the interfaces as in [7][8][9].

Finite-volume discretization
The scheme we consider in this paper is based on two-point flux approximation (TPFA) finitevolumes. Hence, it is subject to some restrictions on the mesh [26,33]. We first review the requirements on the mesh in §2.1. Next, we construct the upstream mobility finite-volume scheme for Richards' equation in §2.2. The main mathematical results of the paper, which are the wellposedness of the nonlinear system corresponding to the scheme and the convergence of the scheme, are then summarized in §2.3.

Admissible discretization of Q T
Let us start by discretizing w.r.t. space. Definition 2.1. An admissible mesh of Ω is a triplet pT , E , px K q KPT q such that the following conditions are fulfilled: (i) Each control volume (or cell) K P T is non-empty, open, polyhedral and convex, with positive d-dimensional Lebesgue measure m K ą 0. We assume that Moreover, we assume that the mesh is adapted to the heterogeneities of Ω, in the sense that for all K P T , there exists i P t1, . . . , Iu such that K Ă Ω i . (ii) Each face σ P E is closed and is contained in a hyperplane of R d , with positive pd´1qdimensional Hausdorff measure ν d´1 pσq " m σ ą 0. We assume that ν d´1 pσ X σ 1 q " 0 for σ, σ 1 P E unless σ 1 " σ. For all K P T , we assume that there exists a subset E K of E such that BK " Ť σPE K σ. Moreover, we suppose that Given two distinct control volumes K, L P T , the intersection K X L either reduces to a single face σ P E denoted by K|L, or its pd´1q-dimensional Hausdorff measure is 0.
(iii) The cell-centers px K q KPT are pairwise distinct with x K P K, and are such that, if K, L P T share a face K|L, then the vector x L´xK is orthogonal to K|L. (iv) For the boundary faces σ Ă BΩ, we assume that either σ Ă Γ D or σ Ă Γ N . For σ Ă BΩ with σ P E K for some K P T , we assume additionally that there exists x σ P σ such that x σ´xK is orthogonal to σ.
In our problem, the standard Definition 2.1 must be supplemented by a compatibility property between the mesh and the subdomains. By "compatbility" we mean that each cell must lie entirely inside a single subregion. Put another way, This has two consequences. The first one is that, if we define The second one is that the subdomain interfaces Γ i,j for i ‰ j coincide necessarily with some edges σ P E . To express this more accurately, let E Γ " tσ P E | σ Ă Γu be the set of the interface edges, E D ext " tσ P E | σ Ă Γ D u be the set of Dirichlet boundary edges, and E N ext " tσ P E | σ Ă Γ N u be the set of Neumann boundary edges. Then, Γ " For later use, it is also convenient to introduce the subset E i Ă E consisting of those edges that correspond to cells in T i only, i.e., and the subset E int of the internal edges, i.e., To each edge σ P E , we associate a distance d σ by setting We also define d Kσ " distpx K , σq for all K P T and σ P E K . The transmissivity of the edge σ P E is defined by a σ " m σ d σ . (2.5) Throughout the paper, many discrete quantities u will be defined either in cells K P T or on Dirichlet boundary edges σ P E D ext , i.e. u " ppu K q KPT , pu σ q σPE D ext q P X T YE D ext , where X can be either R , ě 1, or a space of functions. Then for all K P T and σ P E K , we define the mirror value u Kσ by (2.6) The diamond cell ∆ σ corresponding to the edge σ is defined as the convex hull of tx K , x Kσ , σu for K such that σ P E K , while the half-diamond cell ∆ Kσ is defined as the convex hull of tx K , σu. Denoting by m ∆σ the Lebesgue measure of ∆ σ , the elementary geometrical relation m ∆σ " d m σ d σ where d stands for the dimension will be used many times in what follows.
Another notational shorthand is worth introducing now, since it will come in handy in the sequel. Let f p¨, xq " be a scalar quantity or a function whose dependence of x P Ω is of the type (1.5). Then, for K P T , we slightly abuse the notations in writing where the index ipKq is defined in (2.1). The last equality in the above equation holds by virtue of the compatibility property. For example, we will have not only φ K " φpx K q, λ K " λpx K q, η K psq " ηps, x K q, S K ppq " Spp, x K q but also e K psq " eps, x K q. Likewise, we shall be writing f Kσ p¨q " f p¨, x Kσ q for the mirror cell without any ambiguity: if σ P E int Y E N ext , then x Kσ is a cell-center; if σ P E D ext , then x Kσ lies on the boundary but does not belong to an interface between subdomains.
The size h T and the regularity ζ T of the mesh are respectively defined by The time discretization is given by pt n q 0ď1ďN with 0 " t 0 ă t 1 ă¨¨¨ă t N " T . We denote by ∆t n " t n´tn´1 for all n P t1, . . . , N u and by ∆t " p∆t n q 1ďnďN .

Upstream mobility TPFA Finite Volume scheme
Given a discrete saturation profile ps n´1 K q KPT P r0, 1s T at time t n´1 , n P t1, . . . , N u, we seek for a discrete pressure profile pp n K q KPT P R T at time t n solution to the following nonlinear system of equations. Taking advantage of the notational shorthand (2.7b), we define The volume balance (1.9a) is then discretized into using the approximation where the mirror values p n Kσ and ψ Kσ are given by (2.6). In the numerical flux (2.11a), the edge permeabilities pλ σ q σPE are set to while the edge mobilities are upwinded according to (2.11c) In practice, the definition of η n σ when ϑ n K " ϑ n Kσ has no influence on the scheme. We choose here to give a symmetric definition that does not depend on the orientation of the edge σ in order to avoid ambiguities.
The boundary condition p D is discretized into whereas the initial condition is discretized into The Dirichlet boundary condition is encoded in the fluxes (2.11a) by setting Bearing in mind the definition (2.6) of the mirror values for σ P E N ext , the no-flux boundary condition across σ P E N ext is automatically encoded, i.e., F n Kσ " 0 for all σ P E K X E N ext , K P T and n ě 1. In what follows, we denote by p n " pp n K q KPT for 1 ď n ď N , and by s n " ps n K q KPT for 0 ď n ď N . Besides, we set p D " ppp D K q KPT , pp D σ q σPE D q.

Main results and organization of the paper
The theoretical part of this paper includes two main results. The first one, which emerges from the analysis at fixed grid, states that the schemes admits a unique solution pp n q 1ďnďN .
With Theorem 2.2 at hand, we define the approximate pressure p T ,∆t by We also define the approximate saturation as The second main result guarantees the convergence towards a weak solution of the sequence of approximate solutions as the mesh size and the time steps tend to 0. Let pT m , E m , px K q KPTm q mě1 be a sequence of admissible discretizations of the domain Ω in the sense of Definition 2.1 such that The rest of this paper is outlined as follows. Section §3 is devoted to the numerical analysis at fixed grid. This encompasses the existence and uniqueness result stated in Theorem 2.2 as well as a priori estimates that will help proving Theorem 2.3. The convergence of the scheme, which is taken up in §4, relies on compactness arguments, which require a priori estimates that are uniform w.r.t. the grid. These estimates are mainly adaptations to the discrete setting of their continuous counterparts that arised in the stability analysis sketched out in §1.2. These estimates are shown in §4.1 to provide some compactness on the sequence of approximate solutions. In §4.2, we show that these compactness properties together with the a priori estimates are sufficient to identify any limit of an approximate solution as a weak solution to the problem.
In §5, we provide some details about the practical numerical resolution by laying emphasis on the switch of variable for selecting the primary unknown and on the mesh refinement at an interface in order to better enforce pressure continuity. Finally, in §6, numerical experiments on two configurations (drying and filling cases) for two capillary pressure models (Brooks-Corey and van Genuchten-Mualem) testify to the relevance of the local refinement strategy as a simple technique to preserve accuracy.
Remark 2.4. Theorem 2.3 only states the convergence of the scheme up to a subsequence. In the case where the weak solution is unique, then the whole sequence of approximate solutions would converge towards this solution. As far as we know, uniqueness of the weak solutions to Richards' equation is in general an open problem for heterogeneous media where x Þ Ñ Spp, xq is discontinuous. Uniqueness results are however available in the one-dimensional setting for a slightly more restrictive notion of solutions, cf. [12], or under additional assumptions on the nonlinearities η i , S i , cf. [11].
3 Analysis at fixed grid 3

.1 Some uniform a priori estimates
In this section, our aim is to derive a priori estimates on the solutions to the scheme (2.9)-(2.13). These estimates will be at the core of the existence proof of a solution to the scheme. They will also play a key role in proving the convergence of the scheme.
The main estimate on which our analysis relies is a discrete counterpart of (1.12). We recall that a σ is the transmissivity introduced in (2.5).
In (3.1), the relationship between σ and K is to be understood as follows. For an inner edge σ P E int , although it can be written as σ " K|L or L|K, only one of these contributes to the sum. For a boundary edge σ P E ext , there is only one cell K such that σ P E K , so there is no ambiguity in the sum.
Proof. Multiplying (2.10) by ∆t n pp n K´p D K q, summing over K P T and n P t1, . . . , N u, and carrying out discrete integration by parts yield The discrete energy density function e K : r0, 1s Ñ R`, defined by means of the notation (2.7) from the functions f i " e i introduced in (1.7), is convex by construction. Consequently, Therefore, the quantity A of (3.3a) can be bounded below by the last inequality being a consequence of the boundedness of e K on r0, 1s. Writing ϑ " p`ψ and expanding each summand of (3.3b), we can split B into It follows from [27, Lemma 9.4] and from the boundedness of η that there exists a constant C depending only on λ, µ, ζ T and Ω such that Thanks to these estimates and to the Cauchy-Schwarz inequality, we have On the other hand, Young's inequality provides Hence, So far, we have not used the upwind choice (2.11c) for the mobilities η n σ . This will be done in the next lemma, where we derive a more useful variant of estimate (3.1a), in which η n σ is replaced by η n σ defined below. In a homogeneous medium, η n σ ě η n σ so that the new estimate (3.8) seems to be stronger than (3.1a).
We begin by introducing the functions q η σ : R Ñ p0, 1{µs defined for σ P E by By virtue of assumptions (1.6), each argument of the minimum function is nondecreasing and positive function of p P R. As a result, q η σ is also a nondecreasing and positive function of p P R. Note that q η σ " η i˝Si for all σ P E i , while for interface edges σ Ă Γ i,j , the mere inequality q η σ ď η i˝Si holds. Next, we consider the intervals with the notations aKb " minpa, bq and aJb " maxpa, bq. At last, we set There exists a constant C 3 depending on the same data as C 1 such that Proof. We partition the set E of edges into three subsets, namely, Invoking q η σ " minpη K˝SK , η Kσ˝SKσ q, we can minorize the left-hand side of (3.1a) to obtain Starting from this inequality and using the boundedness of η i and ψ, we can readily show that there exists a constant C depending on the same data as C 1 such that in which the sum over E n 0 was omitted because all of its summands vanish. Simlarly to what was pointed out in equation 2.9 in [1], we notice that since η σ is nondecreasing w.r.t. p, it is straightforward to check that the definition q η n σ :" exactly amounts to (3.10) Taking advantage of this equivalence, we can transform D 1 into Estimate (3.8) finally follows from the boundedness of η, 1{λ and ψ.
The above lemma has several important consequences for the analysis. Let us start with discrete counterparts to estimations (1.14) and (1.15).
Moreover, there exists two constants C 4 , C 5 depending on the same data as C 1 and additionnally on ? Proof. Consider those edges σ P E i -defined in (2.3a)-corresponding to some fixed i P t1, . . . , Iu, for which q η σ " η i˝Si " |Θ 1 i | 2 and η n σ " max J n σ |Θ 1 i | 2 due to (1.13a). By summing the elementary inequality pΘ i pp n K q´Θ i pp n Kσ qq 2 ď η n σ pp n K´p n Kσ q 2 , over σ P E i , i P t1, . . . , Iu and n P t1, . . . , N u using appropriate weights, we get whose right-hand side is obviously less than C 3 , thanks to (3.8). This proves (3.12a). Similarly, the respective definitions of η n σ and Υ have been tailored so that max J n σ |Υ 1 | 2 ď η n σ for all σ P E . As a consequence, pΥpp n K q´Υpp n Kσ qq 2 ď η n σ pp n K´p n Kσ q 2 .
Summing these inequalities over σ P E and n P t1, . . . , N u with appropriate weights and invoking (3.8), we prove (3.12b). The argument for (3.13a) is subtler. Starting from the basic inequality we apply the discrete Poincaré inequality of [27, Lemma 9.1] -which is legitimate since Γ D has positive measure-followed by [27,Lemma 9.4] to obtain The purpose of the next lemma is to work out a weak estimate on the discrete counterpart of B t s, which will lead to compactness properties in §4.1. For ϕ P C 8 c pQ T q, let Lemma 3.4. There exists a constant C 6 depending on the same data as C 1 such that (3.14) Proof. Multiplying (2.10) by ∆t n ϕ n K , summing over K P T and n P t1,¨¨¨, N u and carrying out discrete integration by parts, we end up with Applying the Cauchy-Schwarz inequality and using (3.1b), we get

Existence of a solution to the scheme
The statements of the previous section are all uniform w.r.t. the mesh and are meant to help us passing to the limit in the next section. In contrast, the next lemma provides a bound on the pressure that depends on the mesh size and on the time-step. This property is needed in the process of ensuring the existence of a solution to the numerical scheme.
Lemma 3.5. There exist two constants C 7 , C 8 depending on T , ∆t n as well as on the data of the continuous model λ, µ, p D , ψ, ζ, Ω, T , φ, S i L 1 pR´q and ? η i˝Si L 1 pR´q , 1 ď i ď I, such that Proof. From (3.13a) and from Υppq " p a min i λ i {µ for p ě 0, we deduce that Hence, the upper-bound C 8 is found by maximizing the right-hand side over K P T and n P t1, . . . , N u.
To show that p n K is bounded from below, we employ a strategy that was developed in [13] and extended to the case of Richards' equation in [1,Lemma 3.10]. From (2.12), (2.14) and the boundedness of p D , it is easy to see that p n σ ě inf xPBΩ p D pxq, @σ P E D ext .
Estimate (3.8) then shows that for all K P T such that E K X E D ext ‰ H, we have p n K ě p n σ´d C 3 ∆t n a σ q η σ pp n σ q ": π n K , @σ P E K X E D ext .
The quantity π n K is well-defined, since q η σ pp n σ q ą 0 for p n σ ą´8, and does not depend on time, as p D does not either. Furthermore, if p n K is bounded from below by some π K , then the pressure in all its neighboring cells L P T such that σ " K|L P E K is bounded from below by p n L ě π n K´d C 3 ∆t n a σ q η σ pπ n K q ": π n L .
Again, π n L is well-defined owing to q η σ pπ n K q ą 0. Since the mesh is finite and since the domain is connected, only a finite number of edge-crossings is required to create a path from a Dirichlet boundary edge σ P E D ext to any prescribed cell K P T . Hence, the lower bound C 7 is found by minimizing π n K over K P T and n P t1, . . . , N u.
Lemma 3.5 is a crucial step in the proof of the existence of a solution p n " pp n K q KPT to the scheme (2.9)-(2.14).
Proposition 3.6. Given s n´1 " ps n´1 K q KPT P r0, 1s T , there exists a solution p n P R T to the scheme (2.9)-(2.14).
The proof relies on a standard topological degree argument and is omitted here. However, we make the homotopy explicit for readers' convenience. Let γ P r0, 1s be the homotopy parameter. We define the nondecreasing functions η pγq i : r0, 1s Ñ R`by setting η pγq i psq " p1´γq{µ`γη i psq for s P r0, 1s, and we seek a solution p pγq " pp pγq K q KPT to the problem where the fluxes F pγq Kσ are defined by with ϑ pγq " p pγq`ψ and using the upwind mobilities (3.17c) At the Dirichlet boundary edges, we still set p pγq σ " p D σ . For γ " 0, the system is linear and invertible, while for γ " 1, system (3.17) coincides with the original system (2.9)-(2.14). A priori estimates on p pγq that are uniform w.r.t. γ P r0, 1s (but not uniform w.r.t. T nor ∆t n ) can be derived on the basis of what was exposed previously, so that one can unfold Leray-Schauder's machinery [20,37] to prove the existence of (at least) one solution to the scheme.

Uniqueness of the discrete solution
To complete the proof of Theorem 2.2, it remains to show that the solution to the scheme is unique. This is the purpose of the following proposition.
By a similar argument, we can show that H n K pp n K Kr p n K , pp n Kσ Kr p n Kσ q σPE K q ě 0, @K P T , (3.21) where aKb " minpa, bq. Subtracting (3.21) from (3.20) and summing over K P T , we find where s n K " S K pp n K q, r s n K " S K pr p n K q and R n σ " η K ps n K Jr s n K qpϑ n K J r ϑ n K´ϑ n σ q`´η K ps n σ qpϑ n σ´ϑ n K J r ϑ n K qὴ K ps n K Kr s n K qpϑ n K K r ϑ n K´ϑ n σ q``η K ps n σ qpϑ n σ´ϑ n K K r ϑ n K q`, (3.23) with s n σ " S K pp n σ q. The top line of (3.23) expresses the upwinded flux of (3.20), while the bottom line of (3.23) is the opposite of the upwinded flux of (3.21). Note that, since p n σ " p D σ is prescribed at σ P E D ext , we have ϑ n σ " ϑ n σ J r ϑ n σ " ϑ n σ K r ϑ n σ . Upon inspection of the rearrangement R n σ " rη K ps n K Jr s n K q´η K ps n K Kr s n K qspϑ n K J r ϑ n K´ϑ n σ qὴ K ps n K Kr s n K qrpϑ n K J r ϑ n K´ϑ n σ q`´pϑ n K K r ϑ n K´ϑ n σ q`s η K ps n σ q " pϑ n σ´ϑ n K K r ϑ n K q`´pϑ n σ´ϑ n K J r ϑ n K q`‰, (3.24) it is trivial that R n σ ě 0. As a consequence, (3.22) implies that R n σ " 0 for all σ P E D ext and that s n K " r s n K for all K P T . At this stage, however, we cannot yet claim that p n K " r p n K , as the function S K is not invertible.
Taking into account s n K " r s n K , the residue (3.24) becomes R n σ " η K ps n K qrpϑ n K J r ϑ n K´ϑ n σ q`´pϑ n K K r ϑ n K´ϑ n σ q`s η K ps n σ q " pϑ n σ´ϑ n K K r ϑ n K q`´pϑ n σ´ϑ n K J r ϑ n K q`‰, (3.25) which can be lower-bounded by R n σ ě minpη K ps n K q, η K ps n σ qq|ϑ n K´r ϑ n K | (3.26) thanks to the algebraic identities a`´p´aq`" a and aJb´aKb " |a´b|. In view of the lowerbound on the discrete pressures of Lemma 3.5, we deduce from (1.6b) that s n K ą 0 and r s n K ą 0. The increasing behavior of η K implies, in turn, that η K ps n K q ą 0 and η K pr s n K q ą 0. Therefore, the conjunction of R n σ " 0 and (3.26) yields ϑ n K " r ϑ n K and hence p n K " r p n K for all cells K having a Dirichlet boundary edge, i.e., E K X E D ext ‰ H. It remains to check that p n K " r p n K , or equivalently ϑ n K " r ϑ n K for those cells K P T that are far away from the Dirichlet part of the boundary. Subtracting (3.19) from (3.18) and recalling that s n K " r s n K , we arrive at Consider a cell K P T where ϑ n K´r ϑ n K achieves its maximal value, i.e., ϑ n K´r ϑ n K ě ϑ n L´r ϑ n L , @L P T .

Convergence analysis
Once existence and uniqueness of the discrete solution have been settled, the next question to be addressed is the convergence of the discrete solution towards a weak solution of the continuous problem, as the mesh-size and the time-step are progressively refined. In accordance with the general philosophy expounded in [27], the proof is built on compactness arguments. We start by highlighting compactness properties in §4.1, before identifying the limit values as weak solutions in §4.2.

Compactness properties
Let us define G Em,∆tm : Q T Ñ R d and J Em,∆tm : for σ P E i,m , 1 ď n ď N m and, respectively, for σ P E m , 1 ď n ď N m . We remind that s Tm,∆tm " Spp Tm,∆tm , xq is the sequence of approximate saturation fields computed from that of approximate pressure fields p Tm,∆tm by (2.15b).
Proposition 4.1. There exists a measurable function p : Q T Ñ R such that Υppq´Υpp D q P L 2 pp0, T q; V q and Θ i ppq P L 2 pp0, T q; H 1 pΩ i qq, 1 ď i ď I, such that, up to a subsequence, Proof. We know from Corollary 3.3 that Θ i pp Tm,∆tm q and Υpp Tm,∆tm q are bounded w.r.t. m in L 2 pQ i,T q and L 2 pQ T q respectively, while G Em,∆tm and J Em,∆tm are respectively bounded in L 2 pQ i,T q d and L 2 pQ T q d . In particular, there exist p Θ i P L 2 pQ i,T q, p Υ P L 2 pQ T q, J P L 2 pQ i,T q d , and J P L 2 pQ T q d such that Establishing that p Θ i P L 2 pp0, T q; H 1 pΩ i qq and p Υ P L 2 pp0, T q; H 1 pΩqq with G " ∇ p Θ i and J " ∇ p Υ is now classical, see for instance [25,Lemma 2] or [16,Lemma 4.4].
The key points of this proof are the identification p Θ i " Θ i ppq and p Υ " Υppq for some measurable p, as well as the proofs of the almost everywhere convergence property (4.3a). The identification of the limit and the almost everywhere convergence can be handled simultaneously by using twice [2,Theorem 3.9], once for Θ i ppq and once for Υppq. More precisely, Lemma 3.4 provides a control on the time variations of the approximate saturation s Tm,∆tm , whereas Corollary 3.3 provides some compactness w.r.t. space on Θ i pp Tm,∆tm q and Υpp Tm,∆tm q. Using further that s Tm,∆tm " S i˝Θ´1 i pΘ i pp Tm,∆tm qq with S i˝Θ´1 i nondecreasing and continuous, then one infers from [2, Theorem 3.9] that s Tm,∆tm ÝÑ Let p " Θ´1 i p p Θ i q. Then, (4.3a) and (4.3b) hold. Proving (4.3a) and (4.3c) is similar, and the properties (4.3) can be assumed to hold for the same function p up to the extraction of yet another subsequence.
Finally, by applying the arguments developed in [6, §4.2], we show that Υppq and Υpp D q share the same trace on p0, T qˆΓ D , hence Υppq´Υpp D q P L 2 pp0, T q; V q.
Since η is bounded, Lebesgue's dominated convergence theorem ensures that the convergence holds in L q pQ T q for all q P r1,`8q. The reconstruction η Em,∆tm of the mobility is also uniformly bounded, so we have just to show that η Tm,∆tm´ηEm,∆tm L 1 pQ T q Ñ 0 as m Ñ`8. Letting ∆ Kσ " K X ∆ σ denote the half-diamond cell, we have η Tm,∆tm´ηEm,∆tm L 1 pQ T q ď Let us define r Em,∆tm pt, xq " |η n K´η n Kσ | " r n σ if pt, xq P pt n´1 m , t n m sˆ∆ σ , then r Em,∆tm is uniformly bounded by η 8 " 1{µ. Therefore, where h Tm is the size of T m as defined in (2.8). Besides, for i P t1, . . . , Iu, η i˝Si˝Θ´1 i is continuous, monotone and bounded, hence uniformly continuous. This provides the existence of a modulus of continuity i : R`Ñ R`with i p0q " 0 such that r n σ :" |η˝S˝Θ´1 i pΘ n K q´η˝S˝Θ´1 i pΘ n Kσ q| ď i p|Θ n K´Θ n Kσ |q (4.7) for σ P E i,m . Therefore, if the function q Ei,m,∆tm pt, xq " for σ P E i,m , 1 ď n ď N m , could be proven to converge to 0 almost everywhere in Q i,T , then it would also be the case for r Em,∆tm and R i,m as m Ñ`8, thanks to Lebesgue's dominated convergence theorem. Now, it follows from (3.12a) and from the elementary geometric relation Therefore, q Ei,m,∆tm Ñ 0 in L 2 pQ i,T q, thus also almost everywhere up to extraction of a subsequence. This provides the desired result.

Identification of the limit
So far, we have exhibited some "limit" value p for the approximate solution p Tm,∆tm in Proposition 4.1. Next, we show that the scheme is consistent with the continuous problem by showing that any limit value is a weak solution. Proof. Let ϕ P C 8 c ptΩ Y Γ N uˆr0, T qq and denote by ϕ n K " ϕpt n m , x K q, for all K P T m and all n P t0, . . . , N m u. We multiply (2.10) by ∆t n m ϕ n´1 K and sum over n P t1, . . . , N m u and K P T m to obtain A m`Bm " 0, m ě 1, (4.9) where we have set Thanks to the regularity of ϕ, the function δϕ Tm,∆tm converges uniformly to B t ϕ on Ωˆr0, T s. Moreover, by virtue of (4.3a) and the boundedness of s Tm,∆tn we can state that and, in view of the definition (2.13) of s 0 Tm and of the uniform convergence of ϕ 0 Tm towards ϕp0,¨q, ż From the above, we draw that Let us now turn our attention to the quantity B m of (4.10b), which can be split into B m " B Consider first the convective term B 2 m . It follows from the definition of the discrete gravitational potential ψ K "´ g¨x K , ψ σ "´ g¨x σ , K P T m , σ P E D ext,m and from the orthogonality of the mesh that Therefore, B 2 m can be transformed into where H Em,∆tm pt, xq " pd{d σ qpϕ n´1 Kσ´ϕ n´1 K qn Kσ if pt, xq P rt n´1 m , t n m qˆ∆ σ .
After [16,Lemma 4.4], H Em,∆tm converges weakly in L 2 pQ T q d towards ∇ϕ, while λ Em and η Em,∆tm converge strongly in L 4 pΩq and L 4 pQ T q towards λ and ηpSpp, xqq respectively (cf. Lemma 4.2). Thus, we can pass to the limit in (4.12) and The capillary diffusion term B 1 m appears to be the most difficult one to deal with. Taking inspiration from [13], we introduce the auxiliary quantity Kσ q.
Analogously to [22], we can define a piecewise-constant vector field H Em,∆tm such that and such that H Em,∆tm converges uniformly towards ∇ϕ on Q T . Under these circumstances, r B 1 where Ω i,m " Ť σPEi,m ∆ σ Ă Ω i . The strong convergence of ? η Em,∆tm in L 2 pQ i,T q towards a η i pS i ppqq directly follows from the boundedness of η i combined with (4.3a). Combining this with (4.3b) results in (4.14) Therefore, to finish the proof of Proposition 4.3, it only remains to check that B 1 m and r B 1 m share the same limit. To this end, we observe that, by the triangle inequality, we have Applying the Cauchy-Schwarz inequality and using Proposition 3.1, we find so R Γ,m Ñ 0 as m Ñ`8. Besides, we also apply the Cauchy-Schwarz inequality to R i,m in order to obtain |R i,m | 2 ď C 1 where we have set Define r η Em,∆tm pt, xq " Reproducing the proof of Lemma 4.2, we can show that Therefore, R i,m Ñ 0 as m Ñ`8. Putting things together in (4.15), we conclude that B 1 m and r B 1 m share the same limit, which completes the proof of Proposition 4.3.

Practical aspects of numerical resolution
We provide some details on the resolution strategy for the discrete problem (2.9)-(2.11c). It is based on a parametrization technique to automatically choose the most convenient variable during the Newton iterations ( §5.1) and on the addition of cells on the interfaces between different rock types ( §5.2) to improve the pressure continuity.

Switch of variable and parametrization technique
A natural choice to solve the nonlinear system (2.9)-(2.11c) is to select the pressure pp K q KPT as primary unknown and to solve it via an iterative method such as Newton's one. Nevertheless, the pressure variable is known to be an inefficient choice for s ! 1 because of the degeneracy of Richards' equation. For dry soils, this strategy is outperformed by schemes in which saturation is the primary variable. On the other hand, the knowledge of the saturation is not sufficient to describe the pressure curve in saturated regions where the pressure-saturation relation cannot be inverted. This motivated the design of schemes involving a switch of variable [21]- [32]. In this work, we adopt the technique proposed by Brenner and Cancès [5], in which a third generic variable τ is introduced to become the primary unknown of the system. Then the idea is to choose a parametrization of the graph tp, Sppqu, i.e., to construct two functions s : I Ñ rs rw , 1´s rn s and p : I Ñ R such that spτ q " Spppτ qq and s 1 pτ q`p 1 pτ q ą 0 for all τ P I Ă R. Such a parametrization is not unique, for instance one can take I " R, p " Id which amounts to solving the system always in pressure, but this is not recommended as explained before. Here, we set I " ps rw ,`8q and where S 1 ppś q denotes the limit as p tends to p s " Sps s q from below of S 1 ppq. Since the switch point s s is taken as the inflexion point of S, both s and p are C 1 and concave, and even C 2 if S is given by the Van Genuchten model. Moreover, for all p P R, there exists a unique τ P ps rw ,`8q such that pp, Sppqq " pppτ q, spτ qq. The resulting system F n pτ n q " 0 made up of N T " CardpT q nonlinear equations admits a unique solution τ n , since it is fully equivalent to (2.9)-(2.11c). More details about the practical resolution of this nonlinear system via the Newton method can be found in [4].

Pressure continuity at rock type interfaces
Physically, the pressure should remain continuous on both sides of an interface between two different rock types. But this continuity is here not imposed at the discrete level. The two-point flux approximation based on the cell unknowns is strongly dependent on the mesh resolution and can induce a large error close to the rock type interface. We here propose a very simple method to improve this continuity condition in pressure. It consists in adding two thin cells of resolution δ around the rock-type interface with δ ! ∆x as shown in Figure 1. The idea is here to add two Figure 1: Mesh refinement on both sides of an interface face for a 2D case.
unknowns in the neighborhood of the interface to have a more precise approximation of the pressure gradient on each side of the faces where changes of rock types occur. In this way, we avoid the introduction of face unknowns in our solver which remains unchanged. Alternatives that strictly impose the pressure continuity will be studied and compared with this approach in a forthcoming work.

Numerical results
In this section, we present the results obtained for different test cases. For all these cases, we consider a two-dimensional layered domain Ω " r0m, 5msˆr´3m, 0ms made up of two rock types denoted by RT0 and RT1 respectively, RT0 being less permeable than RT1. Using these two lithologies, the domain Ω is partitioned into three connected subdomains: Ω 1 " r1m, 4msˆr´1m, 0ms, Ω 2 " r0m, 5msˆr´3m,´2ms and Ω 3 " Ω z pΩ 1 Y Ω 2 q, as depicted in Figure 2. The Brooks-Corey [10] and van Genuchten-Mualem [42] petro-physical models are used to model the flow characteristics of both rock types. In these models, the water saturation and the capillary pressure are linked pointwise by the relation s " Sppq where S : R Ñ r0, 1s is nondecreasing and satisfies Sppq " 1´s rn if p ě p b and Sppq Ñ s rw as p Ñ´8, s rw being the residual wetting saturation, s rn the residual non-wetting saturation and p b the entry pressure. More precisely, we have, • for the Brooks-Corey model, • for the van Genuchten-Mualem model, where ηp¨q " k r p¨q{µ, µ " 10´3 Pa¨s being water viscosity, is the relative permeability. The parameters used for both rock types are given in Table 1 for the Brooks-Corey model and in Table 2 for the van Genuchten-Mualem model. With these choices of parameters, water is more likely to be in RT1 than in RT0, in the sense tha, at a given pressure, the water saturation is higher in RT1 than in RT0, as it can be seen on the plots of the capillary-pressure functions depicted in Figures 3-4 for these two petro-physical models. Figures 3-4 also show the relative permeability functions. Note the non-Lipschitz character of the relative permeability in the van Genuchten-Mualem framework. For the numerical tests, in order to avoid infinite values for the derivative of k r psq when s Ñ 1´s rn , we approximate it for s P rs lim , 1´s rn s using a second degree polynomial r k r psq. Such a polynomial satisfies the following constraints: k r ps lim q " r k r ps lim q and r k r p1´s rn q " 1. The value s lim corresponds to s eff " 0.998.

Configurations of the test cases
For both petro-physical models, we consider two configurations further referred as filling and drainage cases, which are described in the following.

Filling case
The filling test case has already been considered in [15,32,36,39]. Starting from an initially dry domain Ω, whose layers' composition is reported in Figure 5, water flows from a part of the top boundary during one day. A no-flow boundary condition is applied elsewhere. More precisely, the initial capillary pressure is set to´47.088¨10 5 Pa and the water flux rate to 0.5m/day through Γ N " tpx, yq | x P r1m, 4ms, y " 0mu. For this simulation an homogeneous time-step ∆t " 1000s is prescribed for the test using the Brooks-Corey model and ∆t " 500s for the one using the van Genuchten-Mualem model. The test case follows the following dynamics. Water starts invading the void porous space in Ω 1 . When it reaches the interface with Ω 3 , capillarity involves a suction force on water from Ω 1 to Ω 3 . Since clay (RT1) has low permeability, water encounters difficulties to progress within Ω 3 .
This yields a front moving downward in Ω 1 which is stiffer for the Brooks-Corey model than for the van Genuchten-Mualem one. In both cases, the simulation is stopped before water reaches the bottom part corresponding to Ω 2 . In Figure 6 we can observe the evolution of the saturation profile during the simulation performed on a 50ˆ30 cells mesh with Brooks and Corey model, whereas the evolution corresponding to van Genuchten-Mualem nonlinearities is depicted in Figure 7.

Drainage case
This test case is designed as a two-dimensional extension of a one-dimensional test case proposed by [38] and addressed in [15,39]. We simulate a vertical drainage starting from initially and boundary saturated conditions during 105¨10 4 s. At the initial time, the pressure varies with depth with p 0 pzq "´ρgz. A Dirichlet boundary condition p D " 0 Pa is imposed on the bottom of the domain, more precisely on Γ D " tpx, yq | x P r0m, 5ms, y "´3mu. The layers' composition of Ω is reported in Figure 8. For this simulation an homogeneous time-step ∆t " 2000s is used for the test with the Brooks-Corey model and ∆t " 800s for the one with the van Genuchten-Mualem model. At the top interface between Ω 1 and Ω 3 , capillarity acts in opposition to gravity and to the evolution of the system into a dryer configuration. The interface between Ω 2 and Ω 3 acts in the reverse way: suction accelerates the gravity driven drainage of RT0.
In Figure 9 we can observe the evolution of the saturation profile during the simulation performed on a 50ˆ30 cells mesh with Brooks and Corey model, whereas the evolution corresponding to van Genuchten-Mualem nonlinearities is depicted in Figure 10.

Comparisons of the numerical treatments of the interfaces
For each petro-physical model and configuration, a numerical convergence analysis is carried out for the schemes with (method A) or without (method B) thin cells, whose thickness is fixed to δ " 10´6m, at rock type interfaces. Five structured meshes with the following resolutions are considered for this analysis: 50ˆ30, 100ˆ60, 200ˆ120, 400ˆ240, 800ˆ480. The evolution of the error is measured using the L 2 pr0, T s, Ωq-norm of the relative difference between the saturations obtained on a given mesh and a reference solution obtained with Method A and the mesh 800ˆ480. The number of Newton iterations obtained with both methods is also compared.

Brooks-Corey model: drainage case
For the drainage case with the Brooks-Corey model, the convergence error is given in Figure 11. First we notice that, for all meshes, the error is smaller with method A than with method B and that we have a linear rate of convergence with the first one whereas this rate is smaller with the latter one. The total, average and maximal number of Newton iterations are also given in Table  3. Method A appears to be slightly more expensive.  Figure 11: L 2 pQ T q relative error in saturation for the drainage case using Brooks and Corey model.
Let us now evaluate the saturation absolute error between results obtained with Method A and Method B. In Figure 12 we plot the absolute-error distribution over the domain at three different 7 total 7 avg 7 max Method A 2038  3  29  Method B 1927  3  29   Table 3: Newton's iterations for the mesh 200ˆ120 for the drainage case using Brooks and Corey model.
times: when the cells line in Ω 1 above the interface between Ω 1 and Ω 3 starts drying, when the cells line in Ω 2 below the interface between Ω 3 and Ω 2 starts drying and at final time.

Brooks-Corey model: filling case
For the filling case with the Brooks-Corey model, the convergence error is given in Figure 13. As for the previous case, Method A enables to recover a linear convergence rate. Except for the first two meshes where the error obtained with Method A is slightly larger, for all other meshes, this error is smaller than the one obtained with method B. The total, average and maximal number of Newton iterations are given in Table 4. The algorithm behaves here in the same way as before.  Figure 13: L 2 pQ T q relative error in saturation for the filling case using Brooks and Corey model.
Let us now evaluate the saturation absolute error between results obtained with Method A and Method B. In Figure 14 we plot the absolute-error distribution over the domain at three different times: when water crosses the interface between Ω 1 and Ω 3 , when cells around this interface are almost saturated and at final time.

Van Genuchten-Mualem model: filling case
For the filling case with the Van Genuchten model, the convergence error is given in Figure 15. Both methods exhibit a linear rate of convergence. On the other hand, the error is slightly larger with method A than with method B. The total, average and maximal number of Newton iterations are given in Table 5.  Figure 15: L 2 pQ T q relative error in saturation for the filling case using the Van Genuchten model. Figure 16 shows the localization of the difference for the numerical solutions provided by the two methods A and B. Unsurprisingly, the difference is located in the neighborhood of the interfaces. Moreover, as suggested by Figures 13 and 15, the influence of the introduction of additional interface unknowns (method A) has a lower impact for van Genuchten-Mualem nonlinearities than for Brook-Corey nonlinearities.  Table 5: Newton's iterations for the mesh 200ˆ120 for the filling case using the Van Genuchten model. Figure 16: Saturation absolute error between Method A and Method B for the filling case with the Van Genuchten model at t P t27.5¨10 3 s, 45¨10 3 s, 86¨10 3 su.

Van Genuchten-Mualem model: drainage case
For the drainage case with the Van Genuchten model, the convergence error is given in Figure  17. Both methods exhibit a linear rate of convergence. Moreover, the error is slightly larger with method B than with method A. The total, average and maximal number of Newton iterations are given in Table 6. Let us now evaluate the saturation absolute error between results obtained with Method A and Method B. In Figure 18 we plot the absolute-error distribution over the domain at three different times: when the cells line in Ω 1 above the interface between Ω 1 and Ω 3 starts drying, when the cells line in Ω 2 below the interface between Ω 3 and Ω 2 starts drying and at final time.  Table 6: Newton's iterations for the mesh 200ˆ120 for the filling case using the Van Genuchten model.

Influence of the parameter δ
Let us now analyze how the thickness of the thin cells employed in Method A affects the accuracy of the solution obtained with this method. We consider the filling and drainage cases along with the Brooks and Corey model and evaluate the relative L 2 pQ T q error between the solution obtained on the 200ˆ120 cells mesh using δ P t10´2m, 10´4m, 10´6mu with respect to the reference solution obtained on the 800ˆ480 cells mesh using δ ref " 10´6m. As shown in Figure 19, the value of δ does not have a significant influence on the overall error as soon as δ is small enough. We also observe a moderate influence on the robustness of the non-linear solver for the values considered here.  Figure 19: L 2 pQ T q relative error in saturation as a function of the thickness δ of the thin cells with Method A using the 200ˆ120 cells mesh.

Conclusions and perspectives
This article aimed at proving that standard upstream mobility finite volume schemes for variable saturated porous media flows still converge in highly heterogeneous contexts without any specific treatment of the rock type discontinuities. The scheme is indeed shown to satisfy some energy stability which provides enough a priori estimates to carry out its numerical analysis. First, the existence of a unique solution to the nonlinear system stemming from the scheme is established thanks to a topological degree argument and from the monotonicity of the scheme. Besides, a rigorous mathematical convergence proof is conducted, based on compactness arguments. No error estimate can then be deduced from our analysis.
Because of the choice of a backward Euler in time discretization and from the upwind choice of the mobilities, a first order in time and space accuracy is expected in the case of homogeneous computational domains. We show in numerical experiments that without any particular treatment of the interfaces at rock discontinuities, this first order accuracy can be lost, especially in the case of Brooks-Corey nonlinearities. This motivates the introduction of a specific treatment of the interfaces. The approach we propose here is based on the introduction of additional unknowns located in fictitious small additional cells on both sides of each interface. Even though the rigorous convergence proof of this approach is not provided here in the multidimensional setting -such a proof can for instance be done by writing the scheme with the specific treatment of the interface (method A) as a perturbation of the scheme without any particular treatment of the interface (method B)-, the numerical experiments show that it allows to recover the first order accuracy without having major impacts on the implementation and on the behavior of the numerical solver.
For future researches, we suggest to test the so-called method A on a two-phase flow test and to compare it to the approaches presented in [7]. Moreover, in a forthcoming work, we propose two other methods to really impose the pressure continuity condition at interfaces. A comparison between all methods will be shown.