[page 509, §1]
[509.1.1] Consider a onedimensional, homogeneous, isotropic and incompressible porous column filled with two immiscible Newtonian fluids. [509.1.2] In a one dimensional model transversal variations and column thickness are neglected. [509.1.3] Let denote the saturation of wetting fluid (called water), and the saturation of nonwetting fluid (called oil). [509.1.4] Here time is , and is the position in a column of length . [509.1.5] Each of the two fluid phases is considered to consist of a continuous, mobile, percolating subphase, and of a discontinuous, isolated, trapped or nonpercolating subphase as discussed in [23, 28, 25, 24, 26]. [509.1.6] Following the notation of  the percolating phase of water is [page 510, §0] indexed by and its nonpercolating phase is indexed by . [510.0.1] The water saturation is then . [510.0.2] The percolating oil phase is indexed as and its nonpercolating phase by . [510.0.3] The oil saturation is . [510.0.4] The volume fraction of phase is defined as where the volume fraction of the pore space, also called porosity, is . [510.0.5] The volume fraction of the solid matrix is denoted as . [510.0.6] Volume conservation requires
[510.1.1] The mass balance of fluid phase can be expressed in differential form as
where are mass density, volume fraction and velocity of phase as functions of position and time . [510.1.2] is the mass transfer rate from all other phases into phase .
[510.2.1] The momentum balance is written as ()
where is the stress tensor in the th phase, is the body force per unit volume acting on the th phase, is the momentum transfer into phase from all the other phases, and denotes the material derivative for phase . [510.2.2] The stress tensor for the nonpercolating phases is defined as the momentum flux across surfaces in the threedimensional continuum (see  for more discussion).
[510.3.1] For a macroscopically homogeneous porous medium
is assumed. [510.3.2] Incompressible fluids are assumed so that their densities
are independent of .
[510.4.1] The percolating and the nonpercolating phases are able to exchange mass through breakup and coalescence of droplets, ganglia and clusters. [510.4.2] The mass transfer rates must depend on rates of saturation change. [510.4.3] They are here assumed to be
[page 511, §0] where are constants. [511.0.1] The parameters , are defined by
where are limiting saturations for . [511.0.2] In eq. (7) the shorthand is used and
denotes the Heaviside unit step function. [511.0.3] Equation (7) follows from the form used in  for small rates of saturation change. [511.0.4] The mass exchange depends on the sign of . [511.0.5] The sign determines the type of process. [511.0.6] It can switch locally between drainage and imbibition. [511.0.7] This results in hysteresis. [511.0.8] The structure of the mass exchange term was chosen with hindsight such that theoretical results obtained in the residual decoupling approximation agree with experimental measurements of capillary pressure. [511.0.9] The mass exchange terms have recently been further generalized to reproduce not only capillary pressure, but also experimental capillary desaturation curves .
[511.1.1] Turning to the momentum balance, note first that the inertial term will not be neglected in this paper. [511.1.2] The stress tensor for the four phases are specified as
where and are the fluid pressures in the percolating phases. [511.1.3] The constants and exponents are constitutive parameters. [511.1.4] The assumptions for the nonpercolating phases reflect their modified pressure. [511.1.5] This phenomenon seems to have been observed in experiment [2, Fig. 2, p. 233]. a (This is a footnote:) aFig. 2 in  shows that the pressure measured by the pore pressure transducers PPT3 and PPT4 rebounds after the end of the infiltration, i.e. when the DNAPL has passed and the transducers measure the pressure of water and residual PCE [511.1.6] In applications, the parameters and are determined by measuring capillary pressure curves (see below). [511.1.7] The body forces are assumed to be given by gravity and capillarity. [511.1.8] They are specified as
with constitutive constants and exponents . [511.1.9] The angle is the angle between the direction of the column and the direction of gravity with [page 512, §0] corresponding to alignment. [512.0.1] In applications the parameters and are determined by measuring capillary pressure curves. [512.0.2] The capillary body forces in eqs. (10) reflect the wetting properties of the medium. [512.0.3] They oppose gravity and reduce buoyancy driven flows of the disconnected phases. [512.0.4] This is illustrated in the figures below.
[512.1.1] Finally, the momentum transfer terms are assumed to be given by linear viscous drag characterized by constitutive resistance coefficients through the equations
where and was used because there is no common interface and hence no direct viscous interaction between these phase pairs. [512.1.2] Remember that the index represents the rock matrix. [512.1.3] For more details on these constitutive assumptions the reader is referred to the original papers [24, 25, 26].
[512.2.1] The balance laws (1b), (2) and (3) combined with the constitutive assumptions given above provide 9 equations for the 10 unknowns with . [512.2.2] To close the system of equations the conditions or could be used. [512.2.3] These conditions apply when the nonpercolating phases are immobile as it is often observed in experiment. [512.2.4] It turns out, however, that there exists a less restrictive and, in our opinion, more natural selfconsistent closure.
[512.3.1] The selfconsistent closure condition used in this paper follows naturally from many limiting cases. [512.3.2] One such limit is the residual decoupling approximation close to hydrostatic equilibrium described in detail in [26, Section 5., p. 216ff]. [512.3.3] A second, more general limiting case is the limit of vanishing velocities, i.e. for . [512.3.4] Here we formulate the selfconsistent closure condition in its most general form as
[page 513, §1] [513.2.1] This condition follows selfconsistently from the constitutive theory. [513.2.2] It expresses the experimental observation that the pressure difference depends more strongly on saturations than on velocities, and that it remains nonzero even for vanishing velocities. [513.2.3] Adding eqs. (3) for and and subtracting eqs. (3) with and from the result gives
[513.2.4] In this form the selfconsistent closure has been used in the numerical calculation below.
[513.3.1] The mathematical model defined above was introduced and discussed in [25, 24, 26]. [513.3.2] It was recently extended to include surface tension . [513.3.3] Notation and model formulation in this paper follow . [513.3.4] Approximations and analytical solutions for some special cases were given in [25, 24, 26]. [513.3.5] Here the system of equations will be solved by numerical methods. [513.3.6] To this end initial and boundary conditions are needed and will be discussed next.