Two-phase flow in porous media : dynamic capillarity and heterogeneous media

We investigate a two-phase porous media ﬂow model, in which dynamic eﬀects are taken into account in phase pressure diﬀerence. We consider a one-dimensional heterogeneous case, with two adjacent homogeneous blocks separated by an interface. The absolute permeability is assumed constant, but diﬀerent in each block. This may lead to the entrapment of the non-wetting phase (say, oil) when ﬂowing from the coarse material into the ﬁne material. We derive the interface conditions coupling the models in each homogeneous block. In doing so, the interface is approximated by a thin porous layer, and its thickness is then passed to zero. Such results have been obtained earlier for standard models, based on equilibrium relationship between the capillary pressure and the saturation. Then, oil is trapped until its saturation on the coarse material side of the interface, exceeds an entry value. In the non-equilibrium case, the situation is diﬀerent. Due to the dynamic eﬀects, oil may still ﬂow into the ﬁne material even after the saturation drops under the entry point, and this ﬂow may continue for a certain amount of time that is proportional to the non-equilibrium eﬀects. This suggests that operating in a dynamic regime reduces the account of oil trapped at interfaces, leading to an enhanced oil recovery. Finally, we present some numerical results supporting the theoretical ﬁndings.


Introduction
Two-phase (wetting/non-wetting) flows are widely encountered in various real-life processes.A few prominent examples are water driven oil recovery, or geological sequestration of CO 2 .These processes involve large spatial scales, and (rock) heterogeneities appear naturally.In this paper we consider a simplified situation, where the medium consists of two different homogeneous blocks with different permeabilities (coarse and fine), which are separated by an interface.This makes the transition from one material to another not smooth, and appropriate conditions have to be imposed at the interface for coupling the models written in each of the two blocks.In particular, when the underlying models are involving entry pressures to describe the dependency of the capillary pressure on the phase saturations, the non-wetting phase may remain trapped in the coarse block at the interface.This situation has been analyzed in [12,14], but for the case when the phase pressure difference depends on the, say, wetting phase phase saturation and the medium itself.These are standard models, for which the dependency between various quantities are determined under equilibrium conditions.Therefore, such models are also called equilibrium models.In the paper mentioned above, regularization arguments (i.e.approximating the interface by a thin porous layer ensuring a smooth transition between the two homogeneous blocks) are employed to derive appropriate coupling conditions between the models in the two sub-domains.The resulting conditions are the flux continuity and an extended pressure condition.We refer to [2,5,6] for the mathematical analysis of such models, where the existence of weak solutions has been analyzed.Further, the case of many layers is studied in [13,10], where homogenization techniques are applied for deriving an effective model.Such kind of models are also studied in [1,4,17,18], where appropriate numerical schemes are studied.
Various experiments [3,9] have invalidated the equilibrium assumptions, and motivated considering non-equilibrium approaches.Here we focus on models involving dynamic effects in the phase pressure difference, as proposed in [19].In this case, we follow the ideas in [12,14] and derive the coupling conditions at the interface separating the two homogeneous blocks.When compared to the equilibrium case, a striking difference appears.In the former the non-wetting phase can only flow into the fine block if its saturation at the coarse block side of the interface exceeds an entry value, in the latter situation this flow can appear for lower saturation values.This is due to the dynamic effects i the phase pressure difference and can reduce the amount of non-wetting phase that remains trapped at the interface.
The models including dynamic effects in the phase pressure difference lead to so-called pseudoparabolic problems.For such models, but posed in homogeneous domains, existence and uniqueness of weak solutions are obtained in [7,16,24].The case of vanishing capillary effects and the connection to hyperbolic conservation laws is studied in [15,11].For dynamic capillarity models in the heterogeneous case, but in the absence of an entry pressure, numerical schemes are discussed in [20].This situation is similar to the case analyzed in [8], where the interface is replaced by a discontinuity in the initial conditions.Also, variational inequality approaches have been considered in [22] for situations including an entry pressure.However, the conditions are simply postulated and no derivation is presented.
In Section 2 we present the mathematical model.For simplicity we consider the case when only the absolute permeability is different in the two blocks, all other parameters being the same.In Section 3, we derive the coupling conditions at the interface.These are the flux continuity and an extended pressure continuity.In the last section, we discuss different numerical approaches, and present numerical experiments that confirm the analysis in Section 3.

Mathematical model
We consider the flow of two immiscible and incompressible phases in a one-dimensional heterogeneous porous medium.Letting s w denote the saturation of the wetting fluid, and s n the saturation of the non-wetting fluid, one has 0 ≤ s w , s n ≤ 1.The porous medium is assumed to be saturated by the two fluids, Mass balance holds for each phase (see [25,21]), where φ is the porosity assumed constant, and q α denotes the volumetric velocity of the phase α.These velocities satisfy the Darcy law where k(x) is the absolute permeability of the porous medium, p α the pressure, µ α the viscosity, and k rα the relative permeability of the α phase.The functions k rα are assumed known.Gravity effects are disregarded, as they have no influence on the interface conditions derived here.Substituting (3) into (2) gives In standard models, the phase pressure difference depends on the saturation, which is determined experimentally.An example in the sense is the Leverett relationship where σ is the interfacial tension, and J a decreasing function.
The relationship in ( 5) is determined by measurements carried out under equilibrium condition.In other words, before measuring the pressure and saturation in a representative elementary volume, fluids have reached equilibrium and are at rest.However, processes of interest may not satisfy this condition, and dynamic effects have to be included.Alternatively to (5), in [19] the following model is proposed The damping factor τ is assumed to be known and constant.Summing the two equations in (4) and using (1), one gets where q = − k(x)krw µw ∂pw ∂x − k(x)krn µ ∂pn ∂x denotes the total flow.In the one-dimension case, this means that q is constant in space.Here we assume it is constant in time as well, and is positive.This allows reducing the two-phase flow model to a scalar equation in terms of, say s = s w .After rescaling the space x with L, the time t with T , and using the reference value K for the absolute permeability and σ φ K for the pressure, one obtains where the F denotes the dimensionless flux of the wetting phase Here q = qT φL > 0, and k(x) = k(x) K .Further, f w is the fractional flow function of the wetting phase, with the mobility ratio M = µ n /µ w , the capillary number N c = σ √ φK µnqL , the dimensionless damping factor τ = τ µ n q σφ 2 and λ(s) = k rn (s)f w (s).For simplicity, we assume here N c = 1, as different values of N c would not have any influence on the conditions derived below.
We consider a simple heterogeneous situation, where two adjacent homogeneous blocks are separated by an interface located at x = 0.For the ease of presentation, we assume that all parameters and functional dependencies except the absolute permeability are the same.For the latter, we have Here k − > k + > 0 are given.Throughout this paper, the non-wetting phase may be oil, or CO 2 or any other phase with non-wetting phase properties.For simplicity below we use as oil non-wetting phase.Also, by pressure we actually mean the phase pressure difference.

Conditions at the interface
The model ( 8)-( 9) is a parabolic equation, where the factor 1 √ k appears under a second order spatial derivative.Since k has a jump discontinuity at the interface, the model is only valid in each of the two blocks, and coupling conditions at x = 0 have to be derived.Commonly, the pressure continuity is taken as a second condition.However, this is shown to be inappropriate for entry-pressure models in the absence of dynamic effects (τ = 0).This statement is made rigorous in [12] by regularizing k.Specifically, the interface is replaced by a thin layer in which k decays continuously from k − to k + .Next to k, we will use the quantity Clearly, and h − > h + > 0.
For given > 0, we approximate the interface x = 0 by the interval [− , ] (the thin layer) and h by a smooth function h , such that h is monotonic in the small interval [− , ].Specifically, the discontinuous function h is now approximated by the smooth function h , such that The function ĥ is smooth and monotonic on [− , ] (see Figure 1).With the given the solution Assuming ∂v ∂t bounded uniformly with respect to , passing to 0 gives lim In fact, this is nothing but the flux continuity at the interface, which is as expected.

Fine block
Figure 2 The strip Ω.
Along the boundary of Ω, we define The goal is to understand how s − (t) and s + (t) are related, as well as p − (t) = J(s − (t)) ∂t .These are nothing but the left and right saturation and pressure at the interface.Note that (11) is an ordinary differential equation in t, where y can be seen as a parameter.To define an initial condition, we let s 0 (y) be a smooth function s 0 : [−1, 1] → R satisfying s 0 (−1) = s − (0) and s 0 (+1) = s + (0).Clearly, the choice of s 0 is not unique.Below we investigate the relation between s − (t) and s + (t), and its dependence in the choice of s 0 and the regularization of h.Since λ(v) > 0, for 0 < v < 1, and λ(0) = λ(1) = 0, from (11), one has In other words • v(y, •) is continuous with respect to t, for t ≥ 0.
For the sake of understanding, we consider some particular cases.

Constant saturation at the coarse side of the interface
We let s − ∈ (0, 1) and assume s − (t) = s − for all t.Let s 0 : [−1, 1] −→ (0, 1) be the given initial value, not necessarily compatible to s − : s 0 (−1) = s − in general.We construct a solution for which the set Ω c = {(y, t) ∈ Ω, 0 < v(y, t) < 1} is connected.From (12), if (y, t) ∈ Ω, v solves the autonomous initial value problem: By (A1) and the assumption on h, P y has a unique solution locally.Let s * be defined by We consider the cases s − > s * and s − ≤ s * , separately.
• Case 1: Since h is decreasing function in y, there exists a unique y * ∈ (−1, 1) such that We study now the long time behavior of v(y, t).We distinguish the following sub-cases.
a) For y ∈ (−1, y * ], define s ∞ (y) as the unique solution of Clearly, s ∞ (y) is the equilibrium point for P y and satisfies s ∞ > s − (see Figure 3).Further, standard phase plane arguments show that, regardless of s 0 (y) ∈ (0, 1), lim This implies the solution v(y, t) increases in t and reaches v = 1 in finite positive time (see Figure 4) The long time behavior of v is summarized in  Proposition 1 Assume s 0 (y) ∈ (0, 1) for all y ∈ [−1, 1] and let s − > s * , where s * is defined in (14).With y * in (15), one has: Corollary 1 In particular, at y = 1, we have for s , as presented in Figure 5.This allows constructing an extended pressure condition, similar to [12].Consider the pressure observe that With the entry pressure since s − > s * , one has By Corollary 1, we obtain the condition: In other words, the pressure remains continuous for t < t * (1) and oil keeps flowing into the fine material although p − is below the entry pressure.This effect is due to incorporating dynamic effects in the phase pressure difference, and would not be possible in their absence (τ = 0).
• Case 2: s − < s * For any y ∈ [−1, 1], with s ∞ (y) defined by ( 16), one has Therefore, there exists an s ∈ (s − , 1), such that which, in view of the monotonicity of h and J, gives s ∞ < s < 1.Consequently, this case is similar to the sub-case y < y * before.
Another special case is when τ 0.Then, the time t * (y) in Proposition 1 and Corollary 1 approaches to 0, and if s − > s * the pressure becomes discontinuous instantaneously.This is, in fact, exactly the behavior in [12] for equilibrium models.

Non-constant saturation at the coarse side of the interface
The results before are obtained for a constant saturation at the coarse side of the interface.Here we generalize these results.We let p ± (t) = J(s ± (t)) ∂t be the phase pressure difference at the two sides of the interface and derive an extended pressure condition similar to (20) and (22).With the entry pressure p + e defined in (19), we assume p − (t) given, and distinguish the following cases.
) for all t > 0, and therefore the pressure remains continuous at both sides of the interface for all t > 0 As in Section 3.1, if the model involves no entry pressure (J(1) = 0), one has p + e = 0. Then p − (t) ≥ p + e = 0, and the pressure is continuous for all t > 0. • Case 2: 0 < p − (t) < p + e for all t > 0 We first assume that an exists such that 0 < p − (t) < p + e − for all t > 0. We have Proposition 4 If an exists such that 0 < p − (t) < p + e − for all t > 0, then there exists a Proof As before, let s 0 : (−1, 1) → R be such that s 0 (y) ∈ (0, 1) for all y.Then, for t small enough, v(y, t) < 1, and v(y, t) is a solution to Since p − (t) < p + e − , by the continuity of h, there exists δ > 0 such that Hence, for y ∈ (1 − δ, 1) and t > 0, we have This implies Clearly, a finite t = t(y) > 0 exists such that, v(y, t(y)) = 1, and v(y, t) for t < t(y).Letting now y = 1, there exists t * > 0, such that s + (t) = v(1, t) < 1 for 0 < t < t * , and s + = v(1, t) = 1 for t ≥ t * , and the proposition is proved.
Observe that, for fixed y ∈ (−1, 1), v(y, t) solving P y and w(y, t) solving ), t > 0, w(y, 0) = s 0 (y), are ordered.Specifically, as long as both v and w remain less 1, one has w(y, t) ≤ v(y, t).To see this, we let u = v − w, and define u − = min{0, u}.Then, subtracting the equations in P y and P y , and multiplying the result by u − gives Here we have used the monotonicity of J. Since at t = 0 one has u − (y, 0) = 0, integrating in time gives u − (y, t) = 0, implying the ordering.Hence w is a lower bound for v, and therefore analyzing w makes sense.Let y * ∈ (−1, 1) be defined by Then for y > y * , a δ > 0 exists such that h(y) = h(y * ) − δ < h(y * ), as long as w(y, t) < 1, one has This shows that w reaches the cut-off value w = 1 at finite time t = t(y).If the function s 0 (y) is constant, it is straightforward to show that t(y) is strictly decreasing to t(1) > 0 (see Figure 6).At y = y * , w = 1 satisfies the equation.Therefore, w ≡ 1 is an equilibrium solution, implying that w(y * , t) solving P y satisfies w(y * , t) < 1 for all t, and thus lim y y * t(y) = ∞.
Based on the analysis before, we discuss particular examples where oil trapping may occur.In the transition region (the blown up interface) Ω = {(y, t) : −1 < y < 1, t > 0}, we have To understand the trapping, we now take s 0 (y) = 1, for −1 < y < 1, and give the pressure at the coarse side of the interface for t > 0. We assume the following behavior (see Figure 7): a) There exists T 1 > 0 such that 0 < p − (t) < p − e for t ∈ (0, T 1 ), where p − e = J(1) h − .In this case we have , assuming that a t ∈ (0, T 1 ) exists such that v(y, t) = 1 for all t < t and v(y, t) < 1 then one has This gives a contradiction since v(y, 0) = s 0 (y) = 1.So we obtain v(y, t) = 1, for −1 ≤ y ≤ 1 and 0 < t < T 1 .
Thus y = ỹ(t) defines a free boundary in the transition region, separating regions where v = 1 from regions where v < 1.
Remark 2 In this context, for t < T 2 one has s + (t) = v(1, t) = 1, and oil remains trapped in the coarse medium at the interface.
c) There exists T 3 > T 2 such that p − > p + e for t ∈ (T 2 , T 3 ).Based on the above analysis, we have v(y, T 2 ) < 1 for all y ≤ 1, since ỹ(T 2 ) = 1.Further, one also has As before, one cannot obtain v(y, t) = 1 for some t > T 2 , since at t, it holds Therefore, v(y, t) < 1 for − 1 < y < 1 and T 2 < t < T 3 .Consequently, we have Remark 3 Next to the pressure continuity, this shows that oil starts flowing into the fine region for t > T 2 .
d) For all t > T 3 , p − < p + e .We assume p − (•) decreasing and lim , we compare the solution v(y, t) and w(y, t) of for all t > T 3 , with v(y, T 3 ) = w(y, T 3 ) < 1.

Remark 4
Compared to the equilibrium case (τ = 0), a striking difference appears.If τ = 0 for a pressure p − (t) behaving as in Figure 7, no oil flows into the fine layer for any t > T 3 .If τ > 0, oil continues flowing for t > T 3 and up to a time T * < ∞, the delay time.This delay appears, as discussed, if lim The following result extends the statement in the remark to a more general situation.
Therefore s + (t) is strictly increasing whenever s + < 1.Furthermore, we have giving for all t such that s + (t) < 1 and For convenience, define Hence, there exists t * < ∞ for which f (t * ) = 1 − s 0 .Consequently, there exists T * < t * for which s + (T * ) = 1, and since s + is increasing, one has s 0 < s + (t) < 1 for all t < T * .

Comparison of extended pressure conditions with static case
In this section, we show the difference in the pressure conditions appearing in the equilibrium and non-equilibrium models between static case and dynamic case.As proved in [12], if s − ≥ s * , then one has and no oil flows into the fine medium.Further, if s − < s * , then s + < 1, but the pressure continuous: [p] = p − (t) − p + (t) = 0.
Then oil flows into the fine medium.Conversely, given p − (t) = J(s − (t)) the pressure at the interface on the coarse side and assuming that ), and therefore s − ≥ s * .In this case, one also has s + = 1.
In other words, a delay (T * − T 3 ) appears in the non-equilibrium case before trapping occurs.
In fact, this delay can be infinite an extreme situation, when s(t) < 1 for all t > 0, can be constructed.To see this, we assume J : (0, 1] → R locally Lipschitz, and study the behavior of with appropriately chosen p − satisfying p − (t) < p + e for all t > 0. First, note that an L > 0 exists such that for all s ∈ [s 0 , 1] Hence, an upper bound to s + is the solution u of This gives and the upper bound for s + reads we obtain s + (t) < 1 for all t > 0, and consequently, the pressure remains continuous for all t > 0, whereas oil flows into the fine layer.

Numerical method and examples
In this section, we provide some numerical examples to illustrate how the dynamic effects influence the flow of the oil across the interface between two homogeneous blocks.For simplicity, in ( 8)-( 9), we take φ = 1.This gives ∂s ∂t for t > 0, and x ∈ (−l, l).The boundary and initial conditions are given below.

Linear numerical scheme
For the discretization of ( 30) and ( 31) we decompose the interval (−l, l) into 2N + 1 cells : where the grid is uniform with ∆x = 2l 2N +1 , we let x j = j∆x for j ∈ {−N − 1/2, ..., N + 1/2}.The discontinuity of h(x) is at x = 0.With ∆t > 0 a given time step, the fully discrete scheme is: where s n−1 i−1/2 is the approximation of s(x, t) at x = x i−1/2 and at t = t n = n∆t.Since q > 0, if i = 0 the upwind flux F n i at x = i∆x and t = t n is defined as Here h +/− means h − if i < 0, or h + if i > 0. At i = 0, we introduce two saturation unknowns s −,n , s +,n and define the F −,n , and F +,n as At the interface we also define the left and right discretized pressures By using the extended pressure condition discussed before one has and (32) implies either s +,n = 1, or the pressure continuity Obviously, ∂ 1 g > 0 and ∂ 2 g < 0. Further, given s − ∈ (0, 1), one has lim , since g(s − , •) is strictly increasing and continuous.For any C n−1 ∈ (−∞, g(s − , 1)] there exists a unique s Also, note that g(s − , 1) is decreasing in s − , lim and therefore g(•, 1) : (0, 1] → [g(1, 1), +∞) is one to one.Recalling that g(s − , s + ) = C n−1 appears, in fact, in the discretized extended pressure condition (32), for given s − and C n−1 > g(s − , 1) when pressure becomes discontinuous at the interface, one considers the pair (s − , 1).In this way, we have actually constructed the curves in the (0, 1] × (0, 1] square: where D(•) is the inverse of g(•, 1).Below we give a property of the discretized extended pressure condition: Proof If p −,n = p +,n , one has s +,n = 1, and obviously, If p −,n = p +,n , then we have In this case, if s +,n ≥ s −,n , one has Together with (35), we yield which concludes the proof.

Fully implicit scheme
Here a nonlinear, implicit scheme is discussed as an alternative to the linear one.Next to improved stability properties, for this scheme we can prove that s ±,n , the saturation at the interface, remain between 0 and 1, a property that is not guaranteed for the linear scheme.To construct the scheme we first define the decreasing function β : R → R by and rewrite the flux in (31) as As before, for i = 0 we write but now the upwind flux F n i becomes As before, by h ± we mean h − if i < 0, or h + if i > 0. Further, if i = 0, the flux is defined on each side of the interface as . For having a conservative scheme, the two expressions should be equal.Combined with the pressure condition (32), and viewing s n ±1/2 as well as the saturation values s n−1 i+1/2 and s ±,n−1 as known, this provides a nonlinear system with s ±,n as unknowns.Below we show that this system has a unique solution pair in the square [0, 1] 2 .
The condition F −,n = F +,n can be written as where R(s −,n , s Obviously, R is increasing in both arguments and one has Note that, if both s −,n−1 , s +,n−1 take the values 0 or 1, the last terms in (39) vanish, giving Since β is decreasing, in this case one has 0 ≤ B ≤ R(1, 1).Further, if s −,n−1 is not equal to 0 or to 1, with ∆t small enough one gets and analogously for s +,n−1 .From (39), this shows again that 0 ≤ B ≤ R(1, 1).This gives the following: Lemma 1 For a sufficiently small time step ∆t, with C n−1 = C(s −,n−1 , s +,n−1 ) defined in (33) and for R and B in (37) and (39), the system has a unique solution pair (s −,n , s Proof The set Γ (C n−1 ) introduced in (34) contains pairs satisfying the pressure condition (32).
If such a pair exists, it solves the system (40).Since as g is increasing in the first argument and decreasing in the second one, long as both s −,n and s +,n are below 1, the curve Γ (C n−1 ) is a graph of a strictly increasing function (see Figure 9).Similarly, for B ∈ [0, R(1, 1)], the set R(•, •) = B is the graph of a decreasing function in the s −,n − −s +,n plane.Therefore, these two curves have at most one intersection point inside the square [0, 1] 2 , implying that (40) has a unique solution pair.
Note that if the solution pair provided above lies inside [0, 1) 2 , then one has pressure continuity across the interface, p −,n = p +,n .Otherwise, the solution pair lies on the boundary of the square [0, 1] 2 .Moreover, assuming that initially one has s −,0 ≤ s +,0 , implying that at the interface separating the two blocks more oil is present at the coarse material side then at the fine material side, repeating the proof of Proposition 6 one obtains that s −,n ≤ s +,n for all n.This means that pressure discontinuity at the interface can only occur if no oil is present at the fine material side, s +,n = 1, which is precisely the discrete pressure condition in (32).

Remark 6
The construction above assumes that s n −1/2 , s n 1/2 are known.In fact, these are part of the solution computed implicitly, at time step t n .This means that, actually, s n −1/2 , s n 1/2 and consequently B depend on s −,n , s +,n .However, decoupling the calculation of the solution pair s ±,n from the effective time stepping suggests an iterative procedure for the implicit scheme: using, say, the values s ±,n−1 as starting point, compute s n i+1/2 by solving (36) away from the interface, and use s n −1/2 , s n 1/2 to update s ±,n .

Numerical results
For the numerical calculations, we used the following functions and parameters: The tests are done in the interval (−1, 1).Further, we present the results in terms of the oil/nonwettig phase saturation s o = 1 − s, as this is the phase for which trapping may occur.The initial oil saturation is hat shaped (see Figure 10) < x < −0.34, 0.9, −0.34 ≤ x ≤ −0.12 0, −0.12 < x < 1.
At x = ±1 we take homogeneous boundary conditions, ∂ x s(±1, t) = 0.This mimics the situations when an oil blob in the coarse layer is displaced by water.After a certain time, the oil reaches the interface.Note that initially no oil is present in the fine medium.Before discussing the results we recall that the saturation s * defined in ( 14) is the limit saturation allowing for pressure continuity in the equilibrium models (τ = 0).This also defines an entry saturation for the oil, s entry = 1 − s * .For the equilibrium model, oil flows into the fine material only if s o > s entry at the coarse material side of the interface, and it remains trapped if s o ≤ s entry .Figure 11, displaying the results obtained for τ = 0 at t = 0.7, confirm this statement.At the coarse side of the interface, one has s − o < s entry (picture on the left) and the oil flux is 0 there (picture in the middle).This means that no oil enters into the fine medium.Further, the pressure is discontinuous at the interface (picture on the right).
The case τ = 1, presented in Figure 12, shows a different situation.In the left picture, although s − o < s entry , one still has s + o > 0, meaning that oil has already entered in the fine medium.This is also confirmed by the middle picture, displaying a non-zero the oil flux at the interface.Finally, the picture on the right shows that the pressure is continuous.Figure 13, presents the results for τ = 0 and at t = 4. Then oil has flown into the fine medium (picture on the left).The flux and the pressure are both continuous (middle and right pictures).At the same time, but with τ = 1, we observe that more oil has flown into the fine media (left picture of Figure 14).As expected, the flux and pressure are continuous as well.The results above suggest that the amount of oil flowing into the fine material increases with τ .To understand this behaviour we compare the oil saturation obtained for τ = 0, τ = 1 and τ = 10, all at the same time t = 4.The profiles in Figure 15 show that, for τ = 0, little oil has flown into the fine media, and this amount is higher for τ = 1.In both cases, s − o , the oil saturation at the coarse side of the interface, already exceeds the entry saturation s entry .However, for τ = 10, s − o < s entry , but oil has still flown into the fine material.On expects that the oil flow into the fine material will take longer for the largest value of τ , in agreement with Corollary 1.
Finally, we observe that in the equilibrium case τ = 0 one can determine the maximal amount of oil that can be trapped at the interface, see [12].Having this in mind, we choose again a hatshaped initial saturation where the total amount of oil equals the maximal amount that can be trapped for equilibrium models.With this initial data, we compute the numerical solutions for three values of τ , namely 0, 10 and 30. Figure 16 shows the results at t = 400, when practically all solutions have reached a steady state and no oil flow is encountered anymore.The left picture shows the result for τ = 0.In this case, the entire amount of oil is trapped in the coarse medium, and the oil saturation s − o is matches the entry saturation s entry .No oil has flown at all into the fine material.As following from the middle picture, for τ = 10, one has s − o < s entry , and the oil remaining trapped in the coarse material is less than in the equilibrium case.The situation becomes more obvious for the solution corresponding to τ = 30.The oil saturation s − o has decayed further, and the trapped oil is less than in the previous cases.This is again in agreement with the analysis in Section 3.2.

Conclusions
We have considered a non-equilibrium model for two-phase flow in heterogeneous porous media, where dynamic effects are included in the phase pressure difference.A simple situation is considered, where the medium consists of two adjacent homogeneous blocks.We obtain the conditions coupling the models in each of the two sub-domains.The first condition is, as expected, flux continuity, whereas the second is an extended pressure condition extending the results in [14] for the standard two-phase flow model.
In the equilibrium case, if an entry pressure model is considered, oil can flow into the fine material only if its saturation exceeds an entry point.In the non-equilibrium case instead, the non-wetting phase may flow even if the oil saturation at the coarse side of the interface is below the entry point, and amount of oil remaining trapped at the interface is less than in the case of equilibrium models.Finally, two different numerical schemes are discussed , and different numerical experiments are presented to sustain the theoretical findings.

Figure 1
Figure 1 The function h .

Figure 4
Figure 4 The solution for s − > s * .

Figure 6
Figure 6 The function t(y).

Figure 8
Figure 8 Saturation inside the thin layer appximating the interface.

Figure 10
Figure 10 Initial oil saturation s 0 o .