Effective equations for two-phase flow with trapping on the micro scale

Effective equations for two-phase flow with trapping on the

micro scale

Citation for published version (APA):

Duijn, van, C. J., Mikelic, A., & Pop, I. S. (2001). Effective equations for two-phase flow with trapping on the micro scale. (RANA : reports on applied and numerical analysis; Vol. 0102). Technische Universiteit Eindhoven.

Document status and date: Published: 01/01/2001

Document Version:

Publisher’s PDF, also known as Version of Record (includes final page, issue and volume numbers)

Please check the document version of this publication:

• A submitted manuscript is the version of the article upon submission and before peer-review. There can be important differences between the submitted version and the official published version of record. People interested in the research are advised to contact the author for the final version of the publication, or visit the DOI to the publisher's website.

• The final author version and the galley proof are versions of the publication after peer review.

• The final published version features the final layout of the paper including the volume, issue and page numbers.

Link to publication

General rights

Copyright and moral rights for the publications made accessible in the public portal are retained by the authors and/or other copyright owners and it is a condition of accessing publications that users recognise and abide by the legal requirements associated with these rights. • Users may download and print one copy of any publication from the public portal for the purpose of private study or research. • You may not further distribute the material or use it for any profit-making activity or commercial gain

• You may freely distribute the URL identifying the publication in the public portal.

If the publication is distributed under the terms of Article 25fa of the Dutch Copyright Act, indicated by the “Taverne” license above, please follow below link for the End User Agreement:

www.tue.nl/taverne

Take down policy

If you believe that this document breaches copyright please contact us at:

openaccess@tue.nl

Trapping on the Micro Scale

C. J. vanDuijn

Department of Mathematics and Computing Science Eindhoven University of Technology

P.O. Box 513

5600 MB Eindhoven, The Netherlands

A.Mikelic

Laboratoire d'Analyse Numerique Universite Lyon 1

69622 Villeurbanne CEDEX, France

I. S. Pop

Department of Mathematics and Computing Science Eindhoven University of Technology

P.O. Box 513

5600 MB Eindhoven, The Netherlands

1 Introduction

A widely used technique to remove oil from reservoirs is water-drive. Through injection wells water is being pumped into the reservoir, driving the oil to the production wells. The presence of rock heterogeneities in the reservoir generally has an unfavorable eect on the recovery rate. For instance, when preferential paths (high permeability regions) exist from injection to production wells, much oil will be by-passed and consequently the oil recovery rate will be small. Conversely, when rock heterogeneities occur perpendicular to ow from injection to production wells (so-called cross-bedded or laminated structures) oil may be trapped at the interface from high to low permeability and become inaccessible to ow, leading again to a reduction in recovery rate. This latter case was studied by Kortekaas 18],van Duijnet al. 13], and more recently byvan Lingen 20] who performed laboratory experiments using a porous column with periodically varying permeability, see Figure 1. In the same context, steady state solutions as well as an averaging procedure were considered byDaleet al. 9], 10].

                         Injection Production

Figure 1: Periodically varying porous medium with high (coarse) and low (ne) permeability layers.

The main purpose of this paper is to derive in a rational way the eective ow equations corresponding to Figure 1, when the ratio of micro scale (periodicity length) and macro scale (column length) is small.

To this end we consider a one dimensional ow of two immiscible and incom-pressible phases (water being the wetting phase, oil the non-wetting phase) through a heterogeneous porous medium, characterized by a constant porosity  and a variable absolute permeabilityk=k(x). The underlying equations are:

mass balance for phases @S

@t

+@q @x

= 0 (=ow)  (1.1)

momentum balance for phases (Darcy law)

q =;k(x) k r ( S )  @p @x  (1.2)

and the complementary conditions

S o+ S w= 1  (1.3) p o ;p w= p c( xS w) : (1.4) Here S q k r

 and p denote, respectively, the saturation, specic

dis-charge, relative permeability, viscosity and pressure of phase . Throughout

we assume that the phase saturations are normalized, i.e. 0S 1.

Condi-tion (1.3) expresses the presence of two phases only. The phase pressures dier due to interface tension on the pore scale. This is expressed by (1.4), wherep

c

denotes the induced capillary pressure. In petroleum engineering it is usually described by the Leverett model, seeLeverett 19] orBear 2],

p c( xS w) =  s  k(x) J(S w)  (1.5)

wheredenotes the interfacial tension andJ the Leverett function. The relative

permeabilitiesk r : 0

1]! 01) and the Leverett functionJ: (01]! 01)

are assumed to be smooth generalizations of power law functions (Corey 8],

A1: k r strictly increasing in 0 1] withk r (0) = 0 A2: J(0+) =1,J(1)>0 andJ 0 <0 in (01],

where the prime denotes dierentiation.

Here we explicitly assume J(1)>0. Physically this means that a certain

pres-sure, the capillary entry pressure given byp c(

x1), has to be exerted on the oil

before it can enter a fully water saturated medium.

Equations (1.1) and condition (1.3) imply that the total specic dischargeq:= q

o+ q

w is constant in space. Throughout this paper we consider it constant

in time as well. With q>0 given, equations (1.1), (1.2) and conditions (1.3),

(1.4) can be combined into a single transport equation for one saturation only. Since we are primary interested in the oil ow, we use the oil saturation for that purpose. In doing so, it is convenient to redene k

rw p c and J in terms of S o. Setting u=S o ( S w= 1 ;u) we now write k rw( u) := k rw(1 ;u) p c( xu) := p c( x1;u) J(u) := J(1;u): In terms of uassumptions A 1;2 become ~A1: ( k rw strictly decreasing in 0 1] withk rw(1) = 0 k ro strictly increasing in 0 1] withk ro(0) = 0 ~A2: J(1;) =1 J(0)>0 andJ 0 >0 in 01):

Remark 1.1

In most cases of practical interest the blow-up of J and J 0 near u = 1 is balanced by the behavior of k

rw near that point, in the sense that k

rw( u)J

0(

u)! 0 as u! 1. The consequence of this behavior and its possible

failure is investigated by van Duijn&Floris 12]. Though important for the well-posedness of the mathematical formulation, no additional assumptions are required for the purpose of this paper.

Applying the scalings

x:=x=L x  t:= tq L x andk:=k=K  (1.6) where L

x is a characteristic macroscopic length scale and

K a characteristic

permeability value, we nd for the oil saturation the balance equation

@u @t

+@F @x

where F =f(u);N c k(x)(u) @ @x p c( xu): (1.7b) Here f(u) = k ro( u) k ro( u) +Mk rw( u) (1.8) denotes the oil fractional ow function, and

(u) =k rw( u)f(u) p c( xu) =J(u)= p k(x): (1.9)

The two dimensionless numbers involved are the capillary number N

c and the

viscosity rationM. They are given by N c=  p K  w qL x andM =  o  w : (1.10)

Remark 1.2

(i) Assumptions ~A1;2 imply

f(0) = 0 f(1) = 1 and f strictly increasing in 01] (0) =(1) = 0 and (u)>0 for 0<u<1:

(ii) Depending on the specic application, the value of the capillary number may vary considerably. For instance, adding surfactants or polymers may substan-tially alteror

w. Likewise the ow rate

qcan have di erent values. Therefore

we investigate in Section 2 the consequences of having a moderate and a small value forN

c.

(iii) Petroleum engineers dene the capillary number (1.10) in the reciprocal way, i.e. N c=  w qL x  p

K. Here we do not adopt this convention because we want

to emphasize the direct proportionality between the capillary number and the interface tension.

Two typical capillary pressuresp c =

p c(

xu) are shown in Figure 2. They relate

to ne (k=k ;) and to coarse ( k=k +) material, where k ; <k +.

We consider equation (1.7) in the domain  =R and fort>0, subject to the

initial condition

u(x0) =u 0(

x) forx2: (1.11)

Whenkis constant andu 0:  ! 01] is such that R u0 0 (s)J 0( s)dsis uniformly

Lipschitz continuous in , Problem (1.7), (1.11) admits a unique weak solution

u:  01)! 01]. This follows from the work of Alt &Luckhaus 1],

van Duijn & Ye 14], Gilding 15], 16], or Benilan & Toure 3]. This weak solution is smooth whenever u 2 (01) and has the usual regularity for

p c u 1 u* 0 fine k = k− coarse k = k+

Figure 2: Dimensionless capillary pressure in terms of oil saturation: top (bot-tom) reects ne (coarse) material.

degenerate equations across possible free boundaries near u = 0 and u = 1.

Whenkis piecewise constant, in particular k(x) = ( k +  x<0 k ;  x>0 (1.12) equation (1.7) cannot always be interpreted across the interface wherek is

dis-continuous. This is due to a possible discontinuity in capillary pressure. Using a regularization procedure, this was demonstrated by van Duijnet al 13] for equation (1.7), and rigorously proven by Bertsch et al 5] for a simplied equation. Instead one considers equation (1.7) only in the sub-domains where

kis constant, together with matching conditions acrossk-discontinuities. Fork

given by (1.12), with k ;

<k

+, the matching conditions read for all t>0: (i) F(t)] = 0 (1.13) (ii) u(0+t) p c( t)] = 0 p c( t)]0 (1.14) where F(t)] = F(0+t);F(0;t) and p c(

t)] likewise. The rst condition

expresses continuity of ux and is obvious. The second conditions tells us that the capillary pressure is only continuous if both phases are present on both sides of the k-discontinuity. This is to be expected from Darcy law (1.2) since then

both phase pressures are continuous. If oil is absent for x>0, i.e. in the ne

With reference to Figure 2, the pressure condition (1.14) can be formulated as (ii') 8 > < > : u(0;t)<u implies u(0+t) = 0 u(0;t)u implies J(u(0;t)) p k + = J(u(0+t)) p k ;  (1.15) whereu  is uniquely dened by J(u ) p k + = J(0) p k ; : (1.16)

The occurrence of oil trapping at the transition from coarse to ne material is directly explained by conditions (1.13), (1.15). Let k be given by (1.12) and

consider a steady state ow (u=u(x)) satisfying

u(1) = 0 (1.17)

i.e. injection and production of water, with oil possibly present near x = 0.

Then, by (1.7), F = constant = 0 onR or f(u)  1;N c p k(x)k rw( u)J 0( u) du dx  = 0 onRnf0g

with (1.15) atx= 0. Sincef(u)>f(0) = 0 for u>0, we have u= 0 or du dx = 1 N c p k(x)k r w( u)J 0( u) >0 (1.18) forx2Rnf0g. Sinceu(+1) = 0, we nd u(x) = 0 for allx>0 (1.19) and, by (1.15), u(0;)u  : (1.20)

Using (1.20) as initial condition for (1.18) on (;10), one easily constructs a

family of non-trivial steady states describing the saturation of the trapped oil in the coarse material. The initial condition in an actual displacement process determines which of the steady states is selected. This is discussed byBertsch

et al 5].

Note that the non-uniqueness results from (1.19). Considering u(1) = ^u 2

(01], one nds a unique steady state satisfying (1.15) (continuity of pressure)

atx= 0. Such solutions were considered byYortsos&Chang 25] for smooth k.

We now turn to the problem with micro structure, as indicated in Figure 1, where trapping occurs at all transitions from high to low permeability. As a

result we expect to nd a trapping related threshold saturation (irreducible oil saturation) below which the oil becomes immobile. We consider the case of a periodic as well as a random medium.

In Section 2 we assume a periodic micro structure of coarse (k=k

+) and ne

(k = k

;) material, each of length L

y << L

x, see Figure 3. This leads to a

natural choice of the small expansion parameter " = L y =L x. We outline the k 0 2L y x −L y −2L y Ly k+ k−

Figure 3: Periodic permeability (unscaled coordinate).

homogenization procedure, study the resulting auxiliary problems and derive the eective (up-scaled or averaged) equations for the limit"&0. In doing so,

the magnitude of the capillary numberN

cis important. We work out two cases:

Capillary limit,N

c = 0(1). This case is relative straightforward because the

auxiliary problem only has constant state solutions (compare steady state solu-tions on (;"0)). As a consequence the eective equation is found explicitly. It

is again of convection-diusion type and it involves weighted harmonic means of the fractional ow and capillary terms. Both convection and diusion vanish from the equation if the averaged oil saturation drops below 1

2 u

.

Balance,N c=

O("). This case is much more involved. Now the diusion term

disappears in the homogenization procedure and one is left with a rst order conservation law of Buckley-Leverett type. This follows from a detailed study of the auxiliary problem. We show that the upscaled oil-fractional ow function is dierent from ak-weighted version off and contains quite surprisingly some

satu-ration drops below a certain value. This irreducible oil satusatu-ration is related to a specic solution of the auxiliary problem.

In Section 3 we consider the case of a random micro structure with respect to both the location of the permeability jumps and the value of the permeability. The eective oil ux is obtained only for the capillary limit (N

c =

O(1)) and

involves again the weighted harmonic means of the fractional ow and capillary pressure terms. The homogenized equation has coe cients depending on the realization, but we prove that average saturation, dened by the homogenized parabolic problem, is a deterministic function. Consequently, it is su cient to solve the eective equation for a single realization.

Section 4 contains some numerical results. There we take power law relative permeabilities and a Brooks-Corey capillary pressure. We compute the eective fractional ow and diusivity for the capillary limitN

c=

O(1), and the eective

fractional ow for the balanceN c=

O(").

Some concluding remarks are given in Section 5.

Daleet al. 9] studied a similar multi-phase ow problem. They consider steady state ow in a periodic porous column, allowing each periodicity cell to have more sub-layers with dierent relative permeabilities and Leverett functions. Without using the homogenization ansatz they derive upscaled expressions for the relative permeabilities. In this paper we present a rigorous analysis of the auxiliary problems, resulting in a fairly complete description of the upscaled equations. In particular, the eect of microscopic trapping, as a result of the dierent entry pressures, is investigated explicitly.

2 Homogenization procedure for periodic layers

A simplied version of Problem (1.7), (1.11), involving a single permeability dis-continuity (or trap) only, was studied by Bertschet al 5]. They established existence and uniqueness of a solution satisfying the usual porous-media equa-tion regularity away from the trap. In particular the soluequa-tion is nonnegative and bounded. Moreover the corresponding ux was shown to be continuous in

x, for almost allt>0.

In this paper we silently assume the same properties for the saturation and ux in the case of multiple traps at arbitrary distances. In particular 0 u 1.

In our problem we deal with two natural length scales: a macroscopic length scaleL

xand a microscopic scale (characteristic length scale of layers) L

y. This

disparity in length scales is what provides us with our expansion parameter

" = L y

=L

x. For xed, but small characteristic layer length L

y the solutions

will in general be complicated having a dierent behavior on the two length scales. Closed-form solutions are unachievable and numerical solutions will be nearly impossible to calculate. It is our object to derive a ow equation at

the macro scale, keeping information about the trapping only through some averaged quantities.

To simplify our considerations we now suppose a periodic structure of the traps being located at the pointsf"i:i2Zg. The corresponding permeabilityk

"( x) is dened by k "( x) =k(x="), where k= ( k + on (2 i;12i) k ; on (2 i2i+ 1): (2.1) Without loss of generality we assume 0 <k

; <k

+

<1. We distinguish two

kinds of matching conditions: one going fromk + to

k

; and vice versa, see also

(1.15). Atx= 2i"we impose: ifu(2i";0)<u   thenu(2i"+ 0) = 0 ifu(2i";0)u   then J(u(2i";0)) p k + = J(u(2i"+ 0)) p k ; : (2.2) Atx= (2i+ 1)"we impose: ifu((2i+ 1)"+ 0)u   then J(u((2i+ 1)"+ 0)) p k + = J(u((2i+ 1)";0)) p k ;  ifu((2i+ 1)"+ 0)<u   thenu((2i+ 1)";0) = 0: (2.3) We now replacekbyk

"in equation (1.7a-b). Clearly this equation holds in the

domain "= RnT ", where T "= " S i2Z i. Letu

"be a solution of (1.7a) satisfying

the matching conditions (2.2) and (2.3). Using the uniform L

1 bound for u

",

we consider the following two scale asymptotic expansion

u "( xt) =u 0( xyt) +"u 1( xyt) +" 2 u 2( xyt)::: (2.4)

where y = x=" represents the fast scale. Substituting this expansion into

equation (1.7) and equating terms of the same order of ", gives equations for u

0 u

1

:::. As established for many linear problems containing periodic

non-homogeneities, see for instance, Bensoussan, Lions &Papanicolaou 4] or

Sanchez-Palencia 23], we expect that

U(xt) = 12 +1 Z ;1 u 0( xyt)dy (2.5)

is the weak limit of u

"and that u 0( x x " t) is the approximation tou " in some

problems in non-homogeneous geometries poses di culties as is shown in Hor-nung 17] and Mikelic 21]. Given the nonlinear nature of equation in (1.7) and matching conditions in (2.2) - (2.3), we shall therefore make no attempt at proving convergence as "! 0. The purpose of this paper is merely to derive

the upscaled equations and to study the corresponding auxiliary problems. Clearly our results depend strongly on the scaling of the capillary numberN

c.

The main cases of interest areN c =

O(1) and N c =

O("). We will deal with

them separately.

2.1 Capillary limit:

N c

=O(1)

Introducing the oil ux

F "= f(u ") ;N c p k "( x)D(u ") @u " @x  (2.6) where D(u ") = k rw( u ") f(u ") J 0( u ")  (2.7)

equation (1.7a) becomes

@u " @t +@F " @x = 0 in " (01): (2.8)

We now apply expansion (2.4) toF

", which gives F "= ;N c D(u 0) @u 0 @y p k" ;1 +f(u 0) ;N c p k  D(u 0) @u 1 @y +@u 0 @x +D 0( u 0) u 1 @u 0 @y  +  f 0( u 0) u 1 ;N c p k  D(u 0) @u 2 @y +@u 1 @x +D 0( u 0) u 1 @u 1 @y +@u 0 @x + D 00( u 0)( u 1)2 2 +D 0( u 0) u 2 @u 0 @y  "+O(" 2) =:F 0 " ;1+ F 1+ F 2 "+O(" 2) : (2.9)

Using this in (2.8) results in the following equations:

" ;2: ;N c @ @y p kD(u 0) @u 0 @y = 0 thus, by continuity ofF ", ;N c p kD(u 0) @u 0 @y =F 0= F 0( xt) (2.10)

which holds for every xy2 R and for allt >0: Note that this observation is

expected because of the continuity of the ux.

" ;1 : 0 = @F 0 @x +@F 1 @y = @ @x  ;N c p kD(u 0) @u 0 @y  + @ @y  f(u 0) ;N c p k  D(u 0) @u 1 @y +@u 0 @x +D 0( u 0) u 1 @u 0 @y  (2.11) " 0: @u 0 @t +@F 2 @y +@F 1 @x = 0: (2.12)

We look for y-periodic solutions of (2.10) satisfying (2.2) and (2.3), with x

and t as given parameters. Our goal is to prove that F

0 = 0. We argue by contradiction. SupposeF 0 <0. Let w(y) :=J(u 0( y)) (w) :=k rw( J ;1( w))f(J ;1( w)) and !(w) = w Z J(0) (s)ds

the last function being strictly increasing. Then (2.10) reads

(w) p k dw dy =; F 0 N c =:F>0: Hence, for ;1<y<0, !(w(0;));!(w(;1 + 0)) = F p k +  giving w(0;)w(;1 + 0) + F p k + 1 jjjj 1 : (2.13) Similarly, for 0<y<1, w(1;0)w(0+) + F p k ; 1 jjjj 1 : (2.14)

Now we apply matching conditions (2.2) and (2.3) in terms ofw. First, suppose w(0;)J(u ). Then w(0+) =J(0) and, by (2.14), w(1;0) >J(0). Hence w(;1 + 0)>J(u ) giving by (2.13) -w(0;)>J(u

), which contradicts the

assumption. Next suppose w(0;) > J(u

). In this case we obtain

p k ; =k + w(0;)<w(0;). By (2.14) and (2.13) we have w(1;0)  r k ; k + w(0;) + F p k ; 1 jjjj 1  r k ; k + w(;1 + 0) + F jjjj 1 ( p k ; k + + 1 p k ; ) : (2.15) If w(;1 + 0) > J(u ), then w(;1 + 0) = q k + k ; w(1;0). Substituting this

into (2.15) yields w(1;0) >w(;1;0), which contradicts the periodicity. If w(;1 + 0)  J(u

), then

w(1;0) = J(0), which contradicts (2.14). Hence F

0

0. A similar argument givesF 0

0, implying F

0= 0 :

This conclusion allows us to solve equation in (2.10) with the matching condi-tions. We nd u 0( y) = ( C>u  for ;1<y<0 C for 0<y<1 (2.16) whereC=J ;1 q k ; k + J(C) , or u 0( y) = ( Cu  for ;1<y<0 0 for 0<y<1: (2.17) Now we consider the"

;1-equation (2.11). Since F

0= 0 and the ux is

contin-uous we nd F 1= F 1( xt):

Using (2.16) and (2.17), the local form ofF 1is F 1= f(C);N c p k + D(C)  @C @x +@u 1 @y   (2.18) for;1<y<0, and F 1= 8 < : f(C);N c p k ; D(C)  @C @x +@u 1 @y  forC>u   0 forCu   (2.19) for 0 <y < 1. Clearly we only have to consider the non-trivial case C >u

.

From (2.18) and (2.19) we deduce

@u 1 @y = 8 > > < > > : f(C);F 1 p k + N c D(C) ; @C @x =:B 1( xt) for ;1<y<0 f(C);F 1 p k ; N c D(C) ; @C @x =:B 2( xt) for 0<y<1:

After integration we observe that B 1+

B

2 = 0. Hence we can solve for F 1 to nd F 1= f(C) p k + D(C) + f(C) p k ; D(C) 1 p k + D(C) + 1 p k ; D(C) ;N c @C @x +@C @x 1 p k + D(C) + 1 p k ; D(C) :

Finally we use the"

0- equation in (2.12). Since F

2is continuous in the fast scale,

we nd for the averaged oil saturationU = 1 2(

C+C) the eective

convection-diusion equation @U @t + @ @x  F(U);N c D(U) @U @x  = 0 (2.20)

where ;1<x<1andt>0. One easily veries F(U) = 8 > > > < > > > : 0 for 0U  1 2u  

strictly increasing for 12u  <U <1 1 forU = 1 and D(U) = 8 > > > < > > > : 0 for 0U  1 2u   >0 for 12u  <U <1 0 forU = 1:

In Section 4 we show the graphs of F andD based on power law relative

per-meabilities and a Brooks-Corey capillary pressure.

2.2 Balance:

N c =O(") WritingN c:= N c

", the oil ux (2.6) becomes F "= f(u ") ;N c " p k "( x)D(u ") @u " @x : (2.21)

Clearly expansion (2.9) changes due to the additional"factor. It now takes the

form F "= f(u 0) ;N c p kD(u 0) @u 0 @y + f 0( u 0) u 1 ;N c p k  D(u 0) @u 0 @x +@u 1 @y +D 0( u 0) u 1 @u 0 @y  "+O(" 2) =:F 0+ F 1 "+O(" 2) : (2.22) Using this expansion in (2.8) gives

@u 0 @t + 1 " @F 0 @y +@F 0 @x +@F 1 @y =O(")

resulting in the equations: " ;1: @F 0 @y = 0

or, by the continuity ofF ", f(u 0) ;N c p kD(u 0) @u 0 @y =F 0= F 0( xt) (2.23)

which holds for everyxy2Rand for allt>0 " 0: @u 0 @t +@F 0 @x +@F 1 @y = 0: (2.24)

First equation (2.23) needs to be considered. It leads to the following auxiliary problem. Problem A u: Given F 2R, nd u: ;10)(01]! 01] satisfying f(u);N c p kk rw( u)f(u)J 0( u) du dy =F in (;10)(01) (2.25)

subject to the matching condition (y= 0) 8 > < > : ifu(0;)<u   thenu(0+) = 0 ifu(0;)u   then J(u(0;)) p k + = J(u(0+)) p k ; : (2.26) and the periodicity condition (y=1)

8 > < > : ifu(;1 + 0)<u  then u(1;0) = 0 ifu(;1 + 0)u  then J(u(;1 + 0)) p k + = J(u(1;0)) p k ; : (2.27) This problem is considered in detail in the next sections. We prove existence for 0F 1 and uniqueness forF >0. Moreover, we show monotone

depen-dence and dierentiability of uwith respect to F. After that equation (2.24)

is averaged over the cell (;10)(01) to obtain the upscaled (macroscopic)

transport equation. This equation turns out to be of Buckley-Leverett type.

2.3 Auxiliary problem

To simplify the analysis we introduce, as in Section 2.1, the function w=J(u)

and set (w) =k rw( J ;1( w)) and '(w) =f(J ;1( w)):

In terms ofw, the auxiliary problem A

ProblemA w: Given F 2R, ndw: ;10)(01]! J(0)1) satisfying '(w)  1;N c p k (w) dw dy  =F in (;10)(01) (2.28)

such that (aty= 0) 8 > < > : ifw(0;)<J(u )  thenw(0+) =J(0) ifw(0;)J(u )  thenw(0+) = r k ; k + w(0;) (2.29) and (aty =1) 8 > < > : ifw(;1 + 0)<J(u )  thenw(1;0) =J(0) ifw(;1 + 0)J(u )  thenw(1;0) = r k ; k + w(;1 + 0): (2.30) We rst demonstrate existence and some qualitative properties for 0 =f(0) F f(1) = 1.

Lemma 2.1

LetF >1. Then there are no solutions to Problem A w.

Proof.

Since f is strictly increasing we have F '(w) ;1 F f(1) ;1>0

and consequently, by (2.28), dw=dy < 0 on (;10) (01). Now suppose w(0;) < J(u

). Then

w(0+) = J(0) and thus w < J(0) on (01),

con-tradicting w  J(0) from the denition. If w(0;)  J(u

), then clearly w(;1 + 0)>w(0;)J(u ) yielding w(0+) = r k ; k + w(0;) w(1;0) = r k ; k + w(;1 + 0):

This impliesw(1;0)>w(0+), contradictingdw=dy<0 on (01).

Lemma 2.2

LetF <0. Then there are no solutions to Problem A w.

Proof.

Equation (2.28) now givesdw=dy>0 on (;10)(01). Now suppose w(0;)J(u ). Then w(0+) =J(0) andw(1;0)>J(0). Hence w(;1 + 0)> J(u ), contradicting w(0;)J(u ). Next let w(0;)>J(u ). Then w(0+) = q k ; k + w(0;) and w(1;0) >w(0+) = q k ; k + w(0;)> q k ; k + J(u ) = J(0). Thus w(;1+0)J(u

) and, from the

w-monotonicity, q k ; k + w(0;)> q k ; k + w(;1+0) or w(0+)>w(1;0), contradictingdw=dy>0 on (01).

Corollary 2.3

LetF= 1. Then u= 1 uniquely solves ProblemA u.

Proof.

We use the u-formulation in ProblemA

u. Clearly

u= 1 is a solution.

To show uniqueness, suppose there exists a solution u such that u(y 0)

< 1

for some y 0

2 (;10)(01). Since du=dy < 0 whenever u < 1, we have

the following two possibilities. Either we have u< 1 everywhere and strictly

decreasing, or there exists y 1 < y 0 such that u(y 1) = 1 : The rst possibility

leads to a contradiction using the monotone relations imposed by the matching conditions, sinceu(0+)>u(1) implies u(0;)>u(;1). The second possibility

impliesu(y) = 1 for ally y

1, in particular

u(;1) = 1, which contradicts the

periodicity.

Lemma 2.4

LetF = 0. Then ProblemA

u admits the family of solutions (for

all0lu ): (u(y)) = 8 < :  y N c p k + + (l)  + for ;1<y<0 0 for0<y<1 where (s) = s Z 0 k rw( v)J 0( v)dv:

Proof.

Equation (2.27) implies that any solution must be a combination of

u0 and d dy (u(y)) = 1 N c p k : (2.31)

One immediately deduces that u(y) = 0 for 0 <y <1 is the only possibility.

Any other choice contradicts the periodicity. Then clearly u(0;)  u  and

(2.31) provides the required structure.

Now we consider the case 0 < F < 1. To understand the structure of the

solutions of ProblemA

wwe rst introduce

Denition 2.5

GivenF 2(01), let 0(

F)2(J(0)' ;1(

F)) be the unique root

of 0(F) Z J(0) V(sF)ds= 1 N c p k ;  (2.32) whereV(F) : (J(0)' ;1( F))(' ;1( F)1)!R + is given by V(sF) = (s)'(s) jF;'(s)j : (2.33)

−1 0 1 y w w(0 )− J(0)

ϕ

−1 (F) ξ(F) 0 w(0 )+ k+ k− w(0 )+ =

Figure 4: Sketch of behavior of solutionw.

Clearly 0(0+) = J(0)  0 2C 1((0 1)) andd 0

=dF >0 forF >0. We are now

in a position to prove the following structure for solutions of ProblemA w, see

also Figure 4.

Proposition 2.6

Let0<F <1. Then any solution of ProblemA

w satises: (i) dw dy <0 on (01), with 0( F)w(0+)<' ;1( F) (ii) dw dy >0 w>' ;1( F) on (;10).

Proof.

By a uniqueness argument for equation (2.28) we note that either

w ' ;1(

F) or w6= ' ;1(

F) on the intervals (;10) and (01). Furthermore w7'

;1(

F) impliesdw=dy70. Using this monotonicity and conditions (2.29),

(2.30), the resultw(0+)<' ;1(

F) follows directly, givingdw=dy<0 on (01).

Integrating equation (2.28) on (01) gives w(0+) Z w(1) V(sF)ds= 1 N c p k ; :

Sincew(1)J(0) we nd w(0+) Z J(0) V(sF)ds 1 N c p k ;  implying w(0+)  0(

F). Since w(0+)> w(1), conditions (2.29), (2.30) give w(0;)>w(;1), proving the second statement of the proposition.

We shall now demonstrate the solvability for Problem A

w. We start with the

simplest case where a solution satises w(1) = J(0) and w(0+) =  0(

F). By

Denition 2.5, such local solutions exist on (01). Using again (2.29), (2.30) we

nd for the left interval

w(;1)J(u ) and w(0;) = r k + k ;  0( F): (2.34)

By (ii) of Proposition 2.6 we need' ;1( F)<J(u )  or F <f(u ) 

for such solutions to exist. Integrating (2.28) over (;10) and using (2.34) once

more yields the nonlinear algebraic equation

q k + k ; 0(F) Z w(;1) V(sF)ds= 1 N c p k +  (2.35) whereq k + k ;  0( F)> q k + k ; J(0) =J(u ) :

If this equation can be solved for w(;1) 2 (' ;1(

F)J(u

)) we have found a

solution of Problem A

w satisfying

w(1) = J(0). To investigate the solvability

we dene, for 0F<f(u ), G(F) = q k + k ;  0 (F) Z J(u  ) V(sF)ds: (2.36a)

One easily veries

G(0) = 0 G(f(u )) =

1anddG=dF >0 on (0f(u ))

:

Hence there exists a uniqueF  2(0f(u )) such that G(F ) = 1 N c p k + : (2.36b)

As a consequence, integral equation (2.35) can be uniquely solved, provided 0 < F  F

: the left hand side in (2.35) decreases monotonically in

w(;1),

becomes unbounded when w(;1)&'(F) and attains a value

1 N c p k + when w(;1)%J(u

). Thus we have shown, see also Figure 5,

Theorem 2.7

Let0<F F  <f(u ), where F  is dened by (2.36b). Fur-ther, let  0(

F) be given by Denition 2.5. Then Problem A

w admits a solution w satisfying w(1) =J(0) w(0+) = 0( F) andw(0;) = r k + k ;  0( F): −1 0 1 y w J(0) ϕ−1 (F) J(u )* k+ k− ξ(F) 0 ξ(F) 0 0 < F < F* −1 0 1 y w J(0) −1 0 1 y w f(u ) < F < 1* F < F < f(u )* * ϕ−1 (F) J(u )* ϕ−1 (F)

Figure 5: Sketch of behavior ofwfor the three ranges of F.

Next we consider F  < F < 1. Since now G(F) > 1 N c p k +  there are no

solutions possible in the classw(1) =J(0). For convenience we introduce  :=w(1)2(b'

;1(

where b= max n J(0) p k ; =k + ' ;1( F) o and z:=w(0 +) 2(' ;1( F)):Below

we construct solutions satisfying w(1) > b and w(;1) > J(u

). Then the

problem of existence for Problem A

w is reduced to the system of algebraic

equations (integrating (2.28) on (;10) and (01), and using (2.29), (2.30) and

(2.33)): z Z  V(sF)ds= 1 N c p k ;  (2.38a) q k + k ; z Z q k + k ;  V(sF)ds= 1 N c p k + : (2.38b)

To study the solvability of this system we introduce

: (b' ;1( F))  ' ;1( F) r k + k ; ' ;1( F) ! !R by (v) = 8 > > > < > > > : v R b V(sF)ds forb<v<' ;1( F) q k + k ; ' ;1 (F) R v V(sF)ds for' ;1( F)<v< q k + k ; ' ;1( F):

Note that  is strictly increasing on ;

b' ;1(

F) 

and strictly decreasing on

' ;1( F) q k + k ; ' ;1( F)

, see Figure 6. By the monotonicity of, the function z=z() = ;1 () + 1 N c p k ; (2.39a) is well-dened on (b' ;1(

F)), satisfying dz=d > 0. Now system (2.38a-b)

reduces to study the mapW : (b' ;1( F))!R, given by W() =  r k + k ;  ! ;  r k + k ; z ! ; 1 N c p k + : (2.39b)

We rst formulate the theorem.

Theorem 2.8

For F 

< F <1, there exists a solution to (2.38a-b), i.e. the

auxiliary Problem A

b ζ v

Ψ

ϕ

−1 (F)

ϕ

−1 (F) k+ k− z k− N_c 1

{

Figure 6: Sketch of and construction ofz=z().

Proof.

Since z(' ;1( F);) =' ;1( F) we have W(' ;1( F);) =; 1 N c p k + <0:

To investigate the behavior near  =b, we usez> 0( F) and consider q k + k ; 0(F) Z q k + k ; b V(sF)ds= ( +1 forf(u ) F<1 > 1 Nc p k + for F  <F <f(u ) :

The rst follows fromq k + k ; b=' ;1( F) forFf(u

), the second from q k + k ; b= J(u

) and (2.36a) for F



<F < f(u

). As a consequence we nd

W() >0

for close tob. SinceW is continuous, the equationW() = 0 has at least one

root, which provides the existence for (2.38a-b).

2.4 Continuity, monotonicity and uniqueness

To construct the eective equation for U, we need to show that the solution

F 1. TheF-dependence is denoted by u=u(yF),w=w(yF), or simply u(F)w(F). We treatF 2(0F

) and

F 2(F 

1) rst, and then consider the

behavior nearF = 0+F=F  and

F = 1. F 2(0F

)

Since uniqueness has not yet been demonstrated, we consider here the solution

w(F) given by Theorem 2.7. It satises w(yF) Z J(0) V(sF)ds= 1 ;y N c p k ; for 0 <y1 (2.40a) q k + k ; 0(F) Z w(yF) V(sF)ds=; y N c p k + for ;1y<0: (2.40b) The smoothness of  0 and V(s) implies w(y) 2 C 1((0 F )) for each y 2

;10)(01]. Let(F) = dw=dF. Dierentiating (2.40a) with respect toF

yields ; w(yF) Z J(0) (s)'(s) (F;'(s)) 2 ds+V(w(yF)F)(yF) = 0: Hence (yF)>0 for 0<y<1 (2.41) and (0+F) = d 0 dF >0 (1F) = 0: From (2.40b) we nd q k + k ;  0 (F) R w(yF) (s)'(s) ('(s);F) 2 ds+V  r k + k ;  0( F)F ! r k + k ; d 0 dF = V(w(yF)F)(yF) implying (yF)>0 for ;1y<0 (2.42) with (0;F) = r k + k ; d 0 dF >0:

F 2(F 

1)

Then any solution of ProblemA

w satises w(yF) Z w(1F) V(sF)ds= 1 ;y N c p k ; for 0 <y1 withw(1F)>J(0). Hence ; w(yF) R w(1F) '(s) (s) (F;'(s)) 2 ds+V(w(yF)F)(yF) = V(w(1F)F)(1F) (2.43) implying the statements:

if(1F)>0 then(yF)>0 for 0<y<1 if(0+F)<0 then(1F)<0: (2.44a) Similarly we deduce on (;10): if(0;F)>0 then(yF)>0 for ;1<y<0 if(;1F) = 0 then(0;F)<0: (2.44b) The conditions at y= 0andy=1 translate into

(0;F) = r k + k ; (0+F) (;1F) = r k + k ; (1F) (2.45) Next we combine (2.44a) and (2.45). Suppose there exists ^F 2(F



1) such that (1F^) = 0. Then(;1F^) = 0(0;F^)<0(0+F^)<0 giving(1F^)<0,

a contradiction.

Hence either (1F) > 0 or (1F) < 0 for all F 2 (F 

1). We rule out

the second possibility. By (2.45), (1F) < 0 gives (;1F) < 0, implying

that w(;1F) is strictly decreasing in (F 

1). However Proposition 2.6 gives w(;1F)>'

;1(

F)!1asF !1, a contradiction. Hence(1F)>0 and by

(2.44a)

(yF)>0 for y2 ;10)(01]: (2.46)

Remark 2.9

Note that the monotonicity result (2.46) applies to any solution of Problem A

w satisfying

w(1F)>J(0). We use this to show uniqueness for

Problem A

w and hence for Problem A

Theorem 2.10

The auxiliary problem(A

u) has a unique solution

u(F) for each F 2(01]. We have i) u(1) = 1 ii) u(F) =J ;1( w(F)) wherew(F) is given by w(yF) Z w(1F) V(sF)ds= 1 ;y N c p k ; for 0<y1 q k + k ; w(0+F) Z w(yF) V(sF)ds=; y N c p k + for ;1y<0 with w(1F) =J(0) w(0+F) =  0( F) for 0 <F F , and w(1F) >J(0) satisfyingW(w(1F)F) = 0 for F  <F <1.

Proof.

In Section 2.3 we have shown that forF 

< F <1 no solutions are

possible with w(1F) = J(0). Furthermore for 0 < F  F

, Problem A

w is

uniquely solvable in the class w(1F) = J(0). What remains is to rule out

solutions satisfyingw(1F)>J(0) for 0<F F

 and to show uniqueness for F



<F<1 in the classw(1F)>J(0).

With W given by (2.39b), et us consider the equation

W((F)F) = 0 with(F) =w(1F)>J(0):

Dierentiating with respect toF and denoting@=@ by a prime gives W 0 d dF +@W @F = 0:

Since@=@F>0, as explained in Remark 2.9, we have W

0(

(F)F)<0 (2.47)

whenever@W=@F >0. The denition ofW involves z=z(F), given by (zF) =(F) + 1 N c p k ; : Hence  0( zF) @z @F = @ @F ((F);(zF))

implying@z=@F>0. Using this we nd directly @W @F = @ @F   r k + k ; F ! ; r k + k ; zF !! ; r k + k ;  0  r k + k ; zF ! @z @F >0:

Thus (2.47) holds for any solution of (A w) with

(F) =w(1F)>J(0).

Next we considerW(bF). In Section 2.3 we showedW(bF)>0 forF >F  andW(' ;1( F)F) =; 1 Nc p k ; <0. In fact, forF <f(u ) we have W(bF) = q k + k ; 0(F) R J(u  ) V(sF)ds; 1 N c p k ; =G(F); 1 N c p k ; (see 2.36a). (2.48) Hence W(bF) = 8 > > < > > : >0 forF >F   0 forF =F   <0 forF <F  :

Combining these inequalities with (2.47) gives uniqueness forF >F

 and

non-existence forF F .

Letu: ;10)(01]! 01], dened by (see Lemma 2.4) (u(y)) = 8 < :  y N c p k + + (u )  + for ;1y<0 0 for 0<y 1

denote the maximal solution corresponding toF = 0.

We are now in a position to formulate the following continuity and monotonicity results.

Theorem 2.11

The solutionu(F) satises:

(i) u()2C 1((0 F ) (F  1)) and @u @F (F)>0 on ;10)(01], except for 0<F <F  where @u @F (1F) = 0 (ii) lim F%1 u(yF) = 1 (iii) lim F%F  u(yF) = lim F&F  u(yF) =u(yF ) (iv) lim F&0 u(yF) =u(y).

Proof.

Monotonicity follows directly from the previous results. Therefore we only need to demonstrate the continuity properties (ii)-(iv).

(ii) By Proposition 2.6 we have

w(yF)>' ;1( F) for ;1y<0 and consequently w(yF)w(1F) = r k ; k + w(;1F)> r k ; k + ' ;1( F) for 0 < y  1 and F > F . Since ' ;1( F) ! 1 as F % 1, the uniform convergence ofu(F) follows.

(iii). The result forF %F

 is a direct consequence of the continuity of 

0( F).

To establish the result for F & F

, we consider the function

W(F) for F

nearF  and

 near b=J(0). Direct computation shows W 0( bF) =; r k + k ; f(u ) k rw( u ) f(u ) ;F <0: (2.49) Since W(F) and W 0(

F) are uniformly continuous in f(F) : b    b+F



F F +

gforsu ciently small, we use (2.48) and (2.49) to nd (F) =w(1F)&J(0) asF &F

 :

The uniform convergence on both intervals now follows from the w(yF)

ex-pressions in Theorem 2.10.

(iv). The uniform convergence in (01] results from 0(

F)&0 asF &0. To

establish the result in ;10) we note that monotonicity and boundedness of u(F) imply lim F&0 u(yF) = ~u(y) pointwise in ;10) with ~u(0;) =u . Moreover, since 0<N c p k + k rw( u)f(u)J 0( u) du dy =f(u);F <1

on ;10)u(F) is uniformly continuous inF. Hence, by Dini's Theorem, the

convergence is uniform in ;10) and ~u 2 C( ;10)). Let y 0

2 ;10) with

~

u(y 0)

>0. ForF >0, the integral equation foru(F) can be written as (u(0;F));(u(yF)) + F u(0;F) R u(yF) k rw( s)J 0( s) f(s);F ds=; y N c p k + : Lety=y 0. Then, for F su ciently small, 0<F u(0;F) Z u(y0F) k rw( s)J 0( s) f(s);F ds < F Const u(0;F) Z ~ u(y0) 1 f(s);F ds!0

as F&0. Hence (u ) ;(~u(y 0)) = ; y N c p k +  implying ~u(y 0) = u(y 0).

2.5 The e ective equation

Letu=u(F) denote the unique solution of ProblemA

u. As in Section 2.2 we write F 0= F 0( xt) and set u 0( xyt) =u(yF 0( xt))

forx2Ry2 ;10)(01] andt>0. The equation for the averaged saturation U(xt) = 12 1 Z ;1 u 0( xyt)dy

results from (2.24). Integrating this equation in y and using the continuity of F 1( xt) we nd @U @t +@F 0 @x = 0 for x2R t>0: (2.50)

From here on we drop the superscript and write F =F

0. As a consequence of

Theorem 2.11 we note that the cell-averaged saturationU =U(F) satises U 2C( 01])\C 1((0 F ) (F  1)) with dU dF >0 on (0F ) (F  1): Moreover, U(0+) =U U(1) = 1 where U = 12 0 Z ;1 u(y)dy:

The continuity and monotonicity allows us to dene the inverseF: 01]! 01]

satisfying, withF(U ) = F , F 2C( 01])\C 1(( UU ) (U  1)) and dF dU >0 on (UU ) (U  1):

Further,

F(U) = 0 for 0U U andF(1) = 1:

Taking F = F(U) as nonlinear ux function in (2.50), results in an eective

equation which is a rst order conservation law forU, with U as macroscopic

irreducible oil saturation.

Under additional (but usual) assumptions onk ro

k rwand

J we show that

equa-tion (2.50) is of Buckley-Leverett type in the following sense.

Theorem 2.12

For o  w >1 and >0, let k ro( s) s o = O(1) k rw( s) (1;s) w = O(1) and (1;s)  J(s) =O(1): Then F 2C 1( 0 1]) (implying F 0( U) = 0) andF 0(1) = 0 :

Proof.

We rst consider the behavior near U =U. Writing equation (2.40a)

in terms ofu=J ;1(

w) and dierentiating with respect toF yields @u @F = F;f(u) k rw( u)f(u)J 0( u) u Z 0 k rw( s)f(s)J 0( s) (F;f(s)) 2 ds:

We now use equation (2.25) twice to rewrite this expression into

@u @F =; @u @y 1 Z y 1 F;f(u(sF)) ds:

Next we integrate iny. Setting U +( F) = R 1 0 u(yF)dy anda(F) =J ;1(  0( F)) we nd dU + dF = 1 Z 0 a(F);u(sF) F;f(u(sF)) ds> 1 F (a(F);U +( F)): Thus d dF (F U +( F))>a(F) implying U +( F)> 1 F F Z 0 a(s)ds for 0<F F  : Since U(F)>U+ 12U +( F)

we have U(F)>U+ 12 F F Z 0 a(s)ds for 0<F F  : (2.51)

We need to estimate a(F) =u(0+F) from below. For this we use Denition

2.5, i.e. a(F) Z 0 k rw( s)f(s)J 0( s) F;f(s) ds= 1 N c p k ;  which gives 1 F;f(a(F)) a(F) Z 0 k rw( s)f(s)J 0( s)ds> 1 N c p k ;  and further 0<F;f(a(F))<Ca(F)f(u(F)) for 0<F <F  

where C(here and below) denotes a suitably chosen positive constant.

Now using f(s)=s 0 =

O(1) (implied by the asymptotic behavior of k ro) we nd, for smallF a(F)>CF 1= 0 :

Combining this with (2.51) gives

F(U)<C(U;U) 0

in a right neighbourhood of U.

Next we consider the dierentiability of F(U) at U = U . For

F < F  we

use (2.40a-b). Dierentiating the equations with respect to F and using the

continuity of w(yF) gives directly the existence of (yF  ;) for each y 2 ;10)(01]. For F >F  we rst observe that (1F) is bounded in a right neighbourhood ofF

. This follows from the proof of Theorem 2.11 (iii). Hence,

in equation (2.43)

V(w(1F)F)(1F)!0 asF &F 

and thus, again using (2.43), (yF +) =

(yF 

;) fory 2 (01). A similar

argument holds in (;10). As a consequence,F is dierentiable atU .

To prove F

0(1) = 0, we construct an upper bound for

U(F) nearF = 1. For ;1<y<0 we have, as in (2.43), w(0;F) Z w(yF) (s)'(s) ('(s);F) 2 ds +V(w(0;F)F)(0F) = V(w(yF)F)(yF):

Hence @u @F > f(u);F k rw( u)f(u)J 0( u) u(0;F) Z u(yF) k rw( s)f(s)J 0( s) (f(s);F) 2 ds

which can be written as

@u @F > @u @y 0 Z y 1 f(u(sF));F ds Consequently,U ;( F) = R 0 ;1 u(yF)dysatises dU ; dF > 0 Z ;1 u(sF);u(;1F) f(u(sF);F ds> 1 1;F fU ; ;u(;1F)g which implies U ;( F)< 1 1;F 1 Z F u(;1s)ds: (2.52)

Next we estimate u(;1F) from above nearF = 1. Sinceu(1F)<f ;1(

F),

the periodicity condition implies

u(;1F)<J ;1  r k + k ; J(f ;1( F)) ! : Using 1;f(s) (1;s) w = O(1) and (1;s)  J(s) =O(1) we nd u(;1F)<1;C(1;F) 1  w near F= 1:

Substituting this estimate into (2.52) and usingU +(

F)<1we deduce F(U)>1;C(1;U)

w

in a left neighbourhood ofU = 1.

3 Randomly layered media in the capillary limit

In this section we drop the periodicity assumption and suppose a stationary ergodic geometrical structure. It is characterized by a probability space ("),

with an ergodic dynamical systemT(x)x2R(see, e.g. 22] or 11] for details).

For a-measurable subsetP " we introduceP =P(!)Rby

and we call it a random stationary set.

In our application we suppose thatP(!) has the following form P(!) =  i2Z (y 2i;1  y 2i)  (3.2)

where the random variablesy i

2Rare strictly increasing with respect toi.

A representative example is a Poisson process # inRwith constant rate >0.

In this case the number of points of # in an intervalA= (ab) has expectation (b;a). The number of points of # in any bounded interval is then nite with

probability 1 and # has no nite limit points. On the other hand the number in (0+1) is innite so that the points in (0+1) can be written in order as

0<y 1 <y 2 <y 3 <:::

Similarly the points in (;10) can be written in order as :::y ;3 <y ;2 <y ;1 <0:

These exhaust the points of # since the probability that 0 2 # is equal to 0. y

n are random variables and the subsequences fy

n

n ;1gand fy n

n 1g

are independent, with the same joint distributions. Furthermore, the random variables` 1= y 1, ` n= y n ;y n;1( n2),` ;1= ;y ;1, l ;n= y ;n+1 ;y ;n( n

2) are independent and each has probability densityg(y) = e ; jyj

:The number

of pointsN(0t] of # in (0t] satises the law of large numbers

lim

t!+1

1

t

N(0t] = with probability 1:

Finally, the process of Poisson is ergodic. Another example are hardcore pro-cesses (Gibbs propro-cesses, Mat$ern propro-cesses, :::). We construct them from a

Poisson point process by eliminating all points having a distance to its neigh-bors smaller than a prescribed value. They satisfy the mixing property and the ergodicity is assured.

By Birkho's Ergodic Theorem there exists a density (fraction of high perme-ability layers) of P, given by

' +:= (P) = lim NM!1 1 y 2M+1 ;y ;2N;1 M X ;N jy 2i( !);y 2i;1( !)j (3.3)

for almost all!2", satisfying

0' +

1:

The corresponding random permeability is given by

k(x!) =k(T(x)!) = ( k +( !) ifx2P k ;( !) ifx2RnP (3.4)

y k+(ω) k−(ω) k+(ω) k−(ω) k+(ω) 2i−1(ω) y 2i(ω) y 2i+1(ω) y 2i+2(ω) y 2i+3(ω) y 2i+4(ω) y

Figure 7: Random distribution of layers and it is a stationary random variable. Then

k "( x!) =k(T(x=")!): (3.5) As a consequenceu = u ( !) through (1.16).

Next we turn to the two scale expansion for the saturation and the ux, adapted to the stochastic case. We write

F "= " ;1 F 0+ " 0 F 1+ "F 2+ :::  (3.6) whereF

k are stationary ergodic random elds and u "= u 0+ "u 1+ " 2 u 2+ ::: : (3.7)

From these expansions and (2.8) we obtain directly

dF 0 dy = 0 implyingF 0= F 0( xt) (3.8)

with the random variableu

0satisfying (2.10). We reconsider this equation for a

given realization!, see Figure 7. As before we want to showF

0= 0. Suppose F

0

<0. Introducing wandas in Section 2.1 we obtain again (w) p k dw dy =; F 0 N c =:F >0: (3.9)

We argue below that this inequality does not permit us to construct a global non-negative solution satisfying the matching conditions at the interfaces. Sup-pose w(y

2i

;0)  J(u

). By (3.9), implying strict monotonicity of w, we have w(y 2i;1+ 0) < J(u ) giving w(y 2i;1

;0) =J(0). This contradicts the

monotonicity of w in (y 2i;2

y

2i;1). Next suppose w(y 2i ;0) > J(u ). Then w(y 2i+ 0) = r k ; k + w(y 2i ;0) andw(y 2i+1 ;0)w(y 2i+0)+ F jjjj 1 1 p k ; jy 2i+1 ; y 2i j. Therefore we have w(y 2i+1+ 0) = r k + k ; w(y 2i+1 ;0)  w(y 2i;1 + 0) + F jjjj 1 ( p k + k ; jy 2i+1 ;y 2i j+ 1 p k + jy 2i ;y 2i;1 j ) .

Repeating this reasoning backwards in i shows that w will drop below J(u )

at the right side of a certain transition, yielding again a contradiction. Hence

F 0

0. A similar argument givesF 0

0, soF 0= 0.

This implies foru 0, with i2Z, u 0( y!) = 8 > < > : C(!)>u  for y 2i;1( !)<y<y 2i( !) J ;1 q k ; k + J(C(!)) fory 2i( !)<y<y 2i+1( !) or u 0( y!) = ( C(!)u  for y 2i;1( !)<y<y 2i( !) 0 fory 2i( !)<y<y 2i+1( !) :

Now consider the "

;1-equation (2.11). Since F 0 = 0, the ergodicity of F 1 impliesF 1= F 1( xt) which is given by F 1= f(u 0) ;N c p k(!)D(u 0) @u 1 @y +@u 0 @x : SupposeC(!)u . Then F 1= 0 on ( y 2i;1( !)y 2i( !)) impliesF 1 = 0 for all x2R andt>0. If C(!)>u , then we have on ( y 2i;1( !)y 2i( !)) @u 1 @y = f(C(!));F 1 p k +( !)N c D(C(!)) ; @C(!) @x =:B 1( !): On (y 2i( !)y 2i+1( !)) we have @u 1 @y = f(C(!));F 1 p k ;( !)D(C(!)) ; @C(!) @x =:B 2( !) withCas in (2.16). Since @u 1 @y

is the local representation of a stationary random variable with zero mean, we have that the mean value of

 fk =k + g B 1( !) + fk =k ; g B 2( !)

is zero. Heredenotes the characteristic function. Hence ' + p k +( !) f(C(!));F 1 N c D(C(!)) + 1;' + p k ;( !) f(C(!));F 1 N c D(C(!)) = ' + @C(!) @x + (1;' +) @C(!)) @x : Solving forF

1gives (dropping the

!-dependence) , forC(!)>u , F 1= ' + p k + f(C) D(C) + 1;' + p k ; f(C) D(C) ' + p k + 1 D(C) + 1;' + p k ; 1 D(C) ;N c ' + @C @x + (1;' +) @C @x ' + p k + 1 D(C) + 1;' + p k ; 1 D(C) : (3.10)

Averaging (2.12) and using the ergodicity of F

2 yields the eective transport

equation @U @t +@F 1 @x = 0 for ;1<x<1 t>0 (3.11)

whereU denotes the averaged oil saturation U ='

+

C+ (1;' +)

C:

We note that for' +

u 

U >0,F

1= 0. In the periodic case '

+ = 1

2. Hence

for each realization!we obtain an equation of the `periodic' form (2.20).

We stress that the averaged saturationU is deterministic. This is implied by

equation (3.11) and by the deterministic initial condition. Consequently, it su ces to consider only one realization to determine U and its corresponding

ux F 1.

Remark 3.1

In the case of the balance N c =

O("), we could proceed

analo-gously. Now we should solve the problem (2.25) - (2.26) on the real line, for every realization. The periodicity condition (2.27) is replaced by the condition thatutakes values between0 and 1 onR. We note that the matching condition

is now posed at every point y i

 i2Z. Solving the auxiliary problem A

u in the

stochastic case is much more complicate than in the periodic case. The analysis of the periodic case was already quite lengthy and in the proofs of Proposition 2.6 and Theorem 2.7 periodicity was essential. Also unboundedness of J

com-plicates proofs. Using arguments from this section we are able to conclude that

F 2 01], but the complete construction is still an open problem. We expect to

consider randomly layered media in the limitN c=

O(") in a future publication.

4 Numerical results

Since we have no convergence proof, we are going to verify the homogenization procedure numerically for the periodic case. Both the capillary limit and the balance will be considered. We will use the Leverett model with Corey relative permeabilities and Brooks-Corey capillary pressure ( 8], 6]). Specically, the following functions and parameters are used:

k ro( u) =u 2  k rw( u) = (1; u) 2  J(u) = (1; u) ; 1 2  M = 1 k + = 1  k ;= 0 :5 withN c being either 1 or

", depending on the case.

Tests are done on the interval (;11), i.eL

x= 1. For both cases we compute the

full problem with a periodic micro-structure as shown in Figure 3 i.e. k(x) =k +

in the coarse layers and k(x) =k

; in the ne layers. The thickness L

y of the

layers is related to the number of cells and determines the expansion parameter

"=L y

=L

x. The matching conditions dened in (1.13) and (1.14) or (1.15) are

solution is averaged on each micro cell consisting of two adjacent layers. This average is compared with the numerical solution of the eective equations. In the tests we consider a medium originally saturated by oil (u(x0) = 1x2

(;11)), with water injection from the left (u(;1t) = 0). At x = 1 the

Neu-mann condition is chosen not to aect the ow. For the full problem 80 micro-structure cells are considered, implying"= 1=80.

4.1 Capillary limit (

N c

=O(1)

)

In this case we takeN

c= 1. Inside each layer of constant permeability we apply

a rst order explicit discretization scheme with upwind nite volumes. Withu n i

denoting the approximate oil saturation att n=

n inside the volume centered

at x i= (

i;1=2)h( being the time step andhthe grid size), the solution at

the next time-step follows from

u n+1 i = u n i ;  h  F n i+1=2 ;F n i;1=2  : (4.1) Here F n

i+1=2 approximates the ux at t=t

n and x =x

i+

h=2 = ih, the edge

between the volumes centered in x i and x i+1. Likewise, F n i;1=2 approximates the ux at t=t n and x=x i ;h=2 = (i;1)h.

Ifiis such that x=ihlies inside a homogeneous micro-layer, the computation

of F n

i+1=2 is straightforward. Recalling the notation introduced in (1.8) and

(1.9), we have F n i+1=2= f(u n i) ;N c p k(x i) D u n i + u n i+1 2 u n i+1 ;u n i h  (4.2)

where the permeabilityk(x

i) is either k

+ or k

;, depending on the type of the

material. The diusion coe cient is given byD(u) =(u)J 0( u) and is calculated at the mean ofu n i and u n i+1.

Computing the ux at a position where the permeability and saturation are discontinuous requires more attention. Let us assume that this position is lo-cated at x=ih, thus separating the control volumes centered in x

i and x i+1. Moreover, letk(x i) = k +and k(x i+1) = k

;. As in 7] and 13] we introduce two

sets of dummy variables atx i: u n i+ and u n i;for all n= 012:::. They satisfy

the pressure condition in (1.15),

u n i; <u  implies u n i+= 0  or u n i; u  implies J(u n i;) p k + = J(u n i+) p k ;  (4.3) and they are chosen such that the numerical ux is continuous at x

Given a pairu n

i satisfying (4.3), we show below how to obtain u n+1 i and F n i+1=2. Once F n

i+1=2 is known, we use (4.1) at

i and ati+ 1 to determineu n+1 i and u n+1 i+1. In terms ofF n i+1=2 we rst write u n+1 i; = u n i; ;  h  F n i+1=2 ;F n i;1=2   u n+1 i+ = u n i+ ;  h  F n i+3=2 ;F n i+1=2   (4.4) from which we uniquely determineF

n i+1=2 and u n+1 i satisfying (4.3). To do so we rst elliminateF n

i+1=2. Summing up the equalities (4.4) yields u n+1 i+ = A;u n+1 i;  (4.5) withAdened as A=u n i;+ u n i+ ;  h  F n i+3=2 ;F n i;1=2  :

Note that theCFLrestriction implies 0A2.

p_c u* A 1 u 0 J(0) k+ J(0) k− J(A−u) k− J(u) k+ J(u) k− un+1_i− p_c u A=u * 1 u 0 J(0) k+ J(0) k− J(A−u) k− J(u) k+ J(u) k− n+1 i− Figure 8: Findingu n+1 i+ : A>u  (left) and A<u  (right).

Now we check which of the two situations in (4.3) applies. Let us assume rst that Au

. By contradictiction, since u

n+1

i satisfy the matching conditions,

we get u n+1 i+  u . Obviously, u n+1 i+  u  also implies A  u . Therefore pressure is continuous in (4.3) iAu

. In this case (4.5) gives J(u n+1 i+ ) p k ; =J(A;u n+1 i; ) p k ; = J(u n+1 i; ) p k + :

Monotonicity of J guarantees the existence of a unique u n+1 i;

 u

 satisfying

the last equality above (see also Figure 8). If A < u , since u n+1 i+  0, we have u n+1 i+ < u

 and, due to the matching

conditions, u n+1 i; = 0. As above, u n+1 i+ <u  also implies A<u . Therefore in

this case we uniquely obtainu n+1 i+ =

A<u .

A similar procedure is used at a transition from ne to coarse material. Details are omitted.

As explained in Section 2.1, the eective equation is known explicitly in the cap-illary limit. Figure 9 shows the eective diusivityDand convectionF in terms

of the cell averaged oil saturationU. Here we use the relative permeabilities and

Leverett function as proposed for this section. This equation is of degenerate

0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1

Figure 9: Eective diusion (left) and convection (right) for the Brooks-Corey model. Note thatD(U) =F(U) = 0 for 0U 1=2u

= 0 :25.

parabolic type, since the eective diusionD(U) vanishes for 0U  1 2 u

 and

at U = 1. Several numerical methods can be applied to such kind of problems.

Here we use the explicit upwind scheme (see 24] for the convergence analysis):

U n+1 i = U n i ;  h  F n i+1=2 ;F n i;1=2   F n i+1=2= F(U n i ) ;D U n i + U n i+1 2 U n i+1 ;U n i h :

Figure 10 shows the solution of the eective equation (solid line) and the av-erage of the solution of the full problem (dashed line) at t = 0:3 and t = 0:7.

Since oil is being displaced from the column, both solutions are approaching the macroscopic irreducible oil saturation corresponding to the maximum amount of trapped oil: U = 1=2u

= 0

:25. The transition region travels with the same

speed in both cases.

The solution of the original problem is shown in Figure 11, together with its average. The leftmost part of the interval is enlarged in the picture on the

0 0.2 0.4 0.6 0.8 1 -1 -0.5 0 0.5 1 effective averaged 0 0.2 0.4 0.6 0.8 1 -1 -0.5 0 0.5 1 effective averaged

Figure 10: Averaged and eective solution oil saturation att= 0:3 andt= 0:7.

0 0.2 0.4 0.6 0.8 1 -1 -0.5 0 0.5 1 oscillatory averaged 0 0.1 0.2 0.3 0.4 0.5 0.6 -1 -0.98 -0.96 -0.94 -0.92 -0.9 -0.88 -0.86 -0.84 -0.82 -0.8 oscillatory averaged

Figure 11: Full problem averaged and oscillatory oil saturation att= 0:7, full

(left) and zoomed view (right).

right. Note the good agreement with the theoretical results: the prole is highly oscillatory on the macro scale and quite at within the micro structure. Further note that even though the original problem is of degenerate type, free boundaries seldomly occur inside the homogeneous sub-layers. As a consequence the solution behaves fairly nondegenerate and thus smoothly. Therefore there is no need to consider many points inside the sub-cells. In our computations an interior grid size ofh=L

y

=10"L x

=10 was su cient to produce good results.

However, since the numerical method is explicit, the time step is subject to a CFL condition.

4.2 Balance (

N c =O(")

)

With N c =

", computations for the full problem are done exactly in the same

way as for the capillary limit. However, the eective equation requires more attention. As shown in Section 2.5, this equation is of Buckley-Leverett type, but the fractional ow function is not known explicitly. In this case a table of values for the pairs (UF) has to be constructed, whereF ranges from 0 to

auxiliary problem (A

u) dened in Section 2.2, and calculate its cell-average as

the corresponding U value in the table. For the purpose of this paper we took F

i =

i&F, with &F = 10 ;3 and

i = 010

3. As stated in Lemma 2.4 we take F(U) = 0 for all U 2 0U'], 'U being the average of the maximal steady state

solution (corresponding tol=u ).

To nd accurate solutions of Problem A

u we modify the dierential equation

through the Kirchhof transform

(u) = Z u 0 k rw( v)f(v)J 0( v)dv: (4.6)

Note that  is strictly increasing and smooth due to the properties of k rw,

f

and J. In general, the integration cannot be carried out explicitly. Therefore,

again we need to construct a table of pairs (u(u)). Here we have applied an

adaptive quadrature method. Thus, instead of solving ProblemA

u, we consider the equivalent

Problem A : Given F i, nd : ;10)(01]!R satisfying f( ;1( ));N c p k d dy =F i in ( ;10)(01) (4.7)

together with the corresponding matching and periodicity conditions dened in (2.26) and (2.27).

The matching and periodicity conditions can be viewed as boundary conditions for equation (4.7) on the two subintervals. To nd a solution u(F

i) we have

applied the following shooting procedure. Choose(1)0 and use this value as

initial condition for equation (4.7) on (01). This yields the corresponding(0+)

and, by the matching conditions, (0;). Use this value as initial condition for

(4.7) on (;10). Then adjust(1) so that(1) and(;1) satisfy the periodicity

condition. In carrying out this shooting procedure several technical di culties had to be resolved. We omit the details in this paper.

Figure 12 shows the eective oil fractional ow function for the specic model considered in this section. Observe that indeed we have recovered Buckley-Leverett model in which the fractional ow has only one inection point. Note that the theoretical analysis only resulted inF

0( U) =F

0(1) = 0. No statements

about the inection points could be given. Also note that the upscaled frac-tional ow contains details of the small scale capillary forces. This eect does not appear explicitly, but it is present due to Problem A

u (or, equivalently,

equation (4.7) with the matching and periodicity conditions). Finally notice that the macroscopic irreducible oil saturation U for the model considered is

much smaller than in the capillary limit case. This is to be expected because the capillary forces are now much smaller (O(")).

0 0.2 0.4 0.6 0.8 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

Figure 12: Eective oil fractional ow for the Brooks-Corey model. Note that

F(U) = 0 for 0U U 2:5410 ;2.

Once the eective convection is known, the oil saturation equation is solved by the rst order explicit upwind scheme

U n+1 i = U n i ;  h (F(U n i ) ;F(U i;1)) : 0 0.2 0.4 0.6 0.8 1 -1 -0.5 0 0.5 1 effective averaged 0 0.2 0.4 0.6 0.8 1 -1 -0.5 0 0.5 1 effective averaged

Figure 13: Averaged and eective solution oil saturation att= 1:0 andt= 1:5.

Figure 13 shows the solution of the eective equation (solid line) and the average of the solution of the full problem (dashed line) att= 1:0 andt= 1:5. For both

solutions we see rst a rarefaction wave, where the oil saturation approaches the maximum amount of trapped oil,U =U. This is followed by a shock. Note the

good agreement of the two solutions in the rarefaction part. The small dierence in the shock location can be explained by numerical errors which occur in the computation of the eective fractional ow function.

The solution of the full problem together with its average is shown in Figure 14, with a zoomed view in the picture on the right. Note the highly oscillatory prole on both scales. Free boundaries and large gradients occur inside almost

0 0.2 0.4 0.6 0.8 1 -1 -0.5 0 0.5 1 oscillatory averaged 0 0.1 0.2 0.3 0.4 0.5 0.6 -1 -0.98 -0.96 -0.94 -0.92 -0.9 -0.88 -0.86 -0.84 -0.82 -0.8 oscillatory averaged

Figure 14: Full problem averaged and oscillatory oil saturation att= 1:5, full

(left) and zoomed view (right).

every homogeneous sub-layer, this being a consequence of the smallness of dif-fusion (O(")). In this case the computational grid has to be quite ne to locate

the free boundaries accurately and to compute the macroscopic irreducible oil saturation. In contrast to the capillary limit, we need here many gridpoints inside each sub-layer (h=L

y

=80"L x

=80), but under a less restrictiveCFL

condition.

5 Conclusions

The results of this paper lead to the following conclusions.

 For N c =

O(1) (capillary limit) the eective equation is of degenerate

parabolic type. The diusion and convection vanish up to the macroscopic irreducible oil saturation (1

2 u

).  For N

c =

O(") (balance) the upscaled equation is of Buckley-Leverett

type, with eects of the local capillary forces in the fractional ow function.

 The macroscopic irreducible oil saturation depends strongly on the value

of the capillary number.

 The solution of the auxiliary problem in the capillary limit has two

con-stant states connected by the pressure condition at the interface.

 The solution of the auxiliary problem in the balance is unique and can be

classied completely.

 The choice of the characteristic values in (1.6) is important for deciding

which of the two cases (capillary limit or balance) applies in a real situa-tion.

 Random layers are considered only in the capillary limit. The eective

 The method used in this paper can be applied to heterogeneous media

in which the porosity, relative permeabilities and Leverett function are periodic as well.

Acknowledgements.

The authors wish to thank Dr. J. Molenaar (Shell Ri-jswijk) for the preliminary computations and for useful discussions and sugges-tions.

References

1] H. W. Alt & S. Luckhaus, Quasilinear elliptic-parabolic dierential equa-tions, Math. Z., 183(1983), 311 - 341.

2] J. Bear, Dynamics of uids in porous media, American Elsevier, New York, 1972.

3] P. Benilan & H. Toure, Sur l'equation generale u t =

a(:u'(:u) x)x+

v

dansL

1. II: Le probleme d'evolution, Ann. Inst. Henri Poincare, Anal. Non

Lineaire, 12(1995), 727 - 761.

4] A. Bensoussan & J. L. Lions & G. Papanicolaou, Asymptotic analysis for periodic structures, Studies in Mathematics and its Applications 5, North-Holland, Amsterdam, New York, Oxford, 1978.

5] M. Bertsch & R. Dal Passo & C. J. van Duijn, Analysis of oil trapping in porous media ow, to appear.

6] R. H. Brooks & A. T. Corey, Hydraulic properties of porous media, Hy-drology Papers, Colorado State University, Fort Collins, Colorado, 1964. 7] G. Chavent & J. Jaret, Mathematical models and nite elements for

reser-voir simulations, Studies in Mathematics and its Applications 17, North-Holland, Amsterdam, New York, Oxford, 1986.

8] A. T. Corey, The interrelation between gas and oil relative permeabilities, Producer's Monthly, 19(1954), 38 - 41.

9] M. Dale & S. Ekrann & J. Mykkeltveit & G. Virnovsky, Eective rel-ative permeabilities and capillary pressure for one-dimensional heteroge-neous media, Transport in Porous Media, 26(1997), 229 - 260.

10] M. Dale, Darcy equations in heterogeneous porous media: interface condi-tions and uniqueness of periodic solucondi-tions, Working Papers from Stavanger College, 43(1998).

11] D. J. Daley, D. Vere-Jones, An introduction to the theory of point processes, Springer Series in Statistics, Springer-Verlag, New York, 1988.

12] C. J. van Duijn & F. J. T. Floris, Mathematical analysis of the inuence of power-law uid rheology on a capillary diusion zone, J. Pet. Sc. Eng., 7(1992), 215 - 237.

13] C. J. van Duijn & J. Molenaar & M. J. de Neef, Eects of capillary forces on immiscible two-phase ow in heterogeneous porous media, Transport in Porous Media, 21(1995), 71 - 93.

14] C. J. van Duijn & Q. X. Ye, The ow of two immiscible uids through a porous medium, Nonlin. Anal. Theory Meth. Appl., 8(1984), 353 - 367. 15] B. H. Gilding, The occurrence of interfaces in nonlinear diusion-advection

processes, Arch. Rat. Mech. Anal., 100(1988), 243 - 263.

16] B. H. Gilding, Qualitative mathematical analysis of the Richards equation, Transport in Porous Media, 5(1991), 651 - 666.

17] U. Hornung, ed., Homogenization and porous media, Interdisciplinary Ap-plied Mathematics 6, Springer-Verlag, New York, 1997 .

18] T. F. M. Kortekaas, Water/oil displacement characteristics in crossbedded reservoir zones, Soc. Pet. Eng. J., (1985), 917 - 926.

19] M. C. Leverett, Capillary behavior in porous solids, Trans AIME Petr. Eng. Div., 142(1941), 152 - 169.

20] P. van Lingen, Quantication and reduction of capillary entrapment in cross-laminated oil reservoirs, PhD Thesis, Sub-Faculty of Applied Earth Sciences, Delft University of Technology, 1998.

21] A. Mikeli$c, Homogenization theory and applications to ltration through porous media, in \Filtration in porous media and industrial application", A. Fasano (ed.), Lecture Notes in Mathematics Vol. 1734, Springer-Verlag, Berlin, 2000, 127 - 214.

22] G. C. Papanicolaou, Diusion in random media, in \Surveys in applied mathematics vol. I", J. B. Keller, D. W. McLaughlin & G. C. Papanicolaou (eds.), Plenum Press, New York, 1995, 205 - 254.

23] E. Sanchez-Palencia, Non-homogeneous media and vibration theory, Lec-ture Notes in Physics 127, Springer-Verlag, Berlin Heidelberg New York, 1980.

24] M. Watanabe, An approach by dierence to the porous medium equation with convection, Hiroshima Math. J., 25(1995), 623 - 645.

25] Y. C. Yortsos & J. Chang, Capillary eects in steady-state ow in hetero-geneous cores, Transport in Porous Media, 5(1990), 399 - 420.