This equation is relevant because of its influence on multiphase models. Sudicky and Huyakorn [1991] reviewed developments in this area through 1990. For present purposes, the paper of Celia et al. [1990] is recalled, which advocated the mixed formulation [ Milly, 1985] of Richards' equation, in which storage and transport terms are expressed in terms of water content and head, respectively. This avoids mass-conservation errors of the traditional head-based form and problematic degeneracies that appear in the content-based form as saturated conditions are approached. A maximum principle, assuring absence of nonphysical oscillations in results, was enforced by lumping the finite-element storage matrix in analogy with finite-difference approximations at the cost of some artificial smearing of the solution. Based on this experience, Celia and Binning [1992a] extended the mixed form to a mass-conservative model of two mobile phases, each consisting of one component with no mass transfer. The nonlinear implicit equations were solved for two phase pressures (p-p formulation) by a modified Picard iteration. This can be viewed as a Newton-Raphson linearization with incomplete Jacobian, including primary-variable dependencies of spatial and temporal derivatives, but not of nonlinear coefficients. Time-step sizes were selected so as to expect 10 to 15 iterations for nonlinear convergence. Another study [ Celia and Binning, 1992b] noted that, if an infiltrating moisture front is moving at near-constant speed, it is important to acknowledge the wave-like nature of the two-phase flow problem and base the model on one pressure and one saturation (p-s formulation) with a fractional-flow function of saturation. This permitted application of an Eulerian-Lagrangian approach [ Dahle et al., 1992] that approximates advection-dominated transport more accurately than standard Eulerian methods. Kueper and Frind [1991a,b] pointed out the ability of the p-s approach to handle infiltration into contaminant-free soil, avoid fictitious saturations, and specify boundary conditions in a straightforward way. Newton-Raphson linearization was used for their fully implicit equations.
The head- or pressure-based formulation was studied extensively by Abriola and Rathfelder [1993] for two mobile phases, seeking to control mass-balance errors in p-p formulations. With an unlumped finite-element storage matrix, they had difficulties because the chord-slope approximation of the capacitance coefficient [ Cooley, 1983] would not conserve mass, and elaborate schemes were presented to attempt to eliminate the resulting errors. A further problem was the degeneration toward zero of the capacitance coefficient as capillary pressure became small.
The fundamental difficulties of the widely-used head-based and p-p formulations [e.g., Osborne and Sykes, 1986; Kaluarachchi and Parker, 1989] are avoided in a straightforward fashion by moving to mixed and p-s forms. From a conservation point of view, the transport terms are naturally expressed in terms of potential gradients or differences, while the storage terms naturally involve contents or saturations. Passing to a head-based formulation requires the unnatural conversion of a change in saturation into the product of a capacitance coefficient and a change in head. Evaluation of this coefficient by a chord-slope formula makes the product equal to the change in saturation, effectively returning to the mixed form, as long as the storage matrix is lumped so that distinct capacitances at adjacent nodes do not enter; thus, the approach of Cooley [1983] was equivalent to the mixed form.
Because pressure typically behaves in a diffusive manner, as represented by a parabolic partial differential equation, a p-p formulation hides the wave-like or hyperbolic parts of the overall behavior of multiphase flow. Essentially, this is an ill-conditioned choice of primary variables when the two pressures are closely related, i.e., when capillary effects are small relative to the system at hand. This manifests itself in the difficulties cited above, because the wave-like behavior is present but cannot be well represented. The p-s formulation does not suffer from this drawback. Presumably, the p-p approach grew out of the head-based form of Richards' equation. With one dynamic phase, capillary pressure, or equivalently head, can be an effective surrogate for saturation. With two, especially two liquid phases, ill-conditioning is likely. Petroleum reservoir simulation, based on the different tradition of Buckley-Leverett theory, has long used analogues of the mixed and p-s forms [ Peaceman, 1977].