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 $\mathbb{W}$ and oil $\mathbb{O}$ ) inside the pore space (called $\mathbb{P}$ ) of a porous sample $\mathbb{S}=\mathbb{P}\cup\mathbb{M}$ with a rigid [page 2324, §0] solid matrix (called $\mathbb{M}$ ). [2324.0.1] Recall the law of mass balance in differential form

$\frac{\partial(\phi _{i}\varrho _{i})}{\partial t}+\mathbf{\nabla}\cdot(\phi _{i}\varrho _{i}\mathbf{v}_{i})=M_{i}$

(1)

where $\varrho _{i}(\mathbf{x},t),\phi _{i}(\mathbf{x},t),\mathbf{v}_{i}(\mathbf{x},t)$ denote mass density, volume fraction and velocity of phase $i=\mathbb{W},\mathbb{O}$ as functions of position $\mathbf{x}\in\mathbb{S}\subset\mathbb{R}^{3}$ and time $t\in\mathbb{R}_{+}$ . [2324.0.2] Exchange of mass between the two phases is described by mass transfer rates $M_{i}$ 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

$\phi _{i}\varrho _{i}\frac{{\rm D}^{i}}{{\rm D}t}\mathbf{v}_{i}-\phi _{i}\mathbf{\nabla}\cdot\Sigma _{i}-\phi _{i}\mathbf{F}_{i}=\mathbf{m}_{i}-\mathbf{v}_{i}M_{i}$

(2)

where $\Sigma _{i}$ is the stress tensor in the $i$ th phase, $\mathbf{F}_{i}$ is the body force per unit volume acting on the $i$ th phase, $\mathbf{m}_{i}$ is the momentum transfer into phase $i$ from all the other phases, and ${\rm D}^{i}/{\rm D}t=\partial/\partial t+\mathbf{v}_{i}\cdot\mathbf{\nabla}$ denotes the material derivative for phase $i=\mathbb{W},\mathbb{O}$ .

[2324.1.1] Defining the saturations $S_{i}(\mathbf{x},t)$ as the volume fraction of pore space $\mathbb{P}$ filled with phase $i$ one has the relation $\phi _{i}=\phi S_{i}$ where $\phi$ is the porosity of the sample. [2324.1.2] Expressing volume conservation $\phi _{{\mathbb{W}}}+\phi _{{\mathbb{O}}}=\phi$ in terms of saturations yields

${S_{{\mathbb{W}}}}+{S_{{\mathbb{O}}}}=1.$

(3)

[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 $M_{i},\mathbf{m}_{i},\mathbf{F}_{i}$ and $\Sigma _{i}$ .

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

$\phi(\mathbf{x},t)=\phi={\rm const}$

(4)

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

$\displaystyle\Sigma _{{\mathbb{W}}}$	$\displaystyle=$	$\displaystyle-{P_{{\mathbb{W}}}}\mathbf{1}$		(5a)
$\displaystyle\Sigma _{{\mathbb{O}}}$	$\displaystyle=$	$\displaystyle-{P_{{\mathbb{O}}}}\mathbf{1}$		(5b)

where ${P_{{\mathbb{W}}}},{P_{{\mathbb{O}}}}$ are the fluid pressures. [2324.2.4] Realistic subsurface flows have low Reynolds numbers so that the inertial term

$\frac{{\rm D}^{i}}{{\rm D}t}\mathbf{v}_{i}=0$

(6)

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

$\displaystyle\mathbf{F}_{{\mathbb{W}}}$	$\displaystyle=$	$\displaystyle-\varrho _{{\mathbb{W}}}g\mathbf{e}_{z}$		(7a)
$\displaystyle\mathbf{F}_{{\mathbb{O}}}$	$\displaystyle=$	$\displaystyle-\varrho _{{\mathbb{O}}}g\mathbf{e}_{z}$		(7b)

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 $M_{\mathbb{W}}=-M_{\mathbb{O}}=0$ holds. [2324.2.7] Momentum transfer between the fluids and the rigid matrix is dominated by viscous drag in the form

$\displaystyle\mathbf{m}_{\mathbb{W}}$	$\displaystyle=$	$\displaystyle-k^{{-1}}\frac{\mu _{{\mathbb{W}}}\,\phi _{{\mathbb{W}}}^{2}}{{k^{r}_{{\mathbb{W}}}}({S_{{\mathbb{W}}}})}\mathbf{v}_{{\mathbb{W}}}$		(8a)
$\displaystyle\mathbf{m}_{\mathbb{O}}$	$\displaystyle=$	$\displaystyle-k^{{-1}}\frac{\mu _{{\mathbb{O}}}\,\phi _{{\mathbb{O}}}^{2}}{{k^{r}_{{\mathbb{O}}}}({S_{{\mathbb{W}}}})}\mathbf{v}_{{\mathbb{O}}}$		(8b)

[page 2325, §0] where $\mu _{{\mathbb{W}}},\mu _{{\mathbb{O}}}$ are the constant fluid viscosities, $k$ is the absolute permeability tensor, and ${k^{r}_{{\mathbb{W}}}}({S_{{\mathbb{W}}}}),{k^{r}_{{\mathbb{O}}}}({S_{{\mathbb{W}}}})$ are the nonlinear relative permeabilitiy functions for water and oil ^b (This is a footnote:) ^b ${k^{r}_{{\mathbb{W}}}}({S_{{\mathbb{W}}}}),{k^{r}_{{\mathbb{O}}}}({S_{{\mathbb{W}}}})$ account for the fact, that the experimentally observed permability of two immiscible fluids deviates from their partial (or mean field) permeabilities $k{S_{{\mathbb{W}}}},k{S_{{\mathbb{O}}}}$ obtained from volume averaging of the absolute permeability. .

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

$\displaystyle\frac{\partial(\varrho _{{\mathbb{W}}}{S_{{\mathbb{W}}}})}{\partial t}+\mathbf{\nabla}\cdot(\varrho _{{\mathbb{W}}}{S_{{\mathbb{W}}}}\mathbf{v}_{{\mathbb{W}}})$	$\displaystyle=$	$\displaystyle 0$		(9a)
$\displaystyle\frac{\partial(\varrho _{{\mathbb{O}}}{S_{{\mathbb{O}}}})}{\partial t}+\mathbf{\nabla}\cdot(\varrho _{{\mathbb{O}}}{S_{{\mathbb{O}}}}\mathbf{v}_{{\mathbb{O}}})$	$\displaystyle=$	$\displaystyle 0$		(9b)

while the momentum balance eq. (2)

$\displaystyle\phi _{{\mathbb{W}}}\mathbf{v}_{{\mathbb{W}}}$	$\displaystyle=$	$\displaystyle-\frac{k}{\mu _{{\mathbb{W}}}}{k^{r}_{{\mathbb{W}}}}({S_{{\mathbb{W}}}})(\mathbf{\nabla}{P_{{\mathbb{W}}}}-\varrho _{{\mathbb{W}}}g\mathbf{e}_{z})$		(10a)
$\displaystyle\phi _{{\mathbb{O}}}\mathbf{v}_{{\mathbb{O}}}$	$\displaystyle=$	$\displaystyle-\frac{k}{\mu _{{\mathbb{O}}}}{k^{r}_{{\mathbb{O}}}}({S_{{\mathbb{W}}}})(\mathbf{\nabla}{P_{{\mathbb{O}}}}-\varrho _{{\mathbb{O}}}g\mathbf{e}_{z})$		(10b)

give the generalized Darcy laws for the Darcy velocities $\phi _{i}\mathbf{v}_{i}$ [3, p. 155]. [2325.1.2] Equations (9) and (10) together with eq. (3) provide 9 equations for 12 primary unknowns ${S_{{\mathbb{W}}}},{S_{{\mathbb{O}}}},\varrho _{{\mathbb{W}}},\varrho _{{\mathbb{O}}},{P_{{\mathbb{W}}}}{P_{{\mathbb{O}}}}$ , $\mathbf{v}_{{\mathbb{W}}},\mathbf{v}_{{\mathbb{O}}}$ . 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]

${P_{{\mathbb{O}}}}-{P_{{\mathbb{W}}}}=\sigma _{{\mathbb{W}\mathbb{O}}}\kappa({S_{{\mathbb{W}}}})={P_{\mathrm{c}}}({S_{{\mathbb{W}}}})$

(11)

where $\sigma _{{\mathbb{W}\mathbb{O}}}$ is the oil-water interfacial tension and $\kappa({S_{{\mathbb{W}}}})$ 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

$\displaystyle\varrho _{{\mathbb{W}}}(\mathbf{x},t)$	$\displaystyle=\varrho _{{\mathbb{W}}}=\mathrm{const}$		(12a)
$\displaystyle\varrho _{{\mathbb{O}}}(\mathbf{x},t)$	$\displaystyle=\varrho _{{\mathbb{O}}}=\mathrm{const}$		(12b)

while in hydrology one thinks of water $\mathbb{W}$ and air $\mathbb{O}$ setting

$\displaystyle\varrho _{{\mathbb{W}}}(\mathbf{x},t)$	$\displaystyle=\varrho _{{\mathbb{W}}}=\mathrm{const}$		(12c)
$\displaystyle\varrho _{{\mathbb{O}}}(\mathbf{x},t)$	$\displaystyle=\frac{{P_{{\mathbb{O}}}}(\mathbf{x},t)}{R_{s}T}$		(12d)

where $R_{s}\approx 287$ J kg ${}^{{-1}}$ K ${}^{{-1}}$ is the specific gas constant and the temperature $T$ is assumed to be constant throughout $\mathbb{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

$\mathbf{q_{\mathbb{W}}}(\mathbf{x},t)+\mathbf{q_{\mathbb{O}}}(\mathbf{x},t)=\mathbf{Q}(t)$

(13)

where the total volume flux $\mathbf{Q}$ is independent of $\mathbf{x}$ and $\mathbf{q}_{i}=\phi _{i}\mathbf{v}_{i}=\phi S_{i}\mathbf{v}_{i}$ with $i=\mathbb{W},\mathbb{O}$ are the volume flux of water and oil. [2325.3.2] Inserting eqs. (10) into eq. (13) and using eq. (11) to eliminate ${P_{{\mathbb{W}}}}$ gives

$\mathbf{Q}=-k\lambda\left[\mathbf{\nabla}{P_{{\mathbb{O}}}}-f_{\mathbb{W}}\mathbf{\nabla}{P_{\mathrm{c}}}-\left(f_{\mathbb{W}}\varrho _{{\mathbb{W}}}+f_{\mathbb{O}}\varrho _{{\mathbb{O}}}\right)g\mathbf{e}_{z}\right]$

(14)

[page 2326, §0] where (with $i=\mathbb{W},\mathbb{O}$ )

$\lambda _{i}=\frac{k^{r}_{i}}{\mu _{i}};\qquad\lambda=\lambda _{\mathbb{W}}+\lambda _{\mathbb{O}};\qquad f_{i}=\frac{\lambda _{i}}{\lambda}$

(15)

are the mobilities $\lambda _{i}$ , total mobility $\lambda$ and fractional flow functions $f_{i}$ , respectively. [2326.0.1] Multiplying eq. (10a) with $f_{\mathbb{O}}$ , eq. (10b) with $f_{\mathbb{W}}$ and subtracting eq. (10b) from eq. (10a), using eq. (13) and $f_{\mathbb{W}}+f_{\mathbb{O}}=1$ to eliminate $\mathbf{q_{\mathbb{O}}}$ gives the result

$\mathbf{q_{\mathbb{W}}}=f_{\mathbb{W}}\left[\mathbf{Q}+k\lambda _{\mathbb{O}}(\varrho _{{\mathbb{W}}}-\varrho _{{\mathbb{O}}})g\mathbf{e}_{z}\right]+k\frac{\lambda _{\mathbb{W}}\lambda _{\mathbb{O}}}{\lambda}\mathbf{\nabla}{P_{\mathrm{c}}}$

(16)

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

$\phi\frac{\partial{S_{{\mathbb{W}}}}}{\partial t}+\mathbf{\nabla}\cdot\left\{ f_{\mathbb{W}}({S_{{\mathbb{W}}}})\left[\mathbf{Q}+k\lambda _{\mathbb{O}}({S_{{\mathbb{W}}}})(\varrho _{{\mathbb{W}}}-\varrho _{{\mathbb{O}}})g\mathbf{e}_{z}+k\lambda _{\mathbb{O}}({S_{{\mathbb{W}}}})\mathbf{\nabla}{P_{\mathrm{c}}}({S_{{\mathbb{W}}}})\right]\right\}=0$

(17)

a nonlinear partial differential equation for the saturation field ${S_{{\mathbb{W}}}}(\mathbf{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]

$\phi\frac{\partial{S_{{\mathbb{W}}}}}{\partial t}+\mathbf{Q}\cdot\mathbf{\nabla}f_{\mathbb{W}}({S_{{\mathbb{W}}}})=0$

(18)

a quasilinear hyperbolic partial differential equation. [2326.0.3] Equation (17) supplemented with a (quasilinear elliptic) equation obtained from $\mathbf{\nabla}\cdot\mathbf{Q}=0$ by defining a global pressure in such a way that the total flux $\mathbf{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 $\mathbb{W}$ corresponds to water and $\mathbb{O}$ to air (as for applications in hydrology) eqs. (12c) and (12d) hold. [2326.1.2] The large density difference $\varrho _{{\mathbb{O}}}\ll\varrho _{{\mathbb{W}}}$ suggests to consider the case $\varrho _{{\mathbb{O}}}\approx 0$ of a very rarified gas or vacuum as a first approximation. [2326.1.3] For $\varrho _{{\mathbb{O}}}=0$ eq. (9b) is identically fulfilled, eq. (12d) implies ${P_{{\mathbb{O}}}}=0$ and then eq. (10b) implies $\mathbf{v}_{{\mathbb{O}}}=0$ . [2326.1.4] In this way the $\mathbb{O}$ -phase vanishes from the problem and one is left only with the $\mathbb{W}$ -phase. [2326.1.5] Inserting eq. (10a) into eq. (9a) and using eq. (11) gives the Richards equation [20]

$\phi\frac{\partial{S_{{\mathbb{W}}}}}{\partial t}+\mathbf{\nabla}\cdot\left\{ k\;\lambda _{\mathbb{W}}({S_{{\mathbb{W}}}})\left[\mathbf{\nabla}{P_{\mathrm{c}}}({S_{{\mathbb{W}}}})+\varrho _{{\mathbb{W}}}g\mathbf{e}_{z}\right]\right\}=0$

(19a)

for saturation or

$\phi\frac{\partial\theta({P_{{\mathbb{W}}}})}{\partial t}=\mathbf{\nabla}\cdot\left[k\;\lambda _{\mathbb{W}}(\theta({P_{{\mathbb{W}}}}))\left\{\mathbf{\nabla}{P_{{\mathbb{W}}}}-\varrho _{{\mathbb{W}}}g\mathbf{e}_{z}\right\}\right]$

(19b)

for pressure after writing ${S_{{\mathbb{W}}}}(\mathbf{x},t)={P_{\mathrm{c}}}^{{-1}}(-{P_{{\mathbb{W}}}}(\mathbf{x},t))=\theta({P_{{\mathbb{W}}}}(\mathbf{x},t))$ with the help of eq. (11). [2326.1.6] To define the nonlinear function $\theta({P_{{\mathbb{W}}}})$ the typical sigmoidal shape has been assumed for ${P_{\mathrm{c}}}(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 ${S_{{\mathbb{W}}}}$ and global pressure $P$ . [2326.2.3] They differ also from the formulation in terms of ${S_{{\mathbb{W}}}},{S_{{\mathbb{O}}}},\mathbf{v}_{{\mathbb{W}}},\mathbf{v}_{{\mathbb{O}}},{P_{{\mathbb{W}}}},{P_{{\mathbb{O}}}}$ 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.

Authors:	R. Hilfer and R. Steinle
originally published in:	European Physical Journal: Special Topics, Vol. 223, pages 2323-2338 (2014)
originally submitted:	August 18th, 2014