Numerical Study for Two-Phase Flow with Gravity in Homogeneous and Piecewise-Homogeneous Porous Media

In this work we study two-phase flow with gravity either in 1-rock homogeneous media or 2-rocks composed media. These phenomena can be modeled by a non-linear scalar conservation law with continuous flux function or discontinuous flux function, respectively. Our study is essentially from a numerical point of view, we apply the new Lagrangian-Eulerian (LEH) finite difference method developed by Abreu and Pérez [1, 2, 3, 4, 5, 25] and the classical Lax-Friedrichs method to obtain numerical entropic solutions. Comparisons between numerical and analytical solutions show the efficiency of the methods even for discontinuous flux function. Our main contribution, is the comparison and error analysis between the new LEH and the classical Lax-Friedrichs (LF) methods, in order to show the good performance of the LEH scheme for models with discontinuous flux functions.


INTRODUCTION
Many problems in engineering, physics and other areas of sciences lead us to the study of conservation laws. For instance, they model many physical phenomena that appear in aerodynamics, fluid mechanics, traffic flow, groundwater flow, multi-phase flow in porous media and others [7,8,9,16,19,20,26,29,30,32]. In general a scalar conservation law in one dimension takes the form: In the case that the medium is composed of two types of rocks, the flux function has a spatial discontinuity and the Riemann solutions contain in general an additional stationary shock [6,11]. Although, for the model studied in [11], the analytical entropic solutions can be obtained through an extension of Oleinik construction with appropriate connections between the fluxes [6,11], there is no an entropy condition that works in general for all discontinuous flux functions. Thus, the usage of appropriate numerical methods to obtain the Riemann solutions of hyperbolic conservation laws with discontinuous flux function gain importance.
Although there are several good numerical methods to deal with discontinuous flux function in conservation laws theory (see for instance [13,14,21,22] among others), up to now there are only two finite difference schemes whose convergence to entropic solutions have been rigorously proved for general space-discontinuous flux functions [1,13,25]: the classical Lax-Friedrichs scheme [13] and the new Lagrangian-Eulerian (LEH) scheme developed by Abreu and Pérez [1,2,3,4,5,25] which is a shock-capturing and high-resolution method for first order hyperbolic problems (maybe with forcing terms, because LEH method has also the "well-balance property when adapted for balance laws). The new scheme provides high-order resolution in smooth regions and prevents the creation of spurious oscillations near discontinuities. The new Lagrangian-Eulerian scheme satisfies Harten's requirements [26] of being a second-order accurate TVD scheme. All these features turn the LEH method a very interesting option to deal with conservation or balance laws with discontinuous flux functions. Our main contribution in this paper, is the comparison and error analysis between the new LEH and the classical Lax-Friedrichs (LF) methods, in order to show the good performance of the LEH scheme for models with discontinuous flux functions.
Motivated by this, in this work we solve some Riemann problems by both the Lagrangian-Eulerian and the Lax-Friedrichs methods for two-phase flow with gravity in two cases; first we test our codes for the homogeneous rock case and then we obtain the numerical solutions for a piecewise-homogeneous rock (two-rocks composed media). In the homogeneous test-case we compare the numerical solutions with the analytical solutions obtained via the classic Oleinik construction. In the case of two-rocks composed media, we use the analytical entropic solutions appearing in [11] and briefly explained in Section 3 of this paper, to calculate the errors of the methods.

IMMISCIBLE TWO-PHASE FLOW WITH GRAVITY
We study the one-dimensional flow in porous media of two immiscible fluids: water and oil, under the effect of gravity force. Initially there is an impermeable interface, in z = 0, separating the mixture of the fluids and we study the flow near the interface, once the interface disappears. This is equivalent to solve a class of Riemann problems for a scalar conservation law.
We assume also that the flow is isothermal, there are no sources or sinks, the porosity is constant, compressibility and capillary pressure effects are neglected. Equation (2.1) expresses the conservation of mass: Here, φ denotes porosity, s i , v i , respectively denote saturation and seepage velocity of the phase i. The seepage velocity is given by Darcy's law (2.2), in which K is the absolute permeability of the rock, k i , µ i , p i and ρ i are the fluid relative permeability, viscosity, pressure and density for each phase i. Also in (2.2) g denotes the gravitational constant. The viscosities µ i are assumed constants, while the permeabilities k i are assumed to be functions of the saturations.
Under the assumption that the porous medium is totally saturated (s w + s o = 1) the equation (2.1)-(2.2) reduce to the following scalar conservation law.
with flux function (2.4) here we denote µ = µ w /µ o and ρ = ρ o /ρ w , s w is the saturation of the phase water. The dimensionless parameter v is related to the pressure gradients.
For two-rock composed media, the porous medium is piecewise-homogeneous, so we define the permeability of the rock by the following equation which leads to a discontinuous flux function of the form The initial Riemann data for equations (2.3)-(2.5) are denoted by

ANALYTICAL SOLUTIONS FOR SCALAR CONSERVATION LAWS
For non-linear scalar conservation laws, discontinuous or non-smooth solutions commonly appear even for smooth initial data. Thus a weak formulation of the conservation law and the corresponding "weak" concept of solutions are needed to deal with this kind of equations. Specially for Riemann problems the weak solutions are combinations of shocks and rarefaction waves. The shocks (even for smooth initial data) happen when characteristics cross, at the time where the characteristics first cross, the function s(z,t) has an infinite slope -the wave "breaks" and a shock forms [17]. The shock waves joining states s l and s r travel with constant speeds σ which satisfy the Rankine-Hugoniot condition: Unfortunately weak solutions are not unique and additional conditions are needed in order to select the physically correct one (entropic solution). For instance we consider the Oleinik entropy condition which apply to general non convex flux f : the solution s(z,t) is the entropy solution if all discontinuities have the property that [17,23] for all s between s l and s r and σ is gives by equation (3.1). The above exposed, is the basement of the Oleinik's geometric construction which allows to obtain the entropic analytic solutions for any Riemann problem in the one-dimensional scalar case for continuous flux functions (the case of two-phase flow in homogeneous porous medium).
Thus the entropic solution can be constructed by using the convex or concave hulls of the flux function, depending on the relation between the initial Riemann left and right states, see the for a right characteristic shock wave or the zero of the equation for a left characteristic shock wave. Figure 1: Example of flux function f and entropic shock waves as limits of travelling waves. Source: [30].
In [11] Kaasschieter considered a generalization for discontinuous flux functions in space variable z Here f l , f r : [0, 1] → R are assumed to be twice differentiable such that f l (0) = f r (0) and f l (1) = f r (1). The states are the left and right states of the stationary shock always appearing in the solutions. They are called fluxes connections in literature [6,11]. The last expression in (3.5) represents the Rankine-Huguniot condition for stationary shocks.
Kaasschieter obtained entropy conditions for a model of two-phase flow with gravity in two-rocks of the form (3.5) with flux functions in (2.5)-(2.7). We will use the same entropy conditions to construct our analytical solutions in Section 5.2.
To obtain the entropy conditions, Kaasschieter [11] regularised the problem by adding the capillary diffusion term, he considered similar solutions s(z,t) = s(ξ ), with ξ = z/t for the regularized problem, he assumed some hypothesis like S L = S R and p c,L (S L ) = p c,R (S R ) (distinct capillary pressures at the connection states), between others, and after some calculations he proved that: (i) All non-stationary waves must satisfy the classical Oleinik entropy condition for one of the flux functions f l or f r . Of course, waves with negative speeds are obtained through the convex or concave hulls of f l , while waves with positive speeds are obtained through the convex or concave hulls of f r .
(ii) There is a stationary shock separating the non-stationary waves groups. The stationary shock connecting states S L and S R can be obtained geometrically by using a horizontal jump line connecting the fluxes.
(iii) The appropriate (unique) selection of these connection states, requires an entropy condition. For this problem the states S L and S R must satisfy the entropy inequality where the notation f l (S L ) 0 is used if f l (S L ) ≥ 0 and lim s↑S L sign( f l (s)) = 1 or lim s↓S L sign( f l (s)) = 1. Analogously the notation f r (S R ) 0 means that f r (S R ) ≤ 0 and lim s↑S R sign( f r (s)) = −1 or lim s↓S R sign( f r (s)) = −1.
Thus, in order to obtain the entropic solutions for our model with discontinuous fluxes, the key fact is the appropriate selection of the connection states S L and S R of the stationary shock by using a horizontal jump connecting the fluxes while satisfying the entropy inequality in (3.7), the other waves can be obtained through the Oleinik construction for f l and f r , for illustration see

FINITE DIFFERENCE METHODS
Finite difference methods are numerical methods for solving differential equations, in particular very useful for conservation laws. This kind of method consists on replacing the derivatives in the differential equations by finite difference approximations obtained through Taylor expansions of the functions.
In order to preserve the conservation property of the solutions, finite differences schemes in conservative forms [15] are used: The function G is called the numerical flux.

Lax-Friedrichs Scheme
The classical Lax-Friedrichs (LF) scheme for conservation laws with continuous flux functions in conservative form (4.1) has the following numerical flux, see [18]: For this method to be convergent, it must satisfy the CFL condition [18] max For conservation laws with discontinuous flux functions the numerical flux of the method can be adapted to

Lagrangian-Eulerian Scheme
The Lagrangian-Eulerian (LEH) method developed by Abreu and Pérez [1, 2, 3, 4, 5, 25] is a scheme based on a Eulerian central scheme finite volume formulation, but in a space-time Lagrangian-Eulerian framework. These scheme is derived from the divergence forms of the equations [25], and it is an order high-resolution scheme for numerically solving nonlinear conservation law problems. This scheme provides high-order resolution in smooth regions and prevents the creation of spurious oscillations near discontinuities. The convergence of the numerical solutions to the entropic solutions for models with discontinuous flux functions was proved in [25], [1].
This method for hyperbolic conservation laws with continuous flux functions is given by the following numerical formulation, see [1,25] for the details: and the CFL condition for this scheme is given by The scheme (4.3) can be written in the conservative form (4.1), where the numerical flux is, see [1,25] G(U n j ,U n j+1 ) = For conservation laws with discontinuous flux functions the method use the following numerical flux with f in the form (4.2) and CFL (see [1,25] for details):

Two-Phase Flow in One-Rock Media
First we consider a test-case of homogeneous media, i.e., constant absolute permeability of the rock K. We proceed to compare the numerical solutions in the following two cases: (i) Purely gravitational flow (i.e., v = 0); (ii) Mixed flow with dominant gravity (v small and non-zero, for instance v = 0.01).
For both cases we set the following parameters µ = 0.25 and ρ = 0.8.
In order to illustrate the solution for the pure gravitational problem (case i), we consider initial Riemann data s l = 1 and s r = 0, see Figure 2 . The solution consists in a shock wave with negative speed joining s l = 1 with s 1 = 0.4732 traveling upwards in the reservoir. There is also a rarefaction wave with positive and negative speeds joining the states s 1 = 0.4732 with s 2 = 0.274, preceding a last shock wave with positive speed joining s 2 = 0.274 with s r = 0. This last shock moves to the right (bottom in the reservoir).
A numerical refinement study for this case is shown in Figures 2b-2d to illustrate that accuracy increases when the mesh is refined. In the Table 1, we show the errors between the approximate solution obtained by Lagrangian-Eulerian and Lax-Friedrichs methods when compared to the analytical Oleinik solution.
The variation of the error with respect to the time is shown in Figure 3, notice that as time increases the error tends to increase, and this is smaller in refined meshes.
In Figures 4 and 5 the solutions for other Riemann data in the case (i) are illustrated.
In the case (ii), for mixed flow with dominant gravity, the Riemann solutions for initial data s l = 1 and s r = 0, are shown in the Figure 6 and in Figure 7 for initial data s l = 0 and s r = 0.8.

Two-Phase Flow in Two-Rocks Media
For two-rocks composed media the conservation law has a discontinuous flux function (2.5)-(2.7). We consider some distinct combination of values for µ, ρ, K l , K r and we only study the pure gravitational flow, i.e. v = 0.
In the case 1, three simulations were performed for distinct values of µ and initial Riemann data, see the    two groups of waves, the first group has negative speeds and travels upward in the reservoir, it is composed by a shock wave joining the states s l = 0 and s 1 = 0.298 followed by a rarefaction wave up to state S L = 0.4095. The second group consists in a solely shock wave with positive speed travelling downwards in the reservoir connecting the states S R = 0.6659 and s l = 1. There is a stationary shock joining the states S L = 0.4095 with S R = 0.6659 separating the two groups of waves.
The errors of numerical solutions are shown in the Table 2, notice that Lagrangian-Eulerian scheme had better performance (more accurate) than the Lax-Friedrichs method and also that the error decreases when the mesh is refined, see the Figures 8b-8d.
The solution in Figure 9 has a structure similar to the solution in Figure 8. The last simulation corresponding to the case 1, see Figure 10, shows a solution without the negative-speed group of waves and only a stationary shock followed by a positive speed shock.
For the cases 2, 3 and 4 the simulations shown in the Figures 11, 12 and 13 respectively, illustrate again the well performance of the methods in particular the high accuracy of the Lagrangian-Eulerian scheme when compared with Lax-Friedrichs method. For all the cases the analytical solutions were obtained using the entropic condition (3.7) developed in [11] to choose the appropriate connection between the fluxes (stationary shock).       The Lax-Friedrichs and the Lagrangian-Eulerian methods were used to obtain the approximate Riemann solutions for two-phase flow with gravity in both homogeneous and piecewisehomogeneous porous media (two-rocks media). For two-rocks media the flux function has a discontinuity in space variable z so find analytical solutions is not simple. We used the entropy criterion (3.7) developed in [11] to compare with our numerical results. Both methods LF and LEH present a diffusive behavior, however the Lagrangian-Eulerian scheme showed to be less diffusive and more accurate as reflected in the table and figures shown in Section 5. Notice also that the error decreases when the mesh is refined.