CHARACTERIZATION OF THE SPEED OF A TWO-PHASE INTERFACE IN A POROUS MEDIUM

A typical situation of oil reservoir simulation is considered in a porous medium where the resident oil is displaced by water injection. An explicit expression of the speed of the oilwater interface is given in a pseudo-2D case via the resolution of an auxiliary Riemann problem. The explicit 2D solution is then corroborated with numerical simulations by solving the transport equation with a generalized scheme of Harten type.


Introduction
Let Ω be a rectangular domain in R 2 and let ∂Ω be its boundary.In this paper, we will be interested in the resolution of a system of type (P) with subject to the initial condition and the following boundary condition for the variable u: u(t,x) = 1 on Γ in (1.4) together with the conditions for p: Equations (1.1)-(1.3)represent a simplified model for an incompressible flow in some porous medium, analogous to the Buckley-Leverett model used in oil reservoir simulations [9].In this context, the system represents a water-oil displacement model, that is, oil trapped in Ω is pushed towards Γ out by injection of water at Γ in , see Figure 1.1.In our case, the variable u represents the saturation of the water phase, the vector V is the total flux density of the two fluids, and p is the pressure.See next paragraph for the derivation of the model.
The functions f and m are nonlinear functions of u.A typical representation of the fractional flow function f , which characterizes the relative permeability of the fluids, is given for some α > 0 by f (u) = λu α λu α + (1 − u) α  (1.6) (see Figure 1.2).In the above formula, λ is the viscosity ratio of the two fluids (λ = µ o /µ w ).Note that such an f , which is obtained by considering a relative permeability of the polynomial form u α , is an increasing function but which may change concavity.As concerning the function m, it will be supposed to be positive, bounded, that is, there exist m 0 and m 1 such that for u ∈ [0,1].More precisely, m, which represents the total mobility of the two fluids, is given by Multiphase flows in porous media are inherently nonlinear and often numerical simulations happen to be the only practical way to understand their qualitative behaviour.Even for simplified systems analogous to (1.1)-(1.2),usually named the Buckley-Leverett model, there is no general result concerning the existence and uniqueness of a global solution of such a system, see, for example, the book of Gagneux et al. [14].This hyperbolicelliptic-type model results by neglecting gravity and capillary forces in the generalized Darcy's law which is usually used for the expression of the velocity in the case of two immiscible fluids sharing a porous medium (cf.[2,5,10,14]).This leads to the expression (1.2) in our case.The consideration of the capillary forces would result in a parabolicelliptic version studied by various authors, see, for example, [13,14].As for the hyperbolic case, results which exist concern essentially treatment of particular boundary conditions.For example, in [14], the mobility m is taken to be a constant which comes to solving −∆p = 0 for the elliptic part.This together with convenient boundary conditions of Neumann type on Γ in and Robin-type condition on Γ out ensures a unique solution in C([0,T];L 1 (Ω)) for all T > 0. In their paper, Schroll and Tveito [27] are interested in the local existence and stability of solutions of such a system.With Neumann conditions on the whole boundary and appropriate regularity hypotheses, they indeed obtain a classical but local in time solution.
In this paper, we propose to study the well-posedness of the system (1.1)-(1.5)and to construct a weak entropy solution in the pseudo-2D case (i.e., when Γ in and Γ out are the lateral sides of the rectangular domain) via the resolution of an auxiliary 1D problem of Riemann type.
We start by describing the simplications called for in the physical problem which led to our model by departing from the theory of mixtures, see, for example, [2,25].Then in Section 3 we define an entropy solution of (1.1)-(1.5),and we verify that the problem is well-posed in the sense of Bardos-LeRoux-Nédélec (BLN) [6].In Section 4, we solve the corresponding one-dimensional problem which in turn allows us to give an entropy solution in the pseudo-2D case.In the last section, some numerical results are presented, to corroborate the theoretical ones, by using a scheme of Harten type.644 A two-phase interface in a porous medium

The theory of mixtures.
As we said in the introduction, the system (1.1)-(1.5)represents a simplified water-oil displacement model used in oil reservoir simulations for enhanced oil recovery.The latter consists in injecting water at a well (Γ in ) in order to push the oil trapped in the porous rock (Ω) towards the production well (Γ out ).
A complete theoretical model must therefore be able to take into account the various physical characteristics and any eventual interactions of the three phases in presence: the two fluids oil (o) and water (w) and the solid rock (r).The so-called theory of mixtures provides such a framework as it presents a systematic methodology to derive equations of problems such as flow through a porous medium.For an overview of models concerning multiphase flows, we refer to the book on Mechanics of Mixture by Rajagopal and Tao [25], and among other papers those of Atkin and Craine [4], Bowen [8], Truesdell [28], and Allen [2].
Let p be the index which denotes the phase (in our case p = o, w, or r) and let φ p , ρ p , and v p represent, respectively, the volume fraction, the mass density, and the velocity of the phase p. Then the porosity of the medium is given by ( The saturation S p of the fluid phase p is then defined by subject of course to the constraint 3) The governing equations of an isothermal flow will be obtained by writing the mass and momentum balance equations for the three phases.As in our case there is supposed to be no interphase mass transfer and, in particular, the rock is chemically inert, the mass conservation for any phase p will read The velocity field of the rock can be reasonably taken as zero (up to fixing a coordinate system in which v r = 0).In particular, we neglect any deformation which could have been caused by the water injection.As concerning the velocity fields of the fluids, they are obtained via the momentum balance equations which for Newtonian fluids are given by where P p is the mechanical pressure in the fluid p, g the gravitational acceleration, z some depth, and m p the momentum exchange from other phases.This latter quantity is A. Abbassi and G. Namah 645 commonly assumed to take the form where Λ p is the mobility tensor of the phase p.If we further assume that the inertial effects in the fluids are negligible (i.e., Dv p /Dt = 0), (2.5) yields the so-called Darcy's law Constitutive laws for mobility are largely phenomenological, the most common versions having the form where µ p is the dynamic viscosity of fluid phase p, k is the permeability tensor, and the relative permeability K r p is a coefficient describing the effects of the other fluid in obstructing the flow of fluid p.For a two-fluid system with no interphase mass transfer as in our case, K r p is typically a function of the saturation S p .

The Buckley-Leverett model.
The Buckley-Leverett model [9] is a simplified model of two phase flows in porous media but which is of particular relevance in the petroleum industry for enhanced oil recovery.This model is based on some basic assumptions, namely, the incompressibility of the solid and the fluids and the fact that the effects of capillary pressure gradients on the flow fields are negligible as compared to the pressure gradients applied through pumping.
The incompressibility assumption implies that φ is constant in time and that the fluids densities are constant in space and time.One can then simplify the continuity equation (2.4) by ρ p and it becomes If we introduce the flux function V p defined by and use (2.2), we get By summing the equations for the two fluids and using the fact that S o + S w = 1, we then have the system 646 A two-phase interface in a porous medium with V T = V o + V w , the total flux given by, due to (2.7), In practice, the interface tension, which occurs due to the difference of pressures of the two fluids, leads to capillary effects.In our two-phase system, there will be a single capillary pressure P c given by P c = P o − P w . (2.14) Now we invoke the assumption of the Buckley-Leverett model corresponding to the fact that capillarity has negligible effects on the flow field-wide, that is, gradP c 0 so that (2.15) Finally assuming that gravity effects are absent as in the 1D case, the expression of V T collapses to Now define the fractional flow of water, f = f (S w ), by (2.17) By noting that V w = V T f , we thus end up with the hyperbolic-type equation with V T given by (2.16).If we now define the total mobility by m, that is, then one sees that it suffices to replace in the above equations the porosity φ by ε, the unknown S w by u, and the total flux function V T by V to get (1.1)-(1.2).The fractional flow f , which is clearly nonlinear through its dependence on the unknown S w , has typically an "S-shaped" profile such as the one depicted in (1.6).As concerning the total mobility function m, due to (2.8), it is given by (2.21) The expression (1.8) considered in this paper is thus obtained by considering an identity permeability tensor k (which comes to assuming that the rock is homogeneous without any particular isotropy) and a polynomial form of the relative permeability, K rw = S α w .
Let us finally remark that not neglecting the capillary forces would lead to completely different models on a mathematical basis.Usually one ends up with quasilinear parabolic degenerate equations instead of hyperbolic ones-see, for example, Chavent and Jaffre [10] where a fictitious global pressure is introduced as a new unknown to weaken the coupling between S and P.

Definition and well-posedness
Let us first say a few words on the boundary conditions.Notice that the boundary condition for u is given only at Γ in .In fact, as we will see further, this is sufficient for the well-posedness of the problem in the sense of BLN [6].As concerning the boundary conditions on p, we have p = 1 at Γ out corresponding to an atmospheric pressure, which is physically meaningful.But without loss of generality for the analysis, we will consider from now on p = 0 at Γ out .The condition on the rest of the boundary is n being the outward unit normal.As for g, it will be taken as follows: with some positive constant q.In this context, q represents the injection speed of the water at Γ in , the rest of the boundary being impermeable except of course for Γ out .Let us note that such a boundary condition is indeed not a regular one (g is not continuous) and the latter could have been regularized but we will not do so as regularity results are not our aim here.Instead we prefer a context which physically makes sense and deals with the weak formulation below.Define Q = Ω×]0,T[ for some time T > 0 and let BV be the space of functions which have bounded variations.Consider then the following problem for a given u in Of course, such a p will depend on t via the function u(•,t).Note that as u ∈ BV (Q), u(•,t) can be defined for all t > 0 in the sense of trace.
and, as m 0 > 0, Poincaré's inequality which holds on ᐂ leads to We then conclude by Lax-Milgram for the existence and uniqueness of p.
Let then which leads to the fact that p − is almost everywhere constant on Ω.And, as p − |Γ out = 0, we finally have p ≥ 0 a.e. on Ω.
Note that, as p can only be decreasing across Γ out , the above proposition ensures that it suffices to impose a boundary condition for u only on Γ in for the problem to be well-posed in the sense of BLN [6].Indeed the correct way to prescribe the boundary conditions for the hyperbolic part in our case happens to be min This last equality will always hold if we impose γ(u) = 1 on the part of the boundary where V • n < 0, that is, on Γ in .We now proceed by giving the definition of a solution of (1.1)-(1.5).For this sake and without loss of generality, we will shift to homogeneous Dirichlet boundary conditions for the hyperbolic part which is obtained by considering the unknown (u − 1) instead of u.Define then the function h by , the pair (u, p) will be called an entropy solution of (1.1)- (ii) p ∈ L ∞ (0,T;ᐂ) satisfies the weak formulation (3.3) for almost all t ∈]0, T[.
In the above definition, p k is the unique solution of (3.3) when u is replaced by the constant k and γ(u) represents the trace of u in L ∞ (∂Ω×]0,T[) and in L ∞ (Ω) for t = 0.This trace is well defined for BV functions (cf.LeRoux [20]).And finally sg is the usual sign function.
This definition is inspired by that of BLN [6] where it is given for a general flux function f = f (x,u,t).In our case, the corresponding flux is h(u)∇p and by taking the entropy pair (U,F) as we recover the definition of [6] with h(u)∇p replaced by f (x,u,t).Also in the same context of a porous medium flow, a definition is given for such a function in [14] but where m is taken as a constant.Consequently this simplication led to the decoupling of the problem P u from the problem (3.3).

A 1D Riemann problem
A first simplification in the 1D case is that the incompressibility condition div V = 0 implies the independence of V with respect to the space variable x so that we have As the total velocity at x = 0 is equal to q, we therefore end up with Let us then consider the scalar conservation equation with an initial condition of Riemann type The pressure p is then obtained by solving Let us point out that equations of type (4.3)-(4.4)can be completely solved by using hyperbolic techniques (cf., e.g., [15]).But before doing so, let us recall the following results.

A. Abbassi and G. Namah 651
The proof is omitted as it results from tedious but straightforward calculations.In fact, if α < 1, then f changes concavity from concave to convex whereas for α > 1, f is first convex before becoming concave as we go from u = 0 to u = 1.Now let be a line of discontinuity (a shock).Denote by u + (resp., u − ) the value of the solution on the right (resp., left) of .Then the speed (c) of the shock is given by the Rankine-Hugoniot conditions: We then have the following.
Proposition 4.3 (Oleinik).u is an entropy solution of (2. 16)- (2.18) if and only if one of the following three conditions is verified.
Note that if f is strictly convex, u will be an entropy solution if and only if the shock is decreasing, that is, u + < u − .Now we solve for the different cases.
(i) α = 1, λ ≤ 1.In that case, f is affine or strictly convex.Note that the initial discontinuity is admissible and it propagates with speed c (cf. (4.7)) given in our case by The corresponding solution, at time t > 0, is then given by The flow function f being strictly concave, the initial discontinuity is not admissible, so that u − and u + will be connected by a rarefaction wave.The solution at some time t > 0 is given by (cf. Figure 1.2), given by The saturation u * , which corresponds to the point of change of concavity of the flow function f , is known as the Welge saturation in oil reservoir simulations.In the interval (u * ,1), f is strictly concave so that two states u * and 1 can be connected by a rarefaction wave whereas the states u * and 0 will be connected by a shock wave (Proposition 4.3) which propagates with speed We therefore end up with a rarefaction wave followed by a shock (cf.(iv) α < 1.And finally the last case corresponds to a function f which is concave then becomes convex.By again considering the upper concave envelope, we define u * which satisfies Then the states 1 and u * can be connected by a shock (Proposition 4.3) at x = f (u * )t A. Abbassi and G. Namah 653 followed by a rarefaction wave given by This completes the resolution in the 1D case.Note that the value α > 1 seems to be the most realistic one.In that case, we end up with the following result concerning the speed of the oil-water interface.
Theorem 4.4.Let the flow function f be given by Then the water saturation is given by the profile of Figure 4.1 and the speed of the oil-water interface is a constant given by For example, when α = 2 (a value typically considered for models with a nonlinear permeability), the speed is given by In fact, we use the previous resolution for the case α > 1 where u * , given by (4.11), can be explicitly computed and is found to verify (4.17).The discontinuity in our case happens to be the oil-water interface so that the latter propagates with the speed c given by the Rankine-Hugoniot condition and the proof is done.
The profile of p. Recall that the pressure p is given for any time Its regularity therefore depends on that of the function x → m(u(x,t * )) which may be discontinuous at the shock point x = x * .Thus p is a decreasing function and has a C 1 regularity except at most at the shock point x = x * where the post-and preshock values  of the derivative are given, respectively, by where u * is as defined by (4.17).For example, for α = 2, we have so that one has m(u * ) = m(0) if λ = 3.And for λ greater (resp., less) than 3, m(u * ) is greater (resp., less) than m(0).Thus, for λ = 3, p happens to be a smooth decreasing function in the x variable.Figure 4.2 shows a typical profile for the pressure for some time t > 0.

The pseudo-2D case
Consider now a rectangular domain Ω = (0,X) × (0,Y ) with which comes to injecting water on the whole left boundary and recovering the oil on the whole right boundary.We call this configuration the pseudo-2D case because it can be assimilated to a one-dimensional flow as is stated in the following result.
Proof.Let us first remark that, for u = u 1 , p 1 defined by (5.2) is the unique solution of the problem (3.3).Indeed, for ϕ ∈ ᐂ, we have (5.3) Note also that as u is the solution obtained in the previous section by solving the problem (4.3), we have, for all t > 0, u(0,t) = 1 (5.4) and as q > 0 and f is an increasing function, it suffices to impose a boundary condition for entering data only, that is, at x = 0. Thus v 1 = u 1 − 1, where u 1 is the restriction of u on (0, X), is the entropy solution of the initial boundary value problem x ∈ (0,X). (5.5) Now it suffices to check that (v 1 , p 1 ) satisfies the inequality of the problem (P u ).Let and define also the function for all y ∈ [0,Y ], x ∈ [0,X], and t ∈ [0,T[.Such a function φ y is in D([0,X] × [0,T[).Now as v 1 is an entropy solution of the problem (5.5), it satisfies the entropy inequality 656 A two-phase interface in a porous medium in the 1D case (cf., e.g., [6]), that is, (5.8) We now subtract the null quantity and then replace q by −m(v 1 + 1)(∂p 1 /∂x) or by −m(k + 1)(∂p k /∂x).Finally integrating in y over (0, Y ), one obtains (5.10) Now it suffices to note that the third integral can be rewritten as with n = (−1,0) on Γ in and (1,0) on Γ out .And as on the rest of the boundary, the above integral can be rewritten on the whole boundary and the proof is done.

Numerical simulations
In this section, we present some numerical simulations concerning the 2D problem by using a nonconservative scheme of Harten type (cf.[1]) for the resolution of the hyperbolic problem A. Abbassi and G. Namah 657 A better resolution of the shock, that is, of the oil-water interface is expected with this scheme because no interpolation is needed for the velocity terms as is the case in the conservative Harten scheme.Indeed, in the latter case, we need to have the values of the velocity terms at the points where the values of u are calculated-for example, at the center of the grids-whereas they are computed at the interfaces of the grids.So some interpolation is necessary.The nonconservative scheme proposed calls for the values of the velocities (V ), of the saturation (u), and of the pressure (p) where they are calculated, that is, at the interfaces for V and at the centre of the grids for u and p.
Let u n i, j be the approximation of u at (x = iδx, y = jδ y, t = nδt).Define the difference operator ∆ of type then the proposed scheme reads as follows: with the functions C defined as It is formulated just as the second-order Harten scheme (cf. [18]) but with the functions ν and γ defined differently.In fact, here we have with V x the horizontal component of V and λ x = δt/δx.The same type of modification applies to the second-order terms γ.One can refer to [1] for the details concerning the derivation of the above scheme and some properties in view of its convergence.Figure 6.1 shows the profile of the saturation obtained in the pseudo-2D case for the standard set of data given in Table 6.1.
We notice the good resolution for the interface and the fact that the flow is independent of the y− variable so that it can be assimilated to a 1D flow.For this sake, whereas the analytical value, obtained by (4.19), is ,51 × 10 −6 m/s.(6.7) Other comparisons can be done to show the concordance of the numerical results obtained with the analytical ones.We end up with two profiles for the interface speed with respect to the inflow speed q (Figure 6.4) and the porosity ε (Figure 6.5).In the two cases, the computed values coincide with those expected according to the formula (4.19).

Figure 6 . 5 .
Figure 6.5.The interface speed as a function of the porosity for λ = 4.