341 Transport in Porous Media 5: 341-379. 1990. Kluwer Academic Publishers. Printed in the Netherlands. © 1990 Two-Phase Flow in Heterogeneous Porous Media I: The Influence of Large Spatial and Temporal Gradients MICHEL QUINT ARD and STEPHEN WHITAKER * Laboratoire Energetique et Phenomenes de Transfert, VA CNRS 873, Ecole Nationale Superieure d'Arts et Metiers, 33405 Talence Cedex, France (Received: 28 February 1989, revised: 14 September 1989) Abstract. In order to capture the complexities of two-phase flow in heterogeneous porous media, we have used the method of large-scale averaging and spatially periodic models of the local heterogeneities. The analysis leads to the large-scale form of the momentum equations for the two immiscible fluids, a theoretical representation for the large-scale permeability tensor, and a dynamic, large-scale capillary pressure. The prediction of the permeability tensor and the dynamic capillary pressure requires the solution of a large-scale closure problem. In our initial study (Quintard and Whitaker, 1988), the solution to the closure problem was restricted to the quasi-steady condition and small spatial gradients. In this work, we have relaxed the constraint of small spatial gradients and developed a dynamic solution to the closure problem that takes into account some, but not all, of the transient effects that occur at the closure level. The analysis leads to continuity and momentum equations for the {3-phase that are given by a{E{3}* at + v. {(vll)} = 0, {(VIl)} = "'!"K$. (v{(p/3)/3}Il- P/3g) - u/3 a{E/3}* - U/3' v Il-Il at 1 1 1 - -.illl : vv{(p/3)Il}/3 - - ｾＱｬＺｶ＼＾ＯＳ＠ - -<I>Il' 1l-/3 1l-/3 1l-/3 - ah}* at Here {(v/3)} represents the large-scale averaged velocity for the {3-phase, {Ed* represents the largescale volume fraction for the {3-phase and K$ represents the large-scale permeability tensor for the {3phase. We have considered only the case of the flow of two immiscible fluids, thus the large-scale equations for the 'Y-phase are identical in form to those for the {3-phase. The terms in the momentum equation involving a{ E/3}*lat and v iJ{ EIl}*/iJt result from the transient nature of the closure problem, while the terms containing vV{(p/3)Il}/3. v<I>/3 and <1>/3 are the results of nonlinear variations in the largescale field. All of the latter three terms are associated with second derivatives of the pressure and thus present certain unresolved mathematical problems. The situation concerning the large-scale capillary pressure is equally complex, and we indicate the functional dependence of {pJc by { Pc } c -_ CP' ( { Yfr )/3}1l ' a{ E/3}* ) Ell }* ' ( Py _ Pil ) g, "{( v P/3 at ' etc .. Because of the highly nonlinear nature of the capillary pressure-saturation relation, small causes can have significant effects, and the treatment of the large-scale capillary pressure is a matter of considerable *Permanent address: Department of Chemical Engineering, University of California, Davis, CA 95616 U.S.A. 342 MICHEL QUINTARD AND STEPHEN WHITAKER concern. On the basis of the derived closure problems, estimates of ull' UIl , etc., are available and they clearly indicate that the nontraditional terms in the momentum equation can be discarded when IH <;;.':t. Here IH is the characteristic length scale for the heterogeneities and .':t is the characteristic length scale for the large-scale averaged quantities. When IH is not small relative to .':t, the nontraditional terms must be considered and nonperiodic boundary conditions must be developed for the closure problem. Detailed numerical studies presented in Part II (Quintard and Whitaker, 1990) and carefully documented experimental studies described in Part III (Berlin et aI., 1990) provide further insight into the effects of large spatial and temporal gradients. Keywords. Two-phase, heterogeneous media, large-scale averaging, dynamic effective properties. Nomenclature Roman Letters A{3W scalar that maps a{E{3}*lat onto p{3w A{3T) scalar that maps a{E{3}*lat onto p{3T) AWT) interfacial area between the w-region and the 1]-region contained within "11'"" m AW<T interfacial area between the w-region and the O'-region contained within "11'"" m AT)<T 2 interfacial area between the 1]-region and the O'-region contained within "11'"" m a{3w a{3T) b{3w b{3T) B{3w B{3T) c{3w c{3T) C{3w C{3T) D{3w D{3T) D{3w D{3T) E{3w E{3T) E{3w E{3T) fjiw fjiT) g 'Je 2 2 vector that maps (a{ E{3}*lat) onto V(3w, m vector that maps (a{ E{3}*lat) onto V(3T)' m vector that maps (V{(P{3){3}{3 - P{3g) onto fi!3w, m vector that maps (V{(P{3){3}{3- P{3g) onto fi!3T)' m second-order tensor that maps (V{(p{3){3}{3 - P{3g) onto v{3w, m2 second-order tensor that maps (V {(p{3 )(3}{3 - P{3g) onto v(3T)' m2 vector that maps (Va{E{3}*lat) onto fi!3w, m vector that maps (Va{E{3}*lat) onto P{3T)' m second-order tensor that maps (Va{E{3}*lat) onto v{3w, m2 second-order tensor that maps (Va{E{3}*lat) onto v{3T)' m2 third-order tensor that maps (V!l{3) onto v{3w, m third-order tensor that maps (V!l{3) onto v{3T)' m second-order tensor that maps (V!l{3) onto P{3w, m2 second-order tensor that maps (V!l{3) onto P{3T)' m2 third-order tensor that maps (V<I» onto v{3w, m third-order tensor that maps (V<I» onto v(3T)' m second-order tensor that maps (VC() onto P{3w second-order tensor that maps (VC() onto p{3T) Pc = fjiw( E{3w) , capillary pressure relationship in the w-region Pc = fjiT)( E{3T)' capillary pressure relationship in the 1]-region gravitational vector, mls 2 largest of either 'Jew or 'JeT) TWO-PHASE FLOW IN HETEROGENEOUS POROUS MEDIA I Kf3 Kf3w Kf37J {Kf3}f3 Kf3 343 unit base vector in the x-direction unit tensor local volume-averaged .a-phase permeability, m2 local volume-averaged .a-phase permeability in the w-region, m 2 local volume-averaged .a-phase permeability in the 1]-region, m 2 large-scale intrinsic phase average permeability for the .a-phase, m2 Kf3 - {Kf3 }f3, large-scale spatial deviation for the .a-phase permeability, m2 Kf3w - {Kf3 }f3, large-scale spatial deviation for the .a-phase permeability in the w-region, m2 Kf37J - {K f3 }f3, large-scale spatial deviation for the .a-phase permeability in the 1]-region, m2 large-scale permeability for the .a-phase, m 2 characteristic length for the .a-phase at the pore-scale, m characteristic length for the u-phase at the pore-scale, m characteristic length associated with local volume-averaged quantities, m characteristic length associated with large-scale averaged quantities, m i = 1, 2, 3, lattice vectors for a unit cell, m characteristic length for the w-region, m characteristic length for the 1]-region, m characteristic length associated with a local heterogeneity, m {(Kf3 + il{Kf3}f3) . VDf3} + {Kf3}f3 . {VDf3}, m3 unit normal vector pointing from the w-region toward the 1]-region (OW7J = - n 7Jw ) unit normal vector pointing from the w-region toward the u-region (owu = - ouw) pressure in the .a-phase, N/m 2 local volume-averaged intrinsic phase average pressure in the .a-phase, N/m 2 large-scale intrinsic phase average pressure in the capillary region of the .a-phase, N/m 2 local volume-averaged intrinsic phase average pressure for the .a-phase in the w-region, N/m 2 local volume-averaged intrinsic phase average pressure for the .a-phase in the 1]-region, N/m 2 (pf3)f3 - {(Pf3)f3}f3, large scale spatial deviation for the .a-phase pressure, N/m 2 344 Pc {PcY TO Ro r CJJt{3 S{3 ｓｾ＠ fI' V{3 fl'x fI'{3 , fI'-y fl'c Vw V7) V w(13) V 7)(,a) V,,(,a) w MICHEL QUINTARD AND STEPHEN WHITAKER (p{3)'!, - {(P{3){3}{3, large scale spatial deviation for the ,a-phase pressure in the w-region, N/m 2 (p ＨＳＩｾ＠ - {(p {3 ){3}{3, large scale spatial deviation for the ,a-phase pressure in the 7]-region, N/m 2 (p-y)-Y - (P{3){3, capillary pressure, N/m 2 large-scale capillary pressure, N/m 2 radius of the local averaging volume, m radius of the large-scale averaging volume '1/'-x, m position vector, m {(I<{3 + il{K{3}(3) . VE{3} + {K{3}{3· {VE{3}' m E{3IE, local volume-averaged saturation for the ,a-phase {E{3}*/{ E}*, large-scale average saturation for the ,a-phase time, s {(I<{3 + il{K{3}(3) . VA{3} + {K{3}{3 . {VA{3}, m {(I<{3 + il{K{3}(3) . VC{3} + {K{3}{3 . {VC{3}, m 2 ,a-phase velocity vector, mls local volume-averaged phase average velocity for the ,a-phase in the wregion, mls local volume-averaged phase average velocity for the ,a-phase in the 7]region, mls large-scale intrinsic phase average velocity for the ,a-phase in the capillary region of the ,B-phase, m/s large-scale phase average velocity for the ,a-phase in the capillary region of the ,a-phase, m/s (v (3) - {( v (3)}{3, large-scale spatial deviation for the ,a-phase velocity, m/s (v (3) w - {( v(3)}f3, large-scale spatial deviation for the ,a-phase velocity in the w-region, m/s (v{3)7) - {(v{3)}{3, large-scale spatial deviation for the ,a-phase velocity in the 7]-region, m/s local averaging volume, m 3 volume of the ,a-phase in fI', m 3 large-scale averaging volume, m 3 capillary region for the ,a-phase within fl'x, m 3 capillary region for the y-phase within fl'x, m 3 intersection of fI'{3 and fI'-y, m 3 volume of the w-region within fl'x, m 3 volume of the 7]-region within fl'x, m 3 capillary region for the ,a-phase within the w-region, m 3 capillary region for the ,a-phase within the 7]-region, m 3 fl'x - Vw(,a) - V7)(,B), region in which the ,a-phase is trapped at the irreducible saturation, m 3 velocity of the interface between the active ( capillary) region and the inactive (trapped fluid) region, mls TWO-PHASE FLOW IN HETEROGENEOUS POROUS MEDIA I y 345 position vector relative to the centroid of the large-scale averaging volume, m Greek Letters local volume-averaged porosity Ef3 local volume-averaged volume fraction for the f3-phase Ef3w local volume-averaged volume fraction for the f3-phase in the w-region Ef3l] local volume-averaged volume fraction for the f3-phase in the 1]-region Ef3u local volume-averaged volume fraction for the f3-phase in the o--region (This is directly related to the irreducible saturation.) large-scale intrinsic phase average volume fraction for the f3-phase large-scale phase average volume fraction for the f3-phase large-scale spatial average volume fraction for the f3-phase Ef3 - {Ef3}f3, large-scale spatial deviation for the f3-phase volume fraction Ef3w - {Ef3}f3, large-scale spatial deviation for the f3-phase volume fraction in the w-region Ef3l] - {Ef3}f3, large-scale spatial deviation for the f3-phase volume fraction in the 1]-region a generic local volume-averaged quantity associated with the f3-phase mass density of the f3-phase, kg/m3 mass density of the y-phase, kg/m 3 P'Y viscosity of the f3-phase, N s/m 2 I-LI3 viscosity of the y-phase, N s/m 2 interfacial tension of the f3 - y phase system, N/m {(I<f3 + ｾｻｋｦＳｽＩ＠ . ｾｮｦＳｽ＠ + {Kf3}f3 . ｻｾｮｦＳｽＧ＠ N/m /1'13/ '11x, volume fraction of the f3-phase capillary (active) region /1'-/ '1Ix, volume fraction of the y-phase capillary (active) region /l'j/l'x, volume fraction of the w-region (¢w + ¢l] = 1) /l'l]/'1Ix, volume fraction of the 1]-region (¢w + ¢l] = 1) V'{(Pf3)f3}f3 - Pf3g, N/m 3 V'{(p,,)Y}'Y - P'Yg, N/m 3 E 1. Introduction The structure of the problem under consideration has been presented in detail by Quintard and Whitaker (1987, 1988), and here we note only that we are interested in solving two-phase flow problems in heterogeneous porous media such as that illustrated in Figure 1. The local volume average permeability tensors are identified as Kf3w and Kf37) where the subscript f3 refers to the fluid identified as the f3-phase and the subscripts wand 1] refer to the wand 1]-regions, respectively. These permeability tensors are to be used with the local volume averaged momentum equations (Whitaker, 1986b) which we list as 346 MICHEL QUINTARD AND STEPHEN WHITAKER Fig. 1. A stratified porous medium. (V{3)W = - K{3w. (V(p{3)r;, - P{3g) , i-t{3 (V{3)Yj = - K{3Yj. i-t{3 ＨｖｰｻＳＩｾ＠ - P{3g) , in V w (f3), (1.1) in VYj(f3). (1.2) If one is concerned with the problem of two-phase flow in a precisely defined, stratified system, one can solve the boundary-value problem associated with Equations (1.1) and (1.2) provided one is willing to assume that the individual strata are homogeneous. This is the type of analysis presented by Yokoyama and Lake (1981) among many others; however, in real life the strata can rarely be thought of as homogeneous (Giordano et al., 1985). In our original study of this problem (Quintard and Whitaker, 1988), we proposed the following definition (by exclusion) of a heterogeneous porous medium: A porous medium is homogeneous with respect to a given process and a given averaging volume when the effective transport coefficients in the volume-averaged transport equations are independent of position. If a porous medium is not homogeneous, it is heterogeneous. The dependence on the size of the averaging volume should be clear from Figure 1. At the level of local volume averaging the medium is clearly heterogeneous while TWO-PHASE FLOW IN HETEROGENEOUS POROUS MEDIA I 347 it will be homogeneous at the level of large-scale averaging if the stratified system is spatially periodic. The idea that the coefficients in a spatially smoothed transport equation depend on the size of the averaging volume brings us face-to-face with the relativistic point of view put forth by Baveye and Sposito (1984), Cushman (1984) and Whitaker (1986c). The latter study indicates how several different measuring devices impact upon the local volume averaged dispersion equation and the associated closure problem. In Part II of this paper (Quintard and Whitaker, 1990), we present results which indicate that the large-scale permeability tensor, Kill, depends on the size of the averaging volume. This situation invites for the use of weighting functions (Matheron, 1965; MarIe, 1967; Anderson and Jackson, 1967) and an implementation of the philosophy put forth by Baveye and Sposito (1984) and by Cushman (1984, 1987). However, the complexities associated with the problem under investigation suggest that there is much to be learned by following the traditional route (Whitaker, 1967; Slattery, 1967) and avoiding the use of weighting functions. While our numerical and laboratory experiments (Quintard and Whitaker, 1990; Bertin et al., 1990) deal with stratified systems, we are interested in more general systems such as the one illustrated in Figure 2. Furthermore, the theory is not limited to discrete variations in the physical properties of the medium under consideration but also applies to continuous variations. This aspect of the theory has been documented by Quintard and Whitaker (1987) for the case of single phase flow in heterogeneous porous media; however, we find it convenient to Fig. 2. Heterogeneous porous medium ｴｗｍ･ｾ＠ MICHEL QUINTARD ･ｸ｡ｭｾｮ＠ AND STEPHEN W HITAKER . Figures 1 and 2. de III Id become 'onII-phase mo l illustrated Ie)that coueegion 348 . in tenns of the (fo; wc<c expms oue 'deas keep in mind that t ee 2 if the satumtwn n effective pomus It is important to ion iIIustmted m F>gu ation. This leads to;_ hase lakes place the iITedue,blc Ihe flow of the . in which the tmpped in the e W Ie" than m equa shown in Figme h e identified the eeg'o an impenneable dium of the I Figure 3, we be considered aS t be the case . . e ,t ean m in the ,,-eeg,on ;hen ,,-eegion, sme only he II-phase. T h." nee no s omus med,. d phase is tmpped of view of the flow of :tme of a helemgen:o: identified the solid from the pom this means that the stm uee 4 wheee w." :olume as Vu(m Fi e fo; the y-phase a;;:iS idea is iIIustmted in the analys>s is fluid specific. ed phases wdhm I is suggestive of tha'gS6b). In thatstudy, olume of the trapp clature used here dia (Whitaker, 1 hile the rzgld vf t Th nomen ous me hases, w ew in poc ed as the and Y-P thin k of V.( a and V.( y). fl 'ds were IdentIfi. 4 one should . I'lar meamng wo-phase flo I F'gme , ) has a SIm . o . ·ble m the two the if-phas<:. n h e while VA y . lied the actm solid was identIfied ats to the flow of the (3-hP a{3s_p' hase can flow IS in Figure . wa. h '''pee . wh·eh e . volume rigid sobd The eegion m , t -scale avemgmg for the. fo I'-Phase{3' and for the large (1.3) r the -p hase domaIn d as, 4 it can b e ex Presse ｾ［･ｴ＠ ｩｭＢｾＧ＠ ｴｹｾ･＠ a: ｳ［ｴｵｾ･＠ ｾｶ＠ Ｚｾ＠ ｉＢＬｾ･Ｍｳ｡ｬ＠ ｨｯｭｧ･ｮｾｳ＠ ｾ＠ ｾ｟＠ ｡ｶ･ＬＢｧｭｾｳ､＠ ｾ＠ m｡ｾ＠ ｣ｳｾｯｷｮ＠ u f/ = f/= - V ({3). hase is given by • domain for the l'-P . '1 ar If3y, the actlVe SImI f/=f/=-Vu(l') l' (1.4) TWO-PHASE FLOW IN HETEROGENEOUS POROUS MEDIA I 349 .. ':.: :. '://..:..:..>/.;.:/:.: ... . , .... :., .... :'. " , . . .' ' . ' '. ',,:.: :.:: .. :: '. . . . '" .. ,".: .. ", ,'" Fig. 4. Trapped fluid in an heterogeneous porous medium. and the intersection of 'Vf3 and 'Vy is identified in Figure 4 as the capillary region. It is this region which is used in the definition of the large-scale capillary pressure, {pJc and we identify this region as 'Ve . Darcy's law, as expressed by Equations (1.1) and (1.2), makes use of two types of averages. The velocity, (v13), is a phase average while the pressure, (Pf3)f3, is an intrinsic phase average. The situation for two-phase flow in heterogeneous porous media is more complex because of the existence of inactive regions for both the {3 and y-phases, and several different average values need to be used in the theoretical analysis. For some generic function '-f;f3 associated with the {3-phase, we list three of these averages as {'-f;f3}f3 J ｾ＠ = '-f;f3 d V, 'Vf3 ·lIiJ {'-f;f3} = Ｂｴｾ＠ '-f;f3 dV, % {'-f;f3}* = f. 'f3 ｾ＠ x J.'x '-f;f3dV, large-scale intrinsic phase average (1.5) large-scale phase average (1.6) large-scale spatial average. (1.7) Similar averages exist for a function '-f;y associated with the y-phase. In addition, the large-scale capillary pressure is defined in terms of the volume of the capillary region according to 350 MICHEL QUINT ARD AND STEPHEN WHITAKER (1.8) Here Pc is the local capillary pressure defined by (1.9) in which (py)Y and (p{3){3 represent intrinsic phase average pressures for the yand {3-phases respectively. 2. Theory The analysis of two-phase flow in heterogeneous porous media can be carried out in terms of only the {3-phase since the equations for the y-phase will follow by analogy. For the two-region model illustrated in Figure 2, the boundary-value problem can be presented in terms of the local volume-averaged equations according to (2.1) (2.2) (2.3) (2.4) (2.5) (2.6) (2.7) (2.8) Here we have used a subscript to indicate in which region a function is evaluated, i.e., E{3w represents the {3-phase volume fraction in the w-region and ＨｰｻＳＩｾ＠ represents the {3-phase pressure in the 7]-region. The convention that is used for the outwardly directed unit normal vector is such that Dl)<T is a unit vector directed from the 7]-region toward the O'-region. We have used w to represent the velocity of the boundary between the two regions, thus if w . Dl)<T is positive the active 7]region is advancing into the inactive O'-region. This represents a process of imbibition for the {3-phase. The volume of the active region for the {3-phase is given by (2.9) while the volume of the inactive region was defined earlier by Equation (1.3). • TWO-PHASE FLOW IN HETEROGENEOUS POROUS MEDIA I 351 We refer to Equations (2.1) through (2.8) as local volume averaged equations and boundary conditions, since they were obtained from the appropriate point equations and boundary conditions by local volume averaging (Whitaker, 1986b). The process of local volume averaging is shown schematically in Figure 5 where the first change of scale is illustrated by 1j1---7 (1jI). It is the second change of scale that interests us in this work, and this is suggested in Figure 5 by (1jI) ---7 {( 1jI)}. Sometimes, there is a tendency to think of the large-scale process as being dominated by the structure of the heterogeneities and that the small-scale phenomena are lost in the averaging procedure. In reality, the mathematics of the method of volume averaging clearly indicates that physical information at one scale is carried forward to the next scale but filtered in the process. The exact nature of this filtering process is not entirely clear at this time, but it will be in the future loeal- sea Ie averaging '(IJ') large - scale averaging - - {(IJ')} Fig. 5. Changes of scale in the analysis of two-phase flow in heterogeneous porous medium. 352 MICHEL QUINTARD AND STEPHEN WHITAKER as more detailed numerical experiments are carried out and more well-planned laboratory studied are completed. Our starting point for this particular analysis is identical to that for our previous study, and in this new effort we have two principal goals: (1) To develop a solution to the transient closure problem. (2) To develop a theory that is valid for large spatial gradients. In order to accomplish this second goal we need to avoid, as much as possible, the imposition of length-scale constraints. Our previous analysis made repeated use of the length-scale constraints given by (2.10) Here lH represents the length scale associated with the heterogeneities (either lw or IT/ in Figure 2), Ro represents the radius of the averaging volume shown in Figure 4, and :£ represents the length scale associated with an average quantity such as {tfif3}f3. In this development we want to avoid, as much as possible, any length scale constraints, and we want to identify where the necessary length scale constraints appear in the analysis. The large-scale form of the continuity equation is obtained by integrating Equations (2.1) and (2.8) over Vw(P) and VT/(f3), respectively, to obtain (2.11) _1 J 'V x ｖｾＨｦＳＩ＠ J aEf3T/ dV + _1 at /Ix ｖｾＨｦＳＩ＠ V'. (vf31T/ dV = O. (2.12) One then makes use of the general transport theorem (Whitaker, 1981, Sec. 3.4) and the spatial averaging theorem (Howes and Whitaker, 1985) to obtain the following form * a{ Ef3}* + V' . {(v )} = O. at f3 (2.13) The averages here are defined by Equations (1.6) and (1.7), but to be explicit we express the large-scale spatial average volume fraction as ｻｅｦＳｽＪ］ｾ＠ /Ix Vw (f3) ｅｦＳｷ､ｶＫｾ＠ /Ix ｖｾＨｦＳＩ＠ ｅｦＳｔＯ､ｶＫｾ＠ /1= V rr(f3) Ef3a dV- (2.14) It is of some importance to note that Equation (2.13) is restricted by the fact that * In our previous work (Quintard and Whitaker, 1988, Equation (2.21» we used ｅｾ＠ to represent the large-scale spatial average of Ef3. This was an unfortunate oversight, since ｅｾ＠ does not conform to the nomenclature indicated by Equation (1.7). 353 TWO-PHASE FLOW IN HETEROGENEOUS POROUS MEDIA I the {3 and y-phases are immiscible so that there is no possibility of diffusional transport of the {3-phase into the y-phase. That case has been considered elsewhere (Whitaker, 1977). The question to be considered at this time is: what are the length scale constraints associated with Equation (2.13)? To begin with, we require that fix be as large or larger than the local averaging volume, two of which are illustrated in Figure 1. From a practical point of view, that imposes no length scale constraint at all. Veverka (1981) and, more recently, Mis (1987) have pointed out that the derivative of the volume average must exist in order to use the spatial averaging theorem to derive Equation (2.13) from Equations (2.11) and (2.12). Mis concludes that the derivative of the average exists almost everywhere and provides an example of a case for which the derivative does not exist. Mis also points out that the derivative of the average (in this case {(v{3)}) exists everywhere if the radius of curvature of the interface (in this case the surface between the several distinct regions) is always less than the radius of curvature of the averaging volume. Howes and Whitaker (1985) also found conditions for which the derivative of the volume average is not defined on a set of measure zero; however, that does not prevent the use of the spatial averaging theorem to derive Equation (2.13), nor does it prevent the use of Equation (2.13) to determine scalar and vector fields associated with {E{3}* and {(v{3/}. In summary, we think that are no length-scale constraints of importance associated with Equation (2.13). 2.1. THE MOMENTUM EQUATION The treatment of the momentum equation given here will differ from our previous analysis since we wish to avoid imposing length-scale constraints. This will require that we be very careful about identifying where an average quantity is evaluated. To this end, we will use the system illustrated in Figure 6 where the position vector r, is expressed as r = x + y. (2.15) Here the vector x locates the centroid of "fix and y locates any point relative to the centroid. We begin our analysis with Equation (2.2) and integrate over V w({3) to obtain ｾｊ＠ "fix V w ({3) ＨｖｻＳＯｷ､ｶ］Ｍｾｬ＠ J.L{3 "fI= J V w ({3) K{3w.C'V(p{3/e-p{3g) dV l· (2.16) The left-hand side of this result poses no problems at all since it will be added to the comparable integral of (v{3/Tj to produce the desired large-scale velocity, {(v{3/}. The single problem in Equation (2.16) rests with the average of the product on the right-hand side, and we attack this problem by means of the decompositions (2.17) 354 MICHEL QUINTARD AND STEPHEN WHITAKER r x Fig. 6. Location of points within the averaging volume. The decompositions for the 7)-region take the form (2.18) and this means that the average of Kf3 will be zero only when {Kf3}f3 is a constant, i.e., 1f =- 'Vf3 v w(f3) 1f Kf3w dV+Kf37]dV=O. 'Vf3 v ry(f3) A A (2.19) In all of our previous studies of multiphase transport phenomena, we have always assumed that averaged quantities could be treated as constant within the averaging volume and that the average of the spatial deviations was zero. We now wish to avoid these assumptions which were based on a length-scale constraint of the form (Carbonell and Whitaker, 1984; Quintard and Whitaker 1987). (2.20) This result is based on a Taylor series expansion of averaged quantities and is thought to be overly severe. Numerical results presented in Part II of this paper tend to confirm this. Keeping in mind that Equation (2.20) is not a part of this analysis, we substitute Equation (2.17) into Equation (2.16) and express the result as ｾｦ＠ 'V x Vw (f3) (Vf3>w dV TWO-PHASE FLOW IN HETEROGENEOUS POROUS MEDIA I 355 The integration in Equation (2.21) is performed with x held constant, thus any quantity evaluated at x can be removed from the integral. In order to simplify the nomenclature for the steps that follow, we use (2.22) and express Equation (2.21) as (2.23) Here we have used the subscript y in order to indicate that a function is evaluated at the position given by x + y. We can now r use of the definitions {Kf3 }f3 ly = {Kf3}f3lx + ｾｻｋ､ｦＳｉｙＬ＠ !lf3ly = !lf3lx + ｾＡｬｦＳｉｙ＠ (2.24) (2.25) in order to arrange Equation (2.23) in the Here one must be very careful to note that we have adopted the convention that averaged quantities outside an integral sign are evaluated at r = x, whereas all terms inside an integral sign are evaluated at r = x + y. If we set ｾｻｋｦＳｽ＠ and ｾＡｬｦＳ＠ equal to zero in Equation (2.26) we recover the original result given by Quintard and Whitaker (1988, Eq. 2.26). By leaving Equation (2.26) in its present form we avoid imposing any length-scale constraint such as the one represented by Equation (2.20). The penalty that we pay fOI avoiding length-scale constraints is the fact that Equation (2.26) is an integral-differential equation. Refering to Equation (2.22), we see that V'{(Pf3)f3}f3 appears under an integral sign in Equation (2.26) and this will eventually lead to a very complicated computational problem. For 7]-region, we can develop an equation analogous to Equation (2.26) and 356 MICHEL QUINTARD AND STEPHEN WHITAKER when the two are added we obtain = - :{3 [ ( ｾ＠ {K{3}{3 + {K(3 + L1{K{3}{3}) . !l{3 + + {(K{3 + L1{K{3}(3) . V'h} + {K{3}{3· {V'h} J - :{3 [{(K{3 + L1{K{3}(3)· L1!l{3} + {K{3}{3 . {L1!l{3} J (2.27) If we set MK{3}{3 and L1!l{3 equal to zero and note that this is consistent with setting {K{3} equal to zero, we recover Equation (2.30) of Quintard and Whitaker (1988). One should think of the flow of the ,8-phase as being caused by the pressure gradient, !l{3 = V'{(p{3){3}{3 - P{3g, and this suggests that the last term.in Equation (2.27) should be thought of as an integral source of momentum that generates a flow in the ,8-phase. We can simplify the representation of Equation (2.27) if we identify this integral source as (2.28) in order to obtain + {(K{3 + MK{3}(3) . V'h} + {K{3}{3· {V'h}J 1 - - ｾｻＳＮ＠ (2.29) /-L{3 Here we have identified the volume fraction of the capillary region for the ,8phase as (2.30) In order to develop a closed form of Equation (2.29), we need to obtain a representation for the spatial deviation of the pressure, P{3' that was defined earlier by Equations (2.17) and (2.18). This is the objective of the next section. 3. Closure We begin our development of the closure problem with an analysis of the ,8-phase in the w-region, and in that region the spatial deviations for the volume fraction, TWO-PHASE FLOW IN HETEROGENEOUS POROUS MEDIA I 357 the velocity, and the pressure are given by (3.1) (V{3)W = {(v{3)}{3 + v{3w, ＨｰｻＳＩｾ＠ = {(p{3){3}{3 (3.2) + P!3w. (3.3) A comparable set of definitions is available for the 1)-region; however, it is sufficient to derive the w-region equations and then accept the 1)-region equations by analogy. 3.1. CONTINUITY EQUATION Substitution of Equations (3.1) and (3.2) into Equation (2.1) immediately leads to aE{3w+V. V w=_[a{Ef3}{3 +V.{(v )}{31 at {3 at {3 J (3.4) and one should think of the term on the right-hand side of this result as a source for the E{3w-field. To be useful, this source must be represented in terms of the large-scale variables, {E{3}* and {(v{3)}. With the aid of Equations (1.5) through (1. 7), we can express these variables as {E{3}* = <P{3{ E{3}f3 + ｾ＠ f; = J V,,({3) Ef3(T dV, {(Vf3)} = <Pf3{(vf3)}f3. (3.5) (3.6) Use of Equations (3.5) and (3.6) in Equation (3.4) eventually leads to aEf3w r7 A --+v·V at f3w (3.7) provide that we take into account that Ef3u is independent of time. This result corresponds to Equation (3.6) of Quintard and Whitaker (1988) except for the fact that our previous representation failed to include <P73 1 as a multiplier of the last two integrals. The correspondence with our previous work can be easily seen if one takes into account that (3.8) 358 MICHEL QUINT ARD AND STEPHEN WHITAKER At this point, it should be clear that the source on the right-hand side of Equation (3.7) is intimately related to the behavior of the inactive regions. One must remember that we have characterized the inactive regions as regions where the velocity is zero. In addition, our entire analysis is restricted to immiscible fluids so that there is no diffusional transport between the two phases. This means that E{3<r cannot depend on time, and if we further characterize the inactive regions as having a value of E{3u that is independent of position we can express Equation (3.7) as aE{3w - + '1'7v • Vｾ＠ {3w -- <P{3-2 [({}* E{3 - E{3<T ) -a<p{3 + {< v (3 )} • '1'7 v <P{3 J. at at (3.9) While it is unlikely that the value of E{3(T will be the same in both the wand 1]regions, this form of Equation (3.7) clearly indicates that the source in the transport equation for E{3w is dominated by the temporal and spatial gradients of <P{3' That spatial and temporal gradients of <P{3 exist for typical imbibition and drainage problems is well known; however, we are going to delay our investigation of this phenomena until some later date. We believe that this aspect of two-phase flow in heterogeneous porous media is extremely important; however, the theoretical analysis is a bit beyond our reach at this time and the remainder of our development will be restricted by a<p{3 = 0 at (3.10) ' From a practical point of view, these conditions are usually obtained by requiring that <P{3 = 1, and this is what we have done in both our numerical and laboratory experiments which are described in Parts II and III of this paper. 01 Ie basis of Equations (3.10) we write Equation (3.9) as aE{3w '1'7v 'v{3w= ｾ＠ 0 --+ . (3.11) at We are now ready to develop the closure equation for v{3w' 3.2. MOMENTUM EQUATION Our analysis of the momentum equation begins with Equation (2.2) and makes use of the decompositions given by Equations (2.17) and (3.2). This leads to (3.12) As was the case with the continuity equation, we seek a representation in terms of the large-scale variables, thus Equation (3.6) is used in order to incorporate the large-scale phase average velocity into the closure equation for v{3w' In addition, we TWO-PHASE FLOW IN HETEROGENEOUS POROUS MEDIA I 359 again use the simplified nomenclature !lf3 = V{(Pf3)f3}f3 - Pf3g (3.13) so that Equation (3.12) can be expressed as Vf3w = - <Pi3 1 {(Vf3)} - ｾ＠ /J-f3 [{K f3 }f3·!lf3 + Kf3w·!lf3 + Kf3w· vhwl· (3.14) The large-scale averaged momentum equation, given by Equation (2.29) can now be substituted into Equation (3.14) to obtain the w-region velocity deviation equation. vf3w 1 A A = - - [(Kf3w - {Kf3 + il{Kf3}f3}f3) .!lf3 + Kf3w· Vhw/J-f3 - {(Kf3 + il{K f3 }f3)· Vh}f3 - {K f3 }f3· {Vh} f3 l + + -1 <Pf3-1 q,f3. (3.15) /J-f3 At this point, we see that there are two sources in the equation for vf3w; the first of these is the large-scale pressure gradient represented by !lf3 and given explicitly by Equation (3.13) and the second is the integral source, $f3' which is defined by Equation (2.28). The 1)-region continuity and momentum equations can be written by analogy to Equations (3.11) and (3.15), but before doing so we need to consider the boundary conditions for the spatial deviations. 3.3. BOUNDARY CONDITIONS Use of Equation (3.2), along with a comparable decomposition for the 1)-region, in the boundary condition given by Equation (2.3) leads to (3.16) Following a similar procedure with the boundary condition given by Equation (2.4) yields (3.17) On the basis of the restriction given by Equation (3.10) we can set w • OW<T and W· oT)<T equal to zero in Equations (2.5) and (2.6). Use of the appropriate decompositions for the velocities then leads to B.C.3 0W<T· vf3w l.C.4 0WT· vf3T) = {(vf3)}f3, at A w,,, (3.18) 0T)<T· {(vf3)}f3, at AT)(T. (3.19) = - 0W(T· - Use of Equations (2.29) and (3.15), along with a bit of algebra leads to (3.20) 360 MICHEL QUINTARD AND STEPHEN WHITAKER By analogy the boundary condition at the 1] - (J surface takes the form (3.21 ) At this point, we have available to us the spatial deviation form of the continuity and momentum equations, along with an appropriate set of boundary conditions for vJ3w, V13'7' PJ3w, and PJ3'7. In our previous analysis of this problem, we invoked a quasi-steady constraint which allowed us to set the time derivative of EJ3w equal to zero. This provided an attractive simplification of Equation (3.11) along with the comparable equation for the 1]-region. In this treatment, we wish to retain the general form of Equation (3.11) and develop a transient solution to the closure problem. This will be done within the framework of local mechanical equilibrium (Quintard and Whitaker, 1988, Sec. 2.3) which requires that the local volume fractions, EJ3w and EJ3'7' are uniquely determined by the capillary pressure-saturation curves for the w-region and the 1]-region. This requires that the capillary number be small compared to one, and this is a situation that occurs for most cases of two-phase flow in porous media. Our previous treatment of the largescale capillary pressure was based on the simplification that averaged quantities could be treated as constant within the averaging volume, and we now wish to avoid that simplification. This requires that we must reconsider the representation for the large-scale capillary pressure. 3.4. CAPILLARY PRESSURE The local capillary pressure is defined by Equation (1.9) and when we use decompositions of the form given by Equation (3.3), the local capillary pressure takes the form (3.22) This relation is only valid in the capillary region illustrated in Figure 4 and we have made use of the nomenclature suggested by both Figures 4 and 5. In order to develop the connection between {PeY and the pressure gradients in Equation (2.29) and the comparable equation for {(vy)}, we make use of the representation {(py)Y}Yly = {(py)Y}Ylx+ ｦｾＺＫｹ＠ Oy·dr+pyg·y (3.23) along with a comparable expression for the {3-phase pressure. This allows us to express Equation (3.22) as Pely = ({(Py)Y}Y - {(PJ3)J3}J3)x + AYJ3ly + (py - h)y + (py - PJ3)g· y, (3.24) where the scalar AYJ3ly is defined by AYJ3ly = ｦｾＺＫｙ＠ (Oy - Ofj) . dr. (3.25) 361 TWO-PHASE FLOW IN HETEROGENEOUS POROUS MEDIA I We remind the reader once again that 01' and 0f3 are defined as (3.26) Given that the large-scale capillary pressure is defined by (3.27) we can use Equation (3.24) to obtain {PcY Ix = ({(PI') 1'} l' - {(Pf3 )f3}f3)x + {AI'f3 K + + (PI'- Pf3)g· ｻｹｽｾ＠ + {PI'- hK. (3.28) When this result is substituted into Equation (3.24), we obtain the local capillary pressure in terms of the large-scale capillary pressure and the dynamic or flow variables, AI'f3' etc. Pcly = {PcYlx + AI'f3ly - ｻａｉＧｦＳｽｾ＠ + + (PI' - h)y - {PI' - + (PI' - Pf3)g· (y - ｻｹｽｾＩ＠ ｐｴｊｾﾷ＠ (3.29) One must remember that when y locates a point in the w-region, Pc is a function of Ef3w and we express this idea as Pcly = @Pw(Ef3w) = {pX + AI'f3 - {AI'f3}' + + (PI' - Pf3)g· (y - {y}c) + (PI' - h) - {PI' - hY· (3.30) Here it is understood that all averaged quantities are evaluated at r = x while all other quantities are evaluated at r = x + y. The function @Pw( Ef3w) will depend on the process under consideration; however, we will assume that the function is known and that the inverse can be found. This allows us to express the local volume fraction as (3.31) and at this point we can begin to see that the dynamic closure problem is very complex. Since the function, @Pw , is often highly nonlinear, the inverse will also be highly nonlinear and small causes, such as the influence of 0f3 or Pf3' may well give rise to large effects. The spatial decomposition given by Equation (3.1) allows us to write (3.32) and with the help of Equation (3.5), we can express this result in terms of the large-scale spatial average volume fraction ｅｦＳｷｬｹ］Ｍ｣＾ｩＱｻｽＪｉｙＫｾ＠ f/{3 v if(f3) Ef3adV+Ef3wly· (3.33) 362 MICHEL QUINTARD AND STEPHEN WHITAKER Here one must keep in mind that ¢>{3 is neither a function of space nor time on the basis of the severe restriction imposed by Equation (3.10). Later it will become very clear that all averaged quantities that are source terms in the closure problem need to be evaluated at r = x, thus we use a Taylor series expansion to express Equation (3.33) as (3.34) One must be careful to remember that in this particular equation the position vector y locates only points in the w-region as suggested by Equation (3.31). We can now substitute Equation (3.31) into Equation (3.34) in order to express the time rate of change of E{3w as (3.35) A quick review of Equations (3.30) and (3.31) will indicate that the time derivative of ｾ［＠ I evaluated at r = x + y represents a very complex function. Along these lines, it is important to note that the dynamic effects associated with fl{3' Pf3' etc, represent nontrivial contributions in some of our calculations that are presented in Part II of this paper. These contributions result from the coupling between the highly nonlinear capillary pressure-saturation relation given by Equation (3.31) and the highly nonlinear permeability-saturation relation given by (3.36) While the dynamic effects in Equation (3.31) can become important in the determination of Kf3w, we believe that the first two terms on the right-hand side of Equation (3.35) plus the dependence of ｾ［ｬ＠ on {PcY' represent the dominant effects in the determination of aE{3jat. This must be considered as an hypothesis at this time and to make matters perfectly clear we state it as such. Hypothesis I. In developing the transient form of the closure continuity equation, it is sufficient to express the local volume fraction only in terms of the largescale capillary pressure. To be explicit, this means that Equation (3.31) is used in the form (3.37) This implies that Ef3w is constant throughout the w-region and we know that this is not true for some cases of practical interest. It is essential to keep in mind that Equation (3.37) is not used in the resolution of Equation (3.36) but only as an approximation in Equation (3.35) where its use leads to TWO-PHASE FLOW IN HETEROGENEOUS POROUS MEDIA I 363 (3.38) At this point, it seems necessary to comment on the manner in which the {Pc}'" - {E{3}* relation is constructed. One begins by assigning a value to {Pc y, n{3' (Py - P{3)g, etc., along with estimates of P{3 and Pr This determines the local capillary pressure and with a knowledge of f!J' wand f!J' 7J one can determine E{3w and E{37J within an averaging volume located at r = x + y. Knowledge of E(3w and E{37J allows one to determine K{3w and K{37J and thus solve for the closure variables according to the method described in Section 4. When this is done P{3 and py can be determined and the above process repeated in order to produce a large-scale, dynamic capillary pressure as a function of {E{3}* and the other parameters that appear in Equation (3.30). This procedure means that the derivative of {Pc}' with respect to {E{3}* that appears in Equation (3.38) actually contains some of the ､ｹｮ｡ｾＭ ic effects that influence E(3w. However, the only manner in which these dynamic effects can be correctly computed is in terms of the time derivative of f!J';;;1 that appears in Equation (3.35). Obviously, the validity of Hypothesis I remains to be explored. In order to simplify Equation (3.38), we define a scalar function, 'Jew, according to 'Je = (af!J';;;l ) (a{pJC) _ cfJ-1 w a{pJc x a{ E{3}* x {3 (3.39) so, that Equation (3.38) takes the form ='Je wa{E{3}*/_cfJjily.v a{E{3}*I. M{3w/ at at y x at (3.40) x We are now ready to return to the closure continuity equation given by Equation (3.11), and with the help of Equation (3.40) we obtain 'Je a{E{3}*/_cfJ- 1 y.v a{E{3}*/ +V'v =0 w at x {3 at x {3w· (3.41 ) Here it is important to think of a{ E{3}* at and V a{ E{3}* at as sources for the v{3w-field and it is important to note that we have arranged the analysis so that these sources are evaluated ,._ i' = x and thus can be considered 364 MICHEL QUINT ARD AND STEPHEN WHITAKER as constants with regard to the solution of the closure problem. So far, we have not done this with the sources, fif3 and «)f3' that appear in the momentum equation given by Equation (3.15) or with the sources that were presented earlier as Equations (3.20) and (3.21). To accomplish this we use the representations fi f3 ly = fif3lx + Y • V'fi f3 lx, (3.42) «)f3ly = cIlf3lx + y. V'«)f3lx (3.43) and then summarize the closure problem for the two-region model as follows: (3.44) Vf3w 1 = - - A A [(Kf3W - {Kf3 + Ll{Kf3}f3}(3)· fif3 + Kf3W· V'hw- fLf3 - {(l<f3 + Ll{Kf3}(3) . V'h}f3 - {Kf3}f3 . {V'h}f3 + + (Kf3W - {Kf3 + Ll{Kf3}f3}(3)y: V'fi f3 ] + 1 -1.... 1 -1 rI +-cPf3 '¥f3+-cPf3 Y' V«)f3' fLf3 fLf3 (3.45) (3.46) (3.47) (3.48) (3.49) Vf3Y/ = - - 1 fL{3 A A [(Kf3Y/ - {K(3 + a{K{3}{3}(3) . fi{3 - {(Kf3 + Ll{K{3}(3) . V'h}{3 - {K{3}{3· {V'h}{3 + + (K{3y/ - {Kf3 + Ll{K{3}{3}(3)y: V'fi{3] + 1 -1 rI + -1 cP{3_/.... '¥f3 + - cP{3 y. V«)f3 fL{3 + K{3y/ . V'hY/- (3.50) fL{3 'Je a{E{3}* - cP- l y. V' a{E{3}* + V' . V = y/ at {3 at {3Y/ o. (3.51) The convention employed here is that every large-scale averaged quantity is evaluated at r = x except for Ll{K{3}{3 which appears under an integral sign and is therefore evaluated at r = x + y. TWO-PHASE FLOW IN HETEROGENEOUS POROUS MEDIA I 365 At this point, we can see that there are six source terms in the closure problem, and we list them as When the closure problem is treated as quasi-steady and the length-scale constraint given by Equation (2.20) is imposed, this development reduces to the one given earlier by Quintard and Whitaker (1988, Sec. 3). In that case, there is only the single source term, 0/3' Removing the length-scale constraint given by Equation (2.20) and including a portion of the transient closure process (remember Equations (3.10» has led us to an extremely complicated problem. It seems likely that some parts of the closure problem may be more important than others; however, extensive calculations will be required in order to learn if simplifications can be made in Equations (3.44) through (3.51). In the next section we present the closure problem in terms of the five required closure variables. This leads to five boundary-value problems of the type illustrated by Equations (3.44) through (3.51). While this appears to present an overwhelmingly difficult computational task, all five problems can be solved using a single computer routine when the la:rge-scale volume averaged flow is one-dimensional. 4. The Closure Variables The boundary-value problem for the large-scale spatial deviations of the volume fraction, the velocity, and the pressure contains six source terms. On the basis of previous studies of closure problems (Ryan et at., 1981; Nozad et al., 1985; Ochoa et al., 1986; Whitaker, 1986a; Quintard and Whitaker, 1987; Quintard and Whitaker, 1988), we propose the following representations in which all averaged quantities are evaluated at r = x. (4.1) (4.3) 366 MICHEL QUINT ARD AND STEPHEN WHITAKER Here we have used boldface italic sans serif letters to represent the third-order tensors, such as Ef3w and Df3Tj. In our previous study of two phase flow in heterogeneous porous media, we were able to prove (Quintard and Whitaker, 1988, Appendix B) that our representations for vf3 and Pf3 were unique if, and only if, Kf3 were spatially periodic, symmetric, and a single principle axis coordinate system existed for the capillary region, r f3. It seems unlikely that these conditions can be satisfied for the problem currently under investigation, thus Equations (4.1) through (4.4) must be thought of as plausible but not yet proved. From Equation (2.29) we can see that we need only determine 'VPf3 in order to complete the closure and determine the coefficients in the large-scale momentum equation. This means that there are five variables b f3 , A f3 , cf3' Df3 and Ef3 , that need to be determined and we refer to these as the closure variables. We begin the process of developing the boundary-value problems for these variables by substituting Equations (4.1) through (4.4) into Equations (3.45) and (3.50). Restricting our presentation to the w-region, we find Bf3w = -[Kf3w· 'Vbf3w - {K f3 }f3· {'Vbf3}f3 - {(l<f3 + Ll{K f3 }f3) . 'Vbf3}f3 + + (Kf3W - {Kf3 + Ll{Kf3}f3}f3)], af3w = - Kf3w . 'VAf3w + {Kf3}f3 . {'VAf3}f3 + {( Kf3 + Ll{K f3 }f3) . 'VA f3 }f3, Cf3w = -Kf3w· 'Vcf3w + {K f3 }f3· {'VCf3}f3 + {(Kf3 + Ll{K f3 }f3)· 'VCf3}f3, (4.5) (4.6) (4.7) Df3w = - [Kf3w· 'VDf3w - {K f3 }f3· {'VDf3}f3 - {(Kf3 + Ll{Kf3 }f3) . 'VDf3}f3 + + (Kf3W - {Kf3 + .Ll{Kf3}f3}f3)y], (4.8) Ef3w = -[Kf3w· 'VEf3w - {K f3 }f3· {'VEf3}f3 - {(Kf3 + Ll{Kf3 }f3) . 'VEf3}f3 + <pi3 1 (ly)] (4.9) Substitution of Equation (4.1) into Equation (3.44) leads to the following set of continuity equations: Yew + V' . af3w = 0, (4.10) - <Pi3 1 y + 'V. Cf3w = 0, :4.11) 'V. Bf3w = 0, (4.12) 'V. Df3w = 0, (4.13) V' . Ef3w = 0. (4.14) Here one must remember that all the averaged terms in the representation for vf3w are evaluated at r = x and are therefore constant with respect to the different- 367 TWO-PHASE FLOW IN HETEROGENEOUS POROUS MEDIA I iation that occurs in Equation (3.44). The continuity equations for the 1]-region follow by analogy, but we will not list them here. In order to obtain the boundary-value problems for the closure variables, we substitute Equations (4.5) through (4.9) into Equations (4.10) through (4.14) and we carry out the analogous procedure for the 1]-region. In addition, we make use of Equations (4.1) through (4.4) in the boundary conditions given by Equation (3.46) through (3.49) in order to develop the boundary conditions for the closure variables. This procedure leads to the following boundary value problems. PROBLEM I (bf3 and the effect of !lf3) V· [Kf3w' Vbf3W - {K f3 }f3· {Vb{3}f3 - {(I<{3 + Il{K f3 }f3) . Vb{3}f3 + + (I<{3w - {Kf3 + Il{Kf3}f3}f3)] = 0, B.C.1. (4.15a) "wTj' Kf3w . Vb f3w = llwTj . Kf3Tj . Vbf3Tj + llwTj . (Kf3Tj - Kf3w ), at A wTj , (4.15b) (4.15c) (4.15d) B.C.4. -llTja-' K{3Tj' Vbf3Tj = "Tja-' Kf3Tj' (4.15e) at ATja-, V· [Kf3Tj' Vbf3Tj - {K{3}f3 . {Vb{3}{3 - {(K{3 + Il{Kf3}(3) . Vb{3}{3 + + (Kf3Tj - {Kf3 + Il{K{3}f3}f3)] = 0. ( 4.15f) PROBLEM II ( Af3 and the effect of a{ ｾｽＪＩ＠ V· [Kf3w . VAf3W - {Kf3}f3 . {VA{3}f3 - {(K{3 + Il{Kf3}(3) . VAf3}f3] = 'Jew, at A wTj , (4.16a) B.C.1. llwTj·K f3w ·VA{3w=llwTj·K{3Tj·VA{3Tj' B.C.2. Af3w=Af3Tj, B.C.3. -llwa-' Kf3w' VAf3w = 0, at Awa-, (4.16d) B.C.4. -llTja-' K{3Tj' VAf3Tj at ATja-, (4.16e) (4.16b) (4.16c) at A wTj , = 0, (4.16f) PROBLEM III (Cf3 and the effect of V a{ ｾｽＪＩ＠ V· [Kf3w . VCf3w - {Kf3}f3 . {VC{3}f3 - {(Kf3 + Il{K{3}f3) . VC{3}f3] B.C.1. "wTj' Kf3w' VC{3w = "wTj' Kf3Tj' VCf3Tj' at A wTj , = - cPi3! y, (4.17a) (4.17b) (4.17c) (4.17d) 368 MICHEL QUINT ARD AND STEPHEN WHITAKER B.C.4. -nT)a' Kf3T)' VCf3T) = 0, at AT)a, (4.17e) V· [Kf3T)' VCf3T) - {Kf3}f3 . {VCf3}f3 - {(I<f3 + MKf3}f3) . VCf3}f3] = - cPi3'y. (4.17f) PROBLEM IV (Df3 and the effect of V!lf3) V· [Kf3w . VDf3w - {Kf3}f3 . {VD f3 }f3 - {(I<f3 + ｾｻｋｦＳｽＩ＠ + (Kf3W - {Kf3 + ｾｻｋｦＳｽＩｹ｝＠ = . VDf3}f3 + (4.18a) 0, B.C.l. nWT)' Kf3w' VDf3w = nWT)' Kf3T)' VDf3T) + nWT)' (Kf3T) - Kf3w)y, B.C.2. Df3w = Df3T)' at AWT)' (4.18b) (4.18c) at AWT)' (4.18d) (4.18e) V· [Kf3T)' VDf3T) - {Kf3}f3 . {VD f3}f3 - {(Kf3 + ｾｻｋｦＳｽＩ＠ + (Kf3T) - {Kf3 + ｾｻｋｦＳｽＩｹ｝＠ = . VDf3}f3 + (4.18f) 0. PROBLEM V (Ef3 and the effect of V<I>f3) V· [Kf3w' VEf3w - {Kf3}f3 . {VE f3 }f3 - {(I<f3 + ｾｻｋｦＳｽＩ＠ . VEf3}f3] = cPi3 11 , (4.19a) (4.19b) (4.19c) (4.19d) B.C.4. -nT)a' Kf3T)' VEf3T) = 0, at AT)if' (4.1ge) V . [Kf3T)' VEf3T) - {Kf3}f3 . {VE f3 }f3 - {(I<f3 + ｾｻｋｦＳｽＩ＠ . VEf3}f3] = cPi3'I. (4.19f) While the form of the boundary-value problem for all five closure variables is identical, thus providing a certain degree of simplification in terms of computer programming, there is an enormous computational effort involved in the solution of Equations (4.15) through (4.19). This results largely from the fact that the coefficients will depend on position through the functional dependence of Kf3w and Kf3T) on Ef3w and Ef3T)' From our discussion in Section 3.4, we know that Ef3w and Ef3T) will be determined by capillary pressure relations such as that given by Equation (3.30). In our previous discussion, we left ftf3 and fty unspecified; however, we are now in a position to use Equations (4.2) and (4.4) to write Pc = g; w( Ef3w) = {pJc + (py - Pf3)g . (y - {y n+ +!ly' [(y + b yw ) - {y + byr] -!lf3' [(y + b f3w ) - {y + bf3}'] + + iV!ly: [(yy + Dyw) - {yy + Dy}C] - iV!lf3: [(yy + Df3w ) - {yy + Df3}C] + TWO-PHASE FLOW IN HETEROGENEOUS POROUS MEDIA I 369 (4.20) In this representation of the local capillary pressure, it is understood that all averaged functions are evaluated at r = x and are therefore constant in terms of the solution of the closure problem. It is possible that computational exploration will reveal some simplification in this representation, and we hope that the simplification indicated by Hypothesis I and Equation (3.37) will survive this exploration. In their present form, the boundary-value problems given by Equations (4.15) through (4.19) are not complete since no boundary conditions have been specified at the entrances and exits of the region under consideration. The traditional approach to this problem is to use a spatially periodic model of the system under consideration (Brenner, 1980) and then impose a periodic condition on the closure function. In Figure 7 we have shown a model of a spatially periodic porous Fig. 7. Spatially periodic. heterogeneous porous medium. 370 MICHEL QUINT ARD AND STEPHEN WHITAKER I w -region Fig. 8. Unit cell for a spatially periodic. heterogeneous porous medium. medium, and a unit cell from that model is illustrated in Figure 8. For the boundary-value problem given by Equations (4.15) we would express the periodicity condition as (4.21) and this would allow us to determine bf3 to within an arbitrary constant. The arbitrary constant can be removed by imposing the following condition on the average value of bf3 (4.22) We shall see in Section 6 that the average values of the closure variables play no role in the determination of the coefficients that appear in the large-scale momentum equation, thus Equation (4.22) is unnecessary. However, it is consistent with the traditional treatment which considers the average of the spatial deviations to be zero, i.e. (4.23) The key point to be made here is that Equations (4.21) and (4.22) are consistent with boundary value problems in which the following conditions are satisfied: TWO-PHASE FLOW IN HETEROGENEOUS POROUS MEDIA I 371 (1) The transport coefficients in the governing differential equations are spatially periodic (2) The average value of the volumetric source terms that appear in the governing differential equations is zero, (3) The average value of the surface source terms that appear in the boundary conditions is zero. In general, these conditions are consistent with Equation (2.20) and that is a length-scale constraint that we do not wish to impose since it would reduce this analysis to that presented previously by Quintard and Whitaker (1988). At this time, we know of no general treatment of nonperiodic boundary conditions; however, we have a partial solution for one-dimensional processes and this is presented in the next section. s. Nonperiodic Boundary Conditions In this development we want to avoid the restrictions imposed by Equations (4.23), and to do so we need to consider the average of the deviation more carefuly. By definition we have {Vf3}=+JV (f3) ｖｦＳｷ､Ｋｾｊ＠ "f; x 'Vx w Vry(f3) vf3l)dV (5.1) and use of the decompositions (5.2) leads us to a general expression for the average of the deviation. {vf3}f3lx = {(Vf3)}Ix - ｾ＠ 'Vx I {(vf3)}ly dV. (5.3) 11'(3 If {(vf3)} is assumed to be constant within the averaging volume, Equation (5.3) reduces to the first of Equations (4.23). If not, it is convenient to use the type of nomenclature illustrated by Equations (2.24) and (2.25) in order to express Equation (5.3) as ｻｖｦＳｽ］Ｍｾｉ＠ 'rx A{(vf3)}dV (5.4) Y{3 If one can compute the large-scale averaged velocity field, then the integral on the right-hand side of Equation (5.4) can be evaluated and the average of the deviation determined. In general, it is best to carry out this calculation in terms of Equation (5.4); however, it is of some value to use a Taylor series expansion to express this result as (5.5) 372 MICHEL QUINTARD AND STEPHEN WHITAKER If "fI'/3 is equal to "fI'x, or if the capillary region is uniformly distributed within "fI'x, the value of {y} is zero, and it is the second derivative of the average value of the velocity that controls the average of the deviation (Whitaker, 1988). In order to calculate {v/3}, we prefer Equation (5.4) to Equation (5.5) since it requires the integration of numerical data rather than differentiation. Given that {v/3} can be evaluated, we can use the representations given by Equations (4.1) and (4.3) to arrive at the following relation V/3 ｽＭｾｻｂ＠ - £lo { /3 } 'U/3+ a/3 }a{€/3}* + {C} /3' at A { JL/3 v a{€/3}* + at 1 1 1 + - {D/3}: Vfi/3 + - {E/3}: ｖｾＯＳ＠ + - <1>/3' (5.6) JL/3 JL/3 JL/3 Here we are confronted with the problem of how to partition the average of v/3 among the five closure variables. Since we believe that fi/3 is the dominant source for the spatial deviations we resolve this problem with the following hypothesis. HYPOTHESIS II. The dominant driving force for the spatial deviations is fi/3' thus we associate the non-zero value of {v /3} entirely with B/3' This means that we replace Equation (5.6) with the following constraints on the averages {B/3} . fi/3 {a/3} = 0, = JL/3{ v /3} - ｾＯＳＧ＠ (5.7) {C/3} = 0, {D/3} = 0, {E/3} = 0, (5.8) In the following paragraphs we will consider only Equation (5.7) in detail since the results for Equations (5.8) can easily be deduced by analogy. We can make use of the representation for B/3w given by Equation (4.5), along with the analogous representation for B/31]' in order to express Equation (5.7) as -[{V· Ｈｋｾ｢ＯＳＩｽ＠ - {(V· - {{(K/3 + ｾｻｋＯＳｽＩ＠ ｋｾＩ｢ＯＳｽ＠ - {{K/3}/3· {Vb/3}/3}/3- . Vb/3}/3}/3 + {K/3 - {K/3 + ｾｻｋＯＳｽ｝＠ . fi/3 (5.9) In going from Equation (5.7) to Equation (5.9), we have made use of the general relation (5.10) and in preparation for an iterative determination of the nonperiodicity of b/3 we express Equation (5.9) as {V· Ｈｋｾ｢ＯＳＩｽﾷ＠ fi/3 = - JL/3{v/3}/3 + ｃｐｩＳＱｾＯ＠ + [{(V . ｋｾＩ｢ＯＳｽ＠ + HK/3 + ｾｻｋＯＳｽＩ＠ - {K/3 - {K/3 + ｾｻｋＯＳｽ｝＠ + + {{K/3}/3· {Vb/3}/3}/3 + . Vb/3}/3}/3 . fi/3' (5.11) 373 TWO-PHASE FLOW IN HETEROGENEOUS POROUS MEDIA I We want to use this result with the unit cell illustrated in Figure 8, and for such a unit cell we can use Equations (4.15b) and (4.15c) to obtain (5.12) Here Awe and ATJe represent the areas of the entrances and exits for the wand Y]regions of the unit cell shown in Figure 8, and "we and "TJe represent the outwardly directed unit normal vectors for these two areas. We can substitute Equation (5.12) in Equation (5.11) and arrange the result as [ｾ＠ I {3 "we . Awe + [{(V· ｋｾＩ｢ｻＳｽ＠ ｋｾｷ｢ｻＳ＠ dA + ｾ＠ V{3 I "TJe . ａｾ･＠ ｋｾｔｊ｢ｻＳ＠ dA 1J . fi{3 + {{K{3}{3· {Vb{3}{3}{3 + + {{(K{3 + d{K{3}{3) . Vb{3}{3}{3 - {K{3 - {K{3 + d{K{3}{3}{3}{3] . fi{3. (5.13) Here one can see the possibility of determining the nonperiodicity of b{3 in some average sense, but it is not clear how Equation (5.13) could be used to obtain general, nonperiodic conditions of the form b{3w(r + t) = b{3w(r) + ､｢ｲｾｷＨＩＬ＠ i = 1, 2, 3, (5.14) b{3TJ(r + l;) = b{3TJ(r) + db{3TJ(r), i = 1,2,3. (5.15) On the other hand, if we consider a one-dimensional flow in the x-direction for the stratified system illustrated in Figure 9, we can determine the desired nonperiodicity in b{3. In order for a flow to occur in a direction perpendicular to a stratified system it is necessary that V {3 be equal to r ¢i3 1 = 1), thus we have % ( {!f;{3}{3 = 1 ｊｘｾ＠ L + <l>ryL ｸｾＭ＼ｉ＾ｷｌ＠ !f;{3 dx. (5.16) For this one-dimensional process, we use the notation (5.17) so that Equation (5.13) takes the form 374 MICHEL QUINT ARD AND STEPHEN WHITAKER y r-----L-----oOof x X:O Fig. 9. Unit cell in a stratified porous medium. Here we have made use of Equation (4.15c) in order to write (5.19) The analogous results for Equations (5.8) are easy to obtain, and here we only list the result that can be derived from the first of Equations (5.8). This is given by (5.20) TWO-PHASE FLOW IN HETEROGENEOUS POROUS MEDIA I 375 Equations (5.18), (5.20) and the others that follow from Equations (5.8) provide nonperiodic conditions for b f3 , A p , etc., that are consistent with Equations (4.15) through (4.19) for a one-dimensional process, and that are consistent with the behavior of the large-scale averaged velocity field as indicated by Equation (5.4). In order to develop this nonperiodic condition, we were forced to adopt Hypothesis II which allowed us to proceed from Equation (5.6) to Equations (5.7) and (5.8). It seems quite possible that the source term on the right hand side of Equation (5.7) does not playa dominant role in the determination of the nonperiodic conditions. This is to say that it is the variation of Kf3 with position that controls the nonperiodicity, and if this is the case our partitioning of Equation (5.6) will prove to be successfull. Implementation of the current theory for two and three-dimensional closure problems requires that either equations of the form of Equations (5.14) and (5.15) be extracted from Equation (5.13) or a demonstration that some average constraint on the nonperiodicity of the closure functions is sufficient. 6. The Large-Scale Problem Two parts of the large-scale problem have already been presented, i.e., the largescale continuity equation is given by Equation (2.13) and the large-scale capillary pressure is implied by Equations (3.28), (3.29) and (4.20). To be explicit, we list the continuity equation for the ,B-phase a{Ef3}* + V . {(vf3)} = 0 at (6.1) and we present the large-scale capillary pressure in terms of the closure variables and Equation (3.28). {pJc = ({(py)YP' - {(pp)f3}P) + (py - P(3)g· {yY + + U y . {y + byY - Up . {y + bpY + + ｾｖｕｹＺｻ＠ + Dy}C - ｾｖｕｦＳＺｻｙ＠ + {,uyA y - + Df3Y + C a{ Ep}* . a{ Ep}* ,uf3Ap} - - + {,uyc y - ,uf3Cp}c . V - - + at at (6.2) In order to develop the large-scale momentum equation for the ,B-phase, we return to Equation (2.29) and make use of the representations for Pf3 which are given by Equations (4.2) and (4.4). This leads to {(vp)} = _l.- K$' (V{(pp)f3}P ,up 1 Pf3g) - uf3 a{Ef3}* - Uf3' V a{Ef3}* at at 1 1 - - AIlf3: VV{(Pf3)P}f3 - - r!llf3: ｖｾｦＳ＠ - - ｾｦＳＧ＠ ,uf3 ,uf3,up (6.3) 376 MICHEL QUINTARD AND STEPHEN WHITAKER Here one must remember that (6.4) and one must keep in mind that .At{3 and rYt{3 represent third-order tensors. The large-scale permeability tensor is defined by K$ = {cP{3{K{3}{3 + {K{3 + MK{3}{3} + {(K{3 + ｾｻｋＳｽＨＩ＠ . V'b{3} + {K{3}{3· {V'b{3} ((l.5) and the remaining coefficients take the form U{3 = {(K{3 + ｾｻｋＳｽＨＩ＠ . V'A{3} + {K{3}{3· {V'A{3}, (6.6) U{3 = {(K{3 + ｾｻｋＳｽＨＩ＠ . V'c{3} + {K{3}{3 . {V'c{3}' (6.7) .M.{3 = {(K{3 + ｾｻｋＳｽＨＩ＠ . V'D{3} + {K{3}{3 . {V'D{3}, (6.8) 1!It{3 = {(K{3 + ｾｻｋＳｽＨＩ＠ . V'E{3} + {K{3}{3 . {V'E{3}. (6.9) The source term «I>{3' given earlier by Equation (2.28) can be written as +{K{3}{3· ｛ｾｦＮ＠ (f!{3ly-f!{3l x )dVJ. "fix (6.10) V{3 in order to indicate the integral-differential nature of Equation (6.3). The secondorder derivative of the pressure is clearly evident in Equation (6.3) and higherorder derivatives are hidden in the integral source terms, V'«I>{3 and <I»{3. It would appear that the term involving M{3: V'V'{(p{3){3}{3 must be discarded for mathematical reasons; however, we have carried it through the analysis in order to clearly indicate the nature of the problem with which we are faced. The origin of this term is in the closure problem, and it appears when we use Equation (3.43) in the closure momentum equation given by Equation (3.15) and in the closure boundary conditions given by Equations (3.20) and (3.21). Since the gradient of f!{3 appears naturally when Equation (3.15) is substituted into Equation (3.4), there seems to be no way to prevent second order derivatives of the pressure from appearing in the closure problem. From the closure problem, they work their way into the large scale momentum equation where they appear as M{3:V'V'{<J7(3){3}{3 in Equation (6.3). Use of the large-scale continuity equation immediately produces the unacceptable third order derivatives in the pressure. 6.1. ORDER OF MAGNITUDE ESTIMATES One can use Equation (6.1) in Equation (6.3) to obtain a form that is convenient for use with order of magnitude estimates. TWO-PHASE FLOW IN HETEROGENEOUS POROUS MEDIA I 1 1 1 - -.tlf3: VV {(Pf3)f3}f3 - -!Jlf3:VcJ>f3 - -cJ>f3' ｾｦＳ＠ • ｾｦＳ＠ 377 (6.11) ｾｦＳ＠ By examining the separate closure problems, one can produce the following order of magnitude estimates: uf3 = o(K;; IH), Uf3 = ｯＨｾＺ＠ ｬｾ＠ ), (6.12) Here IH represents the characteristic length scale for the heterogeneities as illustrated in Figure 2, while Kf3 and Kf3 should be thought of as given by Kf3 = O(Kf3w - Kf3 .,.,), (6.13) Kf3 = O(Kf3w - Kf3.,.,). (6.14) The term 'Je is representative of 'Jew and 'Je.,." with the former given explicitly by Equation (3.39) and the latter given by analogy. The estimate of cJ>f3 is obtained from the definition given by Equation (6.10) and the type of analysis represented by Equations (5.4) and (5.5). We choose :£ to be the characteristic length scale for large-scale averaged quantities, and this leads to (6.15) (6.16) Clearly, these two terms will be negligible compared to {(Vf3)} when IH is small compared to :£ and they are likely to be important when :£ approaches IH. The remaining three terms on the right-hand side of Equation (6.11) can be estimated as (6.17) (6.18) (6.19) 378 MICHEL QUINT ARD AND STEPHEN WHITAKER From these estimates, it seems clear that Equation (6.11) can be simplified to the traditional form K* {(VJ3» = _--.Ii. (V{(PJ3)J3}J3 - PJ3g) ILJ3 (6.20) when the following length-scale constraint is satisfied ｉｈｾ＠ :E. (6.21) It also seems clear that when :E approaches IH all the nontraditional terms in Equations (6.11) may be important. Since the last three terms in Equation (6.11) are related to second derivatives of the pressure we should expect to encounter difficulties when :E approaches IH. 7. Conclusions In this work, we have extended our previous study of two-phase flow in heterogeneous porous media. A portion of the transient closure problem has been resolved, but the influence of spatial and temporal gradients of the volume fraction of the inactive or dead regions remains as on open problem. In this study we have made every effort to avoid the imposition of length scale constraints so that a theory valid for large spatial gradients could be achieved. This has led to the need for nonperiodic boundary conditions and we have been able to resolve this difficulty for only the case of one-dimensional flows at the closure level. The avoidance of length-scale constraints has led to the appearance of second derivatives of the pressure, and if these new terms are important, the mathematical problem becomes ill-posed and revisions in the theory will be necessary. Extensive numerical computation is needed to solve the five boundary-value problems for the closure variables, and both numerical and laboratory experiments are needed to test the theory. Acknowledgement This work was completed while S. Whitaker was a Visiting Scientist at the Laboratoire Energetique et Phenomenes de Transfert - U.A. CNRS 873. Services provided by the laboratory and financial support provided by NSF Grant CBT-86 11610 and CBT 88-12870 are gratefully acknowledged. Financial support provided by Elf, and by ARTEP and PIRSEM (CNRS, EDF, AFME) are gratefully acknowledged. References Anderson, T. B. and Jackson, R., 1967, A fluid mechanical description of fluidized beds, Ind. Eng. Chern. Fundarn. 6, 527-538. ｾ＠ TWO-PHASE FLOW IN HETEROGENEOUS POROUS MEDIA I 379 Baveye, P. and Sposito, G., 1984, The operational significance of the continuum hypothesis in the theory of water movement through soils and aquifers, Water Resour. Res. 20, 521-530. Bertin, H., Quintard, M., Corpel, P. V., and Whitaker, S., 1990, Two-phase flow in heterogeneous porous media III: Laboratory experiments for flow parallel to a stratified system, Transport in Porous Media 5, 543-590. Brenner, H., 1980, Dispersion resulting from flow through spatially periodic porous media. Trans. R. Soc. (London) 297,81-133. Cushman, J. H., 1984, On unifying the concepts of scale, instrumentation and stochastics in the development of multiphase transport theory, Water Resour. Res. 20, 1668-1676. Cushman, J. H., 1987, Stochastic filtering of multi phase transport phenomena, Transport in Porous Media 2, 425-453. Giordano, R. M., Salter, S. J., and Mohanty, K. K., 1985, The effects of permeability variations on flow in porous media, SPE paper 14365 presented at the 60th Annual SPE Meeting, Las Vegas, Nevada. Howes, F. A. and Whitaker, S., 1985, The spatial averaging theorem revisited, Chem. Engng. Sci. 40, 1387-1392. Mis, J., 1987, On the existence of the derivative of the volume average, Transport in Porous Media 2, 615-621. Nozad, I., Carbonell, R. G., and Whitaker, S., 1985, Heat conduction in multiphase systems I: Theory and experiment for two-phase systems, Chem. Engng. Sci. 40, 843-855. Ochoa, J. A., Stroeve, P., and Whitaker, S., 1986, Diffusion and reaction in cellular media, Chem. Engng. Sci. 41, 2999-3013. Quintard, M. and Whitaker, S., 1987, Ecoulements monophasiques en milieu poreux: effet des heterogeneites locales, 1. Meca. Thea. et Appl. 6, 691-726. Quintard, M., and Whitaker, S., 1988, Two-phase flow in heterogeneous porous media: the method of large-scale averaging, Transport in Porous Media 3, 357-413. Quintard, M. and Whitaker, S., 1990, Two-phase flow in heterogeneous porous media II: Numerical experiments for flow perpendicular to a stratified system, Transport in Porous Media 5, 429-472. Ryan, D., Carbonell, R. G., and Whitaker, S., 1981, A Theory of Diffusion and Reaction in Porous Media, AIChE Symposium Series, edited by P. Stroeve and W. J. Ward, Vo!. 77, pp. 46-62. Slattery, J. c., 1967, Flow of viscoelastic fluids through porous media, AIChE 1. 13, 1066-1071. Veverka, V., 1981, Theorem for the local volume average of a gradient revisited, Chem. Engng. Sci. 36, 833-838. Whitaker, S., 1967, Diffusion and dispersion in porous media, AIChE 1. 13(3),420-427. Whitaker, S., 1977, Simultaneous heat, mass and momentum transfer in porous media: a theory of drying, Advances in Heat Transfer, Vol. 13. Academic Press, New York, pp. 119-203. Whitaker, S., 1981, Introduction to Fluid Mechanics, R. E. Kreiger Pub!. Co., Malabar, Fla. Whitaker, S., 1986a, Flow in porous media I: A theoretical derivation of Darcy's law, Transport in Porous Media 1, 3-25. Whitaker, S., 1986b, Flow in porous media II: The governing equations for immiscible, two phase flow, Transport in Porous Media 1, 105-125. Whitaker, S., 1986c, Multiphase transport phenomena: Matching theory and experiment, in G. Papanicolaou (ed.), Advances in Multiphase Flow and Related Problems, SIAM, Philadelphia, pp. 273-295. Whitaker, S., 1988, Comments and corrections concerning the volume-averaged temperature and its spatial derivation, Chem. Engng. Commun. 70, 15-18. Yokoyama, Y. and Lake, L. W., 1981, The effects of capillary pressure on immiscible displacements in stratified porous media, 56th Annual Fall Meeting of the SPE, San Antonio, Oct. 5-7, SPE number 10109.