1 Introduction

[page 2323, §1]   
[2323.1.1] A macroscopic theory of two phase flow inside a rigid porous medium poses not only challenges to nonequilibrium statistical physics and geometry [1], but is also crucial for many applications [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]. [2323.1.2] Despite its popularity the accepted macroscopic theory of two phase flow seems unable to reproduce the experimentally observed phenomenon of saturation overshoot [11].

[2323.2.1] Models for twophase flow in porous media can be divided into macroscopic (laboratory or field scale) models popular in engineering, and microscopic (pore scale) models such as network models [12, 13, 14, 15, 16, 17] that are popular in physics. [2323.2.2] As of today no rigorous connection exists between microscopic and macroscopic models [1, 18]. [2323.2.3] In view of the predominantly non-specialist readership with a physics background it is appropriate to remind the reader of the traditional theory, introduced between 1907 and 1941 by Buckingham, Richards, Muskat, Meres, Wyckoff, Botset, Leverett and others [19, 20, 21, 22, 23]^a (This is a footnote:) ^aThe following introductory paragraphs are quoted from Ref.[24] for convenience of the interdisciplinary readership and following an explicit request from the editor.. [2323.2.4] One formulation of the traditional macroscopic theory starts from the fundamental balance laws of continuum mechanics for two fluids (called water W and oil O) inside the pore space (called P) of a porous sample S=P∪M with a rigid [page 2324, §0]    solid matrix (called M). [2324.0.1] Recall the law of mass balance in differential form

where ϱi⁢x,t,ϕi⁢x,t,vi⁢x,t denote mass density, volume fraction and velocity of phase i=W,O as functions of position x∈S⊂R3 and time t∈R+. [2324.0.2] Exchange of mass between the two phases is described by mass transfer rates Mi giving the amount of mass by which phase i changes per unit time and volume. [2324.0.3] Momentum balance for the two fluids requires in addition

where Σi is the stress tensor in the ith phase, Fi is the body force per unit volume acting on the ith phase, mi is the momentum transfer into phase i from all the other phases, and Di/D⁢t=∂/∂⁡t+vi⋅∇ denotes the material derivative for phase i=W,O.

[2324.1.1] Defining the saturations Si⁢x,t as the volume fraction of pore space P filled with phase i one has the relation ϕi=ϕ⁢Si where ϕ is the porosity of the sample. [2324.1.2] Expressing volume conservation ϕW+ϕO=ϕ in terms of saturations yields

[2324.1.3] In order to obtain the traditional theory these balance laws for mass, momentum and volume have to be combined with specific constitutive assumptions for Mi,mi,Fi and Σi.

[2324.2.1] Great simplification is afforded by assuming that the porous medium is rigid and macroscopically homogeneous

although this is often violated in applications [25]. [2324.2.2] Let us focus first on the momentum balance (2). [2324.2.3] One assumes that the stress tensor of the fluids is diagonal

where PW,PO are the fluid pressures. [2324.2.4] Realistic subsurface flows have low Reynolds numbers so that the inertial term

can be neglected in the momentum balance equation (2). [2324.2.5] It is further assumed that the body forces

are given by gravity. [2324.2.6] As long as there are no chemical reactions between the fluids the mass transfer rates vanish, so that MW=-MO=0 holds. [2324.2.7] Momentum transfer between the fluids and the rigid matrix is dominated by viscous drag in the form

[page 2325, §0]    where μW,μO are the constant fluid viscosities, k is the absolute permeability tensor, and kWr⁢SW,kOr⁢SW are the nonlinear relative permeabilitiy functions for water and oil ^b (This is a footnote:) ^bkWr⁢SW,kOr⁢SW account for the fact, that the experimentally observed permability of two immiscible fluids deviates from their partial (or mean field) permeabilities k⁢SW,k⁢SO obtained from volume averaging of the absolute permeability. .

[2325.1.1] Inserting the constitutive assumptions (4)–(8) into the mass balance eq. (1) yields

while the momentum balance eq. (2)

give the generalized Darcy laws for the Darcy velocities ϕi⁢vi [3, p. 155]. [2325.1.2] Equations (9) and (10) together with eq.  (3) provide 9 equations for 12 primary unknowns SW,SO,ϱW,ϱO,PW⁢PO,vW,vO. Additional equations are needed.

[2325.2.1] Observations of capillary rise in regular packings [26] suggest that the pressure difference between oil and water should in general depend only on saturation [23]

where σW⁢O is the oil-water interfacial tension and κ⁢SW is the mean curvature of the oil-water interface. [2325.2.2] The system of equations is closed with two equations of state relating the phase pressures and densities. [2325.2.3] In petroleum engineering the two fluids are usually assumed to be incompressible

while in hydrology one thinks of water W and air O setting

where Rs≈287J kg-1K-1 is the specific gas constant and the temperature T is assumed to be constant throughout S.

[2325.3.1] When the fluids (water and oil) are incompressible (as in petroleum engineering) eqs. (12a) and (12b) hold. In this case, adding equations (9a) and (9b), using eq.  (3) and integrating the result shows

where the total volume flux Q is independent of x and qi=ϕi⁢vi=ϕ⁢Si⁢vi with i=W,O are the volume flux of water and oil. [2325.3.2] Inserting eqs. (10) into eq.  (13) and using eq.  (11) to eliminate PW gives

[page 2326, §0]    where (with i=W,O)

are the mobilities λi, total mobility λ and fractional flow functions fi, respectively. [2326.0.1] Multiplying eq.  (10a) with fO, eq.  (10b) with fW and subtracting eq.  (10b) from eq.  (10a), using eq.  (13) and fW+fO=1 to eliminate qO gives the result

which can be inserted into eq.  (9a) to give

a nonlinear partial differential equation for the saturation field SW⁢x,t. [2326.0.2] For small k or when gravity and capillarity effects can be neglected the last two terms vanish and eq.  (17) reduces to the Buckley-Leverett equation [27]

a quasilinear hyperbolic partial differential equation. [2326.0.3] Equation (17) supplemented with a (quasilinear elliptic) equation obtained from ∇⋅Q=0 by defining a global pressure in such a way that the total flux Q obeys a Darcy law with respect to the global pressure provides, for incompressible fluids, an equivalent formulation of eqs. (9)-(12b). [2326.0.4] For compressible fluids the situation is different.

[2326.1.1] When W corresponds to water and O to air (as for applications in hydrology) eqs. (12c) and (12d) hold. [2326.1.2] The large density difference ϱO≪ϱW suggests to consider the case ϱO≈0 of a very rarified gas or vacuum as a first approximation. [2326.1.3] For ϱO=0 eq.  (9b) is identically fulfilled, eq.  (12d) implies PO=0 and then eq.  (10b) implies vO=0. [2326.1.4] In this way the O-phase vanishes from the problem and one is left only with the W-phase. [2326.1.5] Inserting eq.  (10a) into eq.  (9a) and using eq.  (11) gives the Richards equation [20]

for saturation or

for pressure after writing SW⁢x,t=Pc-1⁢-PW⁢x,t=θ⁢PW⁢x,t with the help of eq. (11). [2326.1.6] To define the nonlinear function θ⁢PW the typical sigmoidal shape has been assumed for Pc⁢x.

[2326.2.1] The quasilinear elliptic-parabolic Richards equation (19) is the basic equation in hydrology, while the quasilinear hyperbolic Buckley-Leverett equation (18) is fundamental for applications in petroleum engineering. [2326.2.2] Both equations, (18) and (19), differ from the general fractional flow formulation (17) in terms of saturation SW and global pressure P. [2326.2.3] They differ also from the formulation in terms of SW,SO,vW,vO,PW,PO given by (3),(9),(10) and (11). [2326.2.4] These latter equations appropriately supplemented with initial and boundary conditions and spaces of functions resp. generalized functions for the unknowns constitute the traditional theory of macroscopic capillarity in porous media.

[page 2327, §1]    [2327.1.1] The question of domains is important for wellposedness and numerical solution. [2327.1.2] For eq.  (18) it is well known that classical solutions, i.e. locally Lipschitz continuous functions, will in general exist only for a finite length of time [28, 29, 30]. [2327.1.3] Hence it is necessary to consider also weak solutions [31]. [2327.1.4] Weak solutions are locally bounded, measurable functions satisfying eq.  (18) in the sense of distributions. [2327.1.5] Weak solutions are frequently constructed by the method of vanishing viscosity or the theory of contraction semigroups. [2327.1.6] For the Richards equation (19) a domain of definition in the space of Bochner-square-integrable Sobolev-space-valued functions has been discussed in [32]. [2327.1.7] In many engineering applications formulations such as eqs. (17), (18) or (19) with (11) augmented with appropriate initial and boundary conditions are solved by computer programs [33, 34, 35]. [2327.1.8] This concludes our brief introduction into the traditional theory.