Modeling multiphase migration of organic chemicals in groundwater systems--a review and assessment.

Over the past two decades, a number of models have been developed to describe the multiphase migration of organic chemicals in the subsurface. This paper presents the state-of-the-art with regard to such modeling efforts. The mathematical foundations of these models are explored and individual models are presented and discussed. Models are divided into three groups: a) those that assume a sharp interface between the migrating fluids; b) those that incorporate capillarity; and c) those that consider interphase transport of mass. Strengths and weaknesses of each approach are considered along with supporting data for model validation. Future research directions are also highlighted.


Introduction Background
With nearly half of the U.S. population drawing its drinking water from subsurface water supplies, the issue of organic contamination of groundwater has elicited considerable concern from both the public and private sectors. Primary sources of organic chemical pollution of the subsurface include: seepage from unlined lagoons and other surface impoundments; improper landfilling; improper surface and subsurface disposal practices; leaks in pipes, storage tanks, and other processing or transport equipment; and accidental spills (1). Many of the organic liquids emanating from such sources have low water solubilities and may migrate through the subsurface environment as nonaqueous-phase liquids (NAPL) or as organic vapors.
A number of case studies have been cited that attest to the importance of multiphase flow in the distribution and persistence of contaminants in the subsurface. The earliest case histories pertaining to nonaqueous-phase spill events are those that appeared in the European literature in the late 1960s and earlyto mid-1970s (2)(3)(4)(5). These papers document European experience with the accidental and controlled release of fuel oil. Prior to about 1980, there were few documented NAPL case studies in U.S. literature. 'Those case histories that were published also focused on spills of fuel oil (6) or gasoline (7). Subsequent to about 1980, however, a number of field cases have been cited that involve the release of denser-than-water organic liquids. The most detailed *Department of Civil Engineering, The University of Michigan, Ann Arbor, MI 48109-2125. of these are involving a transformer oil spill in Saskatchewan, Canada (8)(9)(10), and NAPL migration at a number of landfill sites near Niagara Falls, NY (11)(12)(13)(14). In these case studies, significant vertical and lateral transport of NAPL is documented. Other representative case studies involving multiphase transport may be found in (11,(15)(16)(17)(18).
NAPL infiltration and migration within the subsurface has been described qualitatively by a number of authors. The earliest description appears to be that of Van Dam (19). He examines a fuel oil spill scenario and analyzes the processes controlling its infiltration and subsequent migration along the water table. Schwille has presented a number of general discussions of the NAPL migration process (20)(21)(22). He explores dissolution and volatilization as well as nonaqueous-phase migration as contaminant transport pathways. In his more recent works, Schwille presents a description of the migration of the heavier-than-water organics. A general overview of the NAPL infiltration process can also be found in Mackay et al. (23). For a detailed presentation of the underlying physical processes and principles governing multiphase flow the reader is referred to the English translation of the work by Albertsen et al. (24).
Consider the spill scenario shown in Figure 1. This cross-sectional schematic depicts the distribution of contaminants in the subsurface for a spill of a lighter-thanwater organic. On introduction into the subsurface, the organic will migrate as a distinct phase downward under the influence of gravity through the unsaturated zone. This vertical migration will also be accompanied, to some extent, by a lateral spreading because ofthe effect of capillary forces. As the spill progresses downward  (111). DIFFUSION ZONE -(soluble components) through the unsaturated zone, it will leave some residual organic trapped in the pore spaces. This entrapment of fluid is due to surface tension effects. In addition to the migration of the nonaqueous phase, some of the organic may volatilize and form a gaseous envelope of organic vapor extending beyond the main zone of contamination.
If the spill is sufficient in size, some ofthe contaminant will eventually reach the saturated zone. Here its behavior is governed chiefly by density. If it is less dense than water, the NAPL will spread laterally along the capillary fringe forming a lens or pancake. It may also depress natural groundwater levels. As it encounters the flowing water phase, soluble components may dissolve to form a contaminant plume that can then migrate under the hydraulic gradients present in the aquifer. If the organic is more dense than water, the water table will offer no impediment to its downward migration. It will displace the water and continue its vertical migration until it encounters some impermeable stratum (Figure 2). Here it will then move under its own pressure and gravity forces along the aquifer bottom. As in the unsaturated zone, some of the organic will be held in the pore space. In the saturated zone, this residual organic will serve as a pollution source to the flowing groundwater. In the unsaturated zone, infiltrating rainwater may dissolve either the organic vapors or the residual organic liquid and transport these organic components to the saturated region.
Over the past two decades, a number of mathematical models have been developed that describe some aspect of the multiphase transport scenario outlined above. This paper examines the state-of-the-art with regard to such modeling efforts. Models developed to describe the nonaqueous-phase transport of organics in the subsur- face as well as those designed to address related issues such as organic phase volatilization or dissolution are reviewed. Those models that deal solely with solute transport are excluded. Of potential interest to readers are two recent reports which present discussions of some of the models reviewed herein (25,26). The first report was written for the Department of Energy and intended to explore the state-of-the-art in the knowledge of NAPL transport. The second study, conducted for the Leaky Underground Storage Tank (LUST) Office of the EPA, evaluated some existing models for their applicability to tank leakage assessment. In this paper, models are presented in a framework that explores their mathematical underpinnings and the assumptions inherent in each approach. The extent to which individual models have been validated is examined and future directions for research are highlighted. The main body of the paper is divided into three major sections: The second section presents the concepts and equations fundamental to all the models reviewed. After a brief overview of related modeling experience in the oil industry, the third section describes in detail individual multiphase flow and transport (MPFT) models. The fourth section highlights the model validation studies which have appeared in the literature to date. Some of the published experimental investigations that are relevant to the verification of significant model assumptions are also presented. organic, or a subset of these. A species can be defined as a specific chemical component that is present in one or more phases or could be considered as a group of such components if some average characteristics of the compounds could be defined. The mass balance equation can be written as: at R (6aP g') + V * (EafpaW,avya) _ V * i = S'a + Ri a (1) Here va is the mass average velocity of the a phase wia is the mass fraction of species i in the a phase Ca is the fraction of volume occupied by the a phase pa is the intrinsic mass density of the a phase

Equation Development
The development of the equations governing multiphase flow and the transport of organic chemicals has its roots in continuum mechanics and mixture theory. A large body of literature exists that pertains to multiphase flows and the development of balance laws. A review of this material, however, is beyond the scope of this work. The interested reader is referred to the comprehensive review article by Bedford and Drumheller (27), which discusses some recent advances in this area. Other relevant papers include (28)(29)(30)(31).
Multiphase contamination of groundwater systems is just one example of a multiphase mixture scenario. In the description of such multiphase flows in porous media, two balance equations are of primary interest. These equations are the mass and momentum balances. Although the other balance equations, including energy and entropy, can be used to extract further information about the system, existing multiphase models focus on these first two balance laws. Model constitutive relations are commonly empirical in nature and have not been derived from any rigorous application of constitutive theory. Mass and momentum balances are presented and discussed below in an effort to present the mathematical framework for existing MPFT models. For consistency, a modified form of the notation employed by Abriola and Pinder (32) will be used throughout this paper.

Mass Balance Equations
The easiest way to develop any of the governing equations currently used in MPFT models is to start with one basic building block, the mass balance equation for a given species i in a phase a. Here a stands for any phase present in the groundwater system. Existing models deal with four such phases: soil, gas, water, and Ji' is the non-advective flux of a species i in the a phase Sia represents the exchange of mass of species i due to interphase diffusion and/or phase change Ric represents an external supply of species i to the a phase The first term in Eq. (1) accounts for the accumulation of the mass of species i in phase a. The second term accounts for the movement of mass due to advection of the phase. Motion due to nonadvective effects (such as diffusion and dispersion) is accounted for by term 3. The first term on the right hand side of Eq. (1) can be thought of as a source/sink term due to a ch,nge of phase. R>, represents the destruction or creation of the species due to chemical reaction or degradation. Eq. (1) is subject to the following constraints: which follow from the definition of mass and volume fraction. It is also apparent that for a massless interface: That is, when mass is lost by one phase due to interphase exchange, an equal amount of mass is gained by another phase.
Based on Eq.
(1), one may develop a mass balance relation for a specific phase by summing over all species present in the phase, or alternatively, one may sum over all phases to get a species mass balance equation. Most MPFT models employ one of these two approaches, although a model consisting of a combination of the two types of equations is also possible. Two classes of model equations will be explored here. The models that employ these equations and more subtle variations between approaches will be discussed in the next section. Immiscible Phase Flow Model Equations. If one is concerned primarily with bulk fluid motion, a set of equations can be developed from Eq. (1) that assumes that there is no mass exchange between phases and no chemical reaction or degradation. Eq. (1) is summed over all species to yield: Here use has been made of constraint [Eq. (2a)]. The nonadvective flux terms also sum to zero because they deal with relative motion of the species within a phase. Equations of the form (3) have generally been applied to model the migration of pure immiscible fluids, that is, fluids that contain primarily one chemical component or fluids whose physical properties can be regarded as spatially invariant.
In general, four equations in the form of Eq. (3) may be written, one for each phase: soil (s), gas (g), water (w), and organic (o). If one assumes that the porous medium is rigid (porosity is invariant in time), the soil equation need not be included. Similarly, if the gas phase is assumed to remain at atmospheric pressure, a common assumption in the modeling of unsaturated flow of water in groundwater systems (33), the gas equation can also be eliminated to yield: n-(a p) + v * (pQs6nsv) = 0 a = w,o (4) Here n is the matrix porosity and sO. is the saturation of the a phase (E^= ns,,a). If fluid compressibilities are neglected, one has finally: a n-gj(s,)+V.(s0,nv)=0 a =W,o (5) Compositional Model Equations. If the modeler is concerned with the interphase transfer of mass (i.e., the formation of a contaminant plume or the transport of organic vapors), it is reasonable to write balance equations for each chemical component of interest. Such an approach is termed compositional in the petroleum industry literature. The species balance equations are obtained by summing equations in the form of Eq. (1) over all phases to yield: {(p6 a,) + V . (pfaV'j') -V . J,Q} = (6) HeeontnEq ]as = g,iS,nco Here constraint [Eq. (6)] has been incorporated. For nonreactive components, the right-hand side of Eq. (6) will be zero. The number of equations in the form of Eq. (6) that are required depend on the number of chemical components of interest. Further simplifications are also possible. If the soil matrix is assumed rigid, as with the immiscible flow case, the soil species equation may be eliminated. It should be recognized that solution of a set ofgoverning equations in the form of Eq. (6) results not only in a description of the physical distribution of each fluid phase but also in the delineation of the composition of this phase at any point in space and time.

The Momentum Balance Equations
In order to solve equations in the form of Eq. (4) or Eq. (6), expressions must be sought for the mass average phase velocities vy. Relations for these velocities can be obtained in a theoretical framework by considering the momentum balances for each phase. Most of the theoretical work in this area has dealt with single fluid phase (saturated flow) systems. Through the application of linear constitutive theory, it can be shown that the momentum balance can be reduced to Darcy's law (34): where vI" is the fluid phase velocity relative to the solid phase (= v" -v) K is a second order symmetric tensor (7a) PO is the fluid phase thermodynamic pressure g is the gravity vector Eq. (7) is valid if inertial as well as (macroscopic) viscous effects can be neglected. These are reasonable assumptions for most groundwater systems.
Through empirical measurements it has been determined that the tensor appearing in Eq. (7) can be expressed as (33): Here k is termed an intrinsic permeability tensor and KU is the fluid dynamic viscosity. The intrinsic permeability tensor depends only on the porous medium and is independent of fluid properties.
Analysis of the momentum balance for the multiphase fluid case is a much more complex task. Recently, Whitaker (35) has investigated two-phase fluid flow in a rigid, porous medium. In his analysis, the fluids are assumed incompressible. Whitaker's development ofthe equations of motion for this case incorporates interfacial tension between the fluid phases. Using a number of 120 vors = -K -(VP' p'g) assumptions and scale arguments, he produces the following form for the momentum balance of phase a: va = Km . (VP" -p'g) + K v (9) Here K, and K. are second-rank tensors and vO is the velocity of the alternate fluid phase. The first term in Eq. (9) is ofthe same form as Eq. (7). A relation between K, and K, however, has yet to be determined theoretically. The second term on the right-hand side of Eq. (9) arises due to the viscous drag of one fluid upon the other. Whitaker suggests that this term might be small for many cases of interest, but a complete investigation of the importance of this term has yet to be undertaken.
All existing MPFT models use a momentum expression in the form of Eq. (9) with the second term dropped. The tensor K, is assumed to be a product of the single phase tensor K and a relative permeability (scalar) function, kr,. Justification for this assumed form of the momentum balance is experimental rather than theoretical. A number of experimental investigations in two-and three-phase fluid systems (36)(37)(38)(39)(40) support the use of this modified form of Darcy's law, which may be written as: a = =-kk ro .(VPa . pag) (10) jAaCor The relative permeability that appears in Eq. (10) is a function of saturation or capillary pressure P. the difference between the two phase pressures (Pa = p' -PO). Relative permeability can be determined in the laboratory or predicted by the application of some semiempircal model. Some of these models are discussed in the fourth section. It is well known that the krc coefficients may also exhibit hysteresis. That is, these coefficients have been found to depend on saturation history; they are not single-valued for a given saturation. Whitaker (35) suggests that the form of his Eq. (9) may be better able to account for such effects. tinuous in space. Due to varying matrix properties, saturation may be discontinuous and solution for such a variable can cause difficulties in a numerical model. Eq. (11) contains a number of undetermined parameters. Porosity n and intrinsic permeability k are assumed known properties of the matrix. Customarily viscosity I,, a weak function of pressure, is assumed constant. This leaves the following parameters for evaluation: 8aI pa, kra (Ila) For the incompressible fluid case, pa will be a constant. In general, the density will depend on the fluid pressure, Pa. Thus, fluid density may be expanded in terms of fluid pressure by incorporation of P', the compressibility of the a phase. Compressibility is defined as: -1 dpax -pa dP(1 For slightly compressible fluids, a is essentially constant.
In MPFT models, saturation is assumed to be some known function of the capillary pressures. This function must be empirically determined. As in the case of relative permeabilities, the saturation function is not single-valued; it exhibits hysteretic behavior.
The first term in Eq. (11) can thus be expanded in terms of the unknown pressures by incorporating the definition [Eq. (12)] and differentiating saturation with respect to capillary pressure(s): a (scp)-n8'pm apa +np 9a 'c n&ãt~p flBaf3apaPC + (13) Note that for a three-phase fluid system, s8 may depend on two capillary pressures.

Constitutive Relations
Modified Darcy's law expressions ofthe form (10) may be substituted into a system of mass balance equations, such as Eq. (4) or Eq. (6), to yield the equations governing the multiphase flow of fluids in a porous medium. As an example, consider the immiscible flow of fluids in a rigid matrix. Incorporation of Eq. (10) into Eq. (4) yields: n a pa) -V k* [paA; * (VP, -pag)] = 0 (11) Eq. (11) can be thought of as a partial differential equation in the unknown P'. Pressure is the most common choice for the unknown parameter since it is con--*. [kk/, * (VPopag)] = 0 (14) at = 0, W In formulating Eq. (14) it has also been assumed that products of density and pressure gradients can be neglected (since both have very small magnitudes in porous media flows). Eq. (14) is the equation governing the flow of the immiscible and slightly compressible fluid, a, in a rigid porous matrix. It is a nonlinear partial differential equation because of the dependency of coefficients s , asjO)PR and k., on the unknown pressures.
A similar equation can be written for each fluid phase. Equations are coupled through the capillary pressure term [2nd term in Eq. (14)] and the dependence of the nonlinear coefficients on capillary pressures.
The compositional model Eq. (6) can be manipulated in an analogous manner to yield the governing equations for multiphase systems with mass exchange. Here, however, there are a number of additional complicating factors because of the dependence of phase properties on composition. For example, fluid viscosity as well as density must be specified as functions of composition. There are also two additional equation terms, which must be evaluated: those accounting for nonadvective flux and chemical reaction.
The nonadvective flux term is commonly assumed to have a Fickian form. This is the form employed in most existing solute transport models. It is intended to account for both molecular diffusion and hydrodynamic dispersion effects. The nonadvective flux is represented as: J, = pa.eDa* Vwic (15) where D" is a macroscopic second-order tensor which is a function of molecular diffusion and fluid velocity. A discussion of the appropriateness of using the form (15) is beyond the scope of this paper. Note, however, that a great deal of research is presently being directed in this area (41)(42)(43)(44)(45)(46)(47).
There is no general functional form which may be specified for the reaction terms which appear in Eq. (6). These must be formulated on a per species basis. Simple decay rates could be directly substituted into Eq. (6) or a system of chemical reaction equilibrium or rate equations could be solved to determine these terms.
Finally, to close the system of governing equations in a compositional model, expressions are needed to relate mass fractions of a given species within all the phases. Existing models employ the assumption of local equilibrium to develop these relations. Local equilibrium implies that within some short-time scale (essentially instantaneously) contiguous phases reach a thermodynamic equilibrium. Thus, the mass fraction of a species in one phase can be related to the mass fractions of this same species in other phases via partition expressions of the form: wia = K-wiP (16) where Kia is called the partition coefficient of species i between the a and ,B phases. In general, partition coefficients should be represented as functions of phase compositions and pressures. These coefficients may be determined from Henry's law constants and solubility relations.
Incorporation ofthe modified form of Darcy's law [Eq.  (17) Eq. (17) can be treated as a partial differential equation in the unknown phase pressures and mass fractions. One equation may be written for each species. The first term in this equation can be expanded in terms of pressures by incorporating compressibilities and saturationpressure derivatives in a similar manner to that employed in developing Eq. (14) for the immiscible flow case. Supplemental relations in the form of Eq. (16) as well as representations of the reaction terms must be available to close the system.
It should be noted that if one considers a region in which only the water phase exists, and the water phase mass balance equation in the form of Eq. (3) is incorporated into Eq. (17) above, one has the familiar solute transport equation: ici +v.Vci-V.(D -Vci)= Riw (18) where ci = p'wx,i is the component i concentration in the water phase. Here, it has also been assumed that the presence of the solute does not significantly affect water phase density. A similar form can also be developed for regions in which only the gas phase is present.

Review of MPFT Models
The Oil Industry Experience In this section, MPFT models which have been developed and/or applied to the prediction oforganic chemical migration in the subsurface are presented and evaluated. Before these models are discussed, however, it will be instructive to briefly review the experience of the oil industry in modeling similar multiphase flow problems. A thorough review of this material is beyond the scope of this work. For more information the interested reader is referred to the article on the stateof-the-art in reservoir modeling by Coats (48) or the excellent texts by Aziz and Settari (49) and Peaceman (50). A recent review of multiphase flow modeling from a broader perspective, including current developments in the numerical area, is given in (51).
Over the last four decades, much of the research pertaining to multiphase flows has focused on describing fluid flows in underground petroleum reservoirs. Such reservoirs may contain natural gas as well as connate water and oil. Thus, the description of oil recovery involves three-phase mass balance equations of the type presented in the last section. Furthermore, the development of enhanced oil-recovery technologies (such as CO2 flooding) has led to interest in the interphase mass transfer of oil components. Species mass balance equations in the form of Eq. (17) govern such processes.
The earliest efforts at modeling oil recovery centered on the two-phase waterflood problem in which water displaces oil within the formation. Buckley and Leverett (52) introduced a simplified model that could describe such a two-phase displacement process. In this model, capillary effects are assumed negligible. This means that gradients of capillary pressure are omitted from the governing equations. In addition, the two fluid phases, as well as the porous matrix, are assumed incompressible. If water is injected at a constant rate, this results in a constant flowrate (q) within the formation (q = qW + q0 = constant at any point). Neglecting gravity effects (assuming a horizontal formation), two-phase flow equations based on Eq. (11) may be written for one-dimensional displacement: Since the total flow rate is known and constant and sO, = 1 -sW, only one of the above equations needs to be solved. Eq. (19b) may be rewritten as: 988w a y-) =0 (20) where fw K,JKQ + K, is called the fractional flow of water. Eq. (20) is known as the Buckley-Leverett equation. It is a hyperbolic, nonlinear, partial differential equation in saturation. The nonlinearity arises because fw depends on the unknown fluid saturations. There is a great deal of literature that focuses on the solution of the Buckley-Leverett problem (49,(53)(54)(55). Experience indicates that a physically realistic solution can only be obtained by introducing a second-order dissipative term or artificial capillarity into the equation. This is achieved in the petroleum literature by a procedure known as upstream weighting. For presentations of the upstream weighting concept see (56,57).
As shown later in this section, the concept of negligible capillarity has been applied to the modeling of organic infiltration into the unsaturated zone and organic spreading on the water table. The one-dimensional nature of the Buckley-Leverett problem and the assumption of constant total fluid flow, however, appear to severely limit its applicability to problems involving immiscible contaminant migrations in groundwater systems.
In order to deal with more complex, multi-dimensional, nonsteady, three-phase flow situations, a class of models commonly called "black-oil simulators" have been developed in the petroleum literature. Examples of such models can be found in (49,50). The roots of this modeling approach lie in the species mass balance equations given in Eq. (17). Three species are modeled: water, stock-tank oil, and stock-tank gas. Whether a component is part of the stock-tank oil or gas is determined by the phase in which it appears at the surface under stock-tank conditions (STC). In pressurized reservoir formations, some of the stock-tank gas commonly partitions into the oil phase. All other interphase par-titioning (such as oil dissolution in the water phase), as well as the nonadvective flux terms or reaction terms appearing in Eq. (17) are neglected. Based on these assumptions, species mass balance equations are manipulated to produce a system of three, nonlinear partial differential equations that describe the flow of threefluid phases (water, oil, and gas) in a petroleum reservoir.
Two methods are commonly employed for the solution of this system. The first, introduced by Sheldon et al. (58) and Stone and Garder (59), is called the implicit pressure-explicit saturation (IMPES) method. This method is implemented by producing a new equation (called the "pressure equation"). The pressure equation is formed from a linear combination of the governing equations that are combined in such a way as to eliminate the time derivatives of saturation from these equations (recognizing that s, + sW + Sg = 1). Here, the critical assumption is that the change in capillary pressures over the time step is negligible, and capillary pressures are treated explicitly in the formulation. The pressure equation is solved implicitly for the organic phase pressure; then phase saturations are computed from the explicit solution ofthe water and oil equations. Although this solution technique is conceptually attractive because it involves the implicit solution of only one equation, it suffers from stability restrictions on the size of the time step which may be employed.
The second technique for the solution of the black-oil governing equations is termed the simultaneous solution (SS) method. As its name suggests, this approach involves the implicit solution of all three equations. It was introduced by Douglas et al. (60) and then extended and further analyzed in (61,62). To solve the equations, the time derivatives are expanded in terms of fluid pressures, as was done to produce Eq. (14). Due to its implicit solution technique, the SS method has greater stability than the IMPES method. Additionally, capillary pressure need not be lagged. It should be noted, however, that the stability of this approach is affected by the treatment of the nonlinear coefficients that appear in the convective flux terms. An implicit treatment of these coefficients, which was first implemented by Blair and Weinang (63), produces the greatest stability.
Models similar to the black-oil simulators have been developed to model the two-and three-phase flow of contaminants in the subsurface. These models are discussed later in this section. Since capillary pressures play an important role in such systems, the SS method has been uniformly employed to solve these equations. Because the black-oil models are based on the assumption that the partitioning of the oil into the water phase can be neglected and that nonadvective fluxes (diffusion/dispersion) are unimportant, these models are not suitable for the examination of pollutant partitioning and transport in groundwater systems.
A final group of models, termed compositional simulators, have been developed by the oil industry to examine the chemically complex scenarios involved in tertiary oil recovery. Early work in this area was performed by Roebuck et al. (64). More recent examples of compositional simulators are presented in (65)(66)(67)(68)(69). Compositional models are based on a system of equations which govern the transport of the chemical species of interest. These equations have the general form given in Eq. (17). Such equations are supplemented with equilibrium partition relations or equations of state. Paralleling the black-oil models, two methods of solution are common: a sequential approach based on the IMPES technique, and the SS method. In these models, instantaneous local equilibrium is assumed between the phases, and oil partitioning into the water phase is neglected. As with the black-oil models, typical formulations neglect diffusion and dispersion of components within a phase. In the sequential approach, capillary terms are evaluated at the previous time level. (Capillary pressure changes over the step are neglected.) With modification, the compositional approach is applicable to problems of groundwater contamination. It is capable of describing contaminant partitioning and transport in the subsurface and lends itself well to the incorporation of chemical reaction and degradation effects. A small number of models which are based on this approach have appeared in the groundwater literature and are discussed at the end of this section.
Before we move on to examine groundwater models, we should mention that the vast majority of the models that have appeared in the petroleum literature employ the finite difference method in solving their equations. This is primarily due to the fact that problems have been encountered in applying finite element techniques to the solution of the immiscible flow equations (54). As some of these computational difficulties are surmounted, the finite element technique is gaining popularity. [For recent examples of the application of this method to reservoir simulation, see (70)(71)(72)(73).]

MPFT Models in Hydrology
In the remainder of this section, the MPFT models that have been developed to examine the multiphase transport of contaminants in groundwater systems are presented and discussed. For coherence, these models have been grouped into three categories. The first group is primarily comprised of analytical and semi-analytic models which were developed in Europe in the early to late 1970s to assess the behavior of fuel oil spills in the subsurface. All of the models within this category treat the immiscible flow of contaminants as a sharp interface problem, much as the Buckley-Leverett approach to the water flooding of petroleum reservoirs discussed in the last section. The second group consists of recent computer models that also deal with the immiscible flow of organics, but in these models, the sharp interface approximation is relaxed (capillarity is incorporated). The last group of models includes those that incorporate the interphase transfer of mass. Simple approaches to the description of aqueous phase and vapor-phase transport are also included in this category. Sharp Interface Models. Van Dam (19) presented the first detailed analysis of the nonaqueous phase organic transport problem. In his paper, the equations governing the three-phase immiscible flow of fluids in a porous medium are presented. The work emphasizes the importance of capillarity to the migration of contaminants, especially in heterogeneous, anisotropic media. Using an analogy between the porous matrix and a collection of capillary tubes, Van Dam discusses capillary rise, fluid wettability, and residual saturations. Assuming that water is the wetting phase in groundwater systems, Van Dam uses considerations of equilibrium potential to develop an approximate expression for the lateral extent of petroleum hydrocarbon (oil) migration on the water table. His equation has its roots in the expression developed by Leverett (74) for the approximation of the capillary rise in a porous medium based on a capillary tube analogy: Apg AApg hc is the average capillary rise Pc is the average capillary pressure between the two fluids B is a wettability constant characteristic of the fluids a is the surface tension between the fluids Ap is the density difference between the fluids Van Dam suggests that at equilibrium the thickness of the oil zone upon the water table will be given by 2(hcX)g where (h,) is the capillary rise of oil in an oil-gas system. Within this oil zone, the oil saturation will vary between zero and its residual value (sro). Thus, an average saturation in the zone may be assumed to be '/2Sro. Combining these assumptions with Eq. (21) gives rise to the following expression for the lateral extent of an oil lens: where Q is the volume of infiltrated oil and A, is the horizontal area of the oil lens. The quantity in brackets is the volume of oil contained per unit vertical crosssectional slice of the lens. Van Dam offers Eqs. (22a,b) as a good approximation for determining the oil extent along a sloping water table. In this case, however, the lens will be displaced in the direction of flow. Eqs. (22a,b), although attractively simple, have a number of limitations. The most obvious is that this type of analysis is restricted to contaminants whose specific gravity is less than one. In addition, the approximation is valid only for the calculation of an equilibrium contaminant distribution in a homogeneous, isotropic me-dium where the water table has been stationary for the duration of the spreading. The analysis is conservative in that it neglects the organic phase which is held up as a residual in the unsaturated zone above the water table.
Following Van Dam's work, a number ofinvestigators developed analytical and semi-analytical approximations to describe the time evolution of a petroleum hydrocarbon phase spill (75)(76)(77)(78). To simplify the analysis, the migration of an oil body was divided into two major stages: vertical infiltration, and lateral spreading along the capillary fringe. Individual equations were developed for each stage. The method of analysis and the major assumptions inherent in the resulting approximations are presented below.
As with the numerical approaches described in the following sections, the semi-analytical approaches to oil spill analysis employ two basic principles: mass conservation and the validity of Darcy's law. In the semianalytical approach, however, the mass balance relation is implicit rather than explicit, i.e., a mass balance equation is not written. Instead, a time-distance oil profile is developed directly from Darcy's law. To make the problem more tractable, the saturation-capillary pressure relation between two fluids is idealized as a rectangle (Fig. 3). The assumption is also made that the oil displaces only a single fluid in its migration (air in the unsaturated zone; water during lateral spreading within the capillary fringe). This displacement of one fluid by another is assumed to behave as a piston flow, that is, the saturation proffle behind the displacement front is assumed uniform (Fig. 4). By making this assumption, one neglects front smearing due to capillary forces. The zone in which both fluids are mobile within the pore space is neglected. One would expect this piston-flow assumption to be more accurate in describing behavior within coarse-textured uniform soils or in situations where gravity or pressure is the predominant driving force.
Consider the vertical infiltration of an oil spill into the unsaturated zone. If the spill area on the surface is fairly large (on the order of square meters), the lateral migration of the oil in the unsaturated zone can be neglected (76 In writing the above equation, a rigid matrix and an incompressible organic fluid have been assumed. Organic phase pressure has been replaced by the sum of the gas pressure and the capillary pressure between the two fluids. If h, = PogIpg is defined as the "capillary head" and the gas phase pressure is assumed constant, one has: Suppose a steady seepage event is considered. Since the piston-flow assumption has been invoked, the saturation behind the infiltrating front is assumed uniform, and the coefficient multiplying the gradient term in Eq. (23b) will be constant. Steady seepage implies v°= constant, which in turn, leads us to conclude that ah,I dz must also be constant. Thus, the capillary head must be linearly distributed in the infiltrating organic body.
A reader familiar with the hydrology literature will recognize that this type of approximation is analogous to the commonly employed Green-Ampt model for water infiltration into the unsaturated zone (79). One can compute dahlaz from two known capillary pressures: the value at the soil surface, (hc)surface, and the (assumed) value at the front, hc (Fig. 3). Recognizing that v°d zldt, Eq. (24) becomes: Ponso case of organic redistribution. Here, however, saturation profile evolution is idealized as shown in Figure 6. Comparison with Figure 5 reveals differences in assumed curve shape. In Illangasekare and Reible's work two fronts-a wetting front and a drainage front-are identified. Behind the drainage front, the organic saturation is assumed to be at a residual level. As in the steady infiltration case, a linear head distribution within the saturated organic zone is assumed. An equation based on conservation of total mass of organic can be written: Eq. (24) is an integrable expression in space and time.
Using similar assumptions, integrable expressions can also be developed for the line and point source cases (77). A model recently presented in (80,81) closely parallels the work described above. An expression identical in form to Eq. (24) is developed and its integration is implemented in a numerical code that permits a variable organic head at the soil surface (variable ponding depth). Corapcioglu and Hossain (82) present a two-dimensional organic infiltration model, which is also based on the sharp interface assumption. Here, the displacement of water by organic is considered. Mass flow rates in the horizontal and vertical directions are assumed constant and independent. Employing these assumptions, the Buckley-Leverett Eq. (20) may be extended to a two-dimensional form. The resultant hyperbolic equation is solved using a modification of a standard method of characteristics (MOC) solute transport code.
Another infiltration case of interest is that of a sudden release of a limited oil volume (the oil redistribution problem). In treating this case, Mull (6,75,76) assumes that the capillary term in Eq. (23b) is negligible, i.e., he assumes that the redistribution process is governed by gravity. The coefficient multiplying the gravity term in Darcy's law will not be constant for this case since sO, as well as kro will decline with time as the oil migrates vertically. Oil saturation proffles are assumed to be rectangular as in Figure 5. Since the total mass of oil must remain constant, the saturation profile is inversely proportional to the distance of migration: files (77). z (24a) If relative permeability is some known function of saturation, it can, thus, also be expressed as a function of z. Hence one has from Eq. (23b): As is Eq. (24), this is an integrable expression. Oil migration is assumed to cease at a point z such that kro = 0 (s = Sro).
Illangasekare and Reible (80,81) also examine the To compute lateral spreading of the oil within the capillary fringe, Darcy's law can again be employed and an integrable expression developed. As in the vertical displacement case, piston flow is assumed. Here the capillary field, however, is axial. The driving forces accounting for fluid spread are the accumulated head of oil in the infiltration zone beneath the source and the water table slope. In this case, the magnitude of the driving force will decrease with time as the lens spreads on the water table. Radial spread is assumed to cease when the lens thickness at the source equals the capillary rise. An analysis of this situation is complicated by the fact that the oil head may also depress natural groundwater levels. The vertical displacement of natural groundwater levels can be computed from equilibrium considerations.
Schiegg (78) developed a computer code based on the above semi-analytical approach to describe an oil spill event. In his model, analytical methods are used to calculate the total time of infiltration and the maximum radius of extent for the infiltration phase. At the completion of this first phase, the oil body is thus assumed to be in the form of a cylinder of a known radius. This cylinder is used as an initial condition for computation of the lateral spreading phase of the spill. His program calculates the time to completion of migration and the maximum dimensions of the oil lens. A lateral spreading model is also presented in (83). This model is based on the principles discussed above. Here, however, the spreading oil body is divided into two regimes: an inner cylinder where pressure forces predominate and an outer cylinder of height h, within which capillary forces control the spread.
In the analytical models discussed above, the soil is assumed homogeneous and isotropic. Although in theory one could extend this approach to analyze layered soil systems, the computations become quite cumbersome. For many types of soils, it is anticipated that the piston-flow assumption will not be satisfactory. As with Van Dam's equilibrium analysis, the lateral spreading models are valid only for contaminants of specific gravity less than one. (Only the vertical infiltration analysis would be valid for heavier organics.) It is also evident that the division of the oil spill evolution into two phases is somewhat artificial since the second phase will begin before the first phase is completed. Despite the many drawbacks enumerated above, this semi-analytical approach can be useful in estimating the extent and rate of the oil spread in the subsurface. Models using this approach require less data than the models discussed in the next sections, and the data required can be more easily estimated. If one wishes to explore the effects of heterogeneities, subsequent infiltration episodes or sub-sequent displacement events (such as infiltrating rainfall or a rising water table), the semi-analytical approach is clearly inadequate. Similarly, semi-analytical models cannot be employed to examine the merits of proposed remedial schemes.
Another way of modeling the lateral migration of oil along the water table is to recognize the analogy between this system and that of a sharp interface salt water intrusion problem. Holzer (6) appears to have been the first to recognize this similarity. In his work, he employs an analytical solution developed by Hantush (84) for a decaying freshwater lens to compute the unsteady decay of an oil lens on the water table. An isotropic aquifer and horizontal water table are assumed. Although Holzer acknowledges the importance of residual saturation to the oil-spreading problem, his solution neglects this capillary effect.
A more sophisticated model based on the salt-water/ fresh-water analogy is presented by Hochmuth and Sunada (85). This is a two-dimensional areal model that describes the lateral migration of oil after reaching the water table. As in the semi-analytical approaches presented above, the capillary pressure/saturation curves are idealized as shown in Figure 3. The model employs these curves along with the principle of hydrostatics and mass and momentum balances to determine oil/air and oil/water interface elevations as functions of space and time. Model formulation and method of solution are described briefly below.
Consider the idealized distribution of fluids shown in Figure 7. If horizontal flow within each layer (organic, water) is assumed, the fluid pressures within each of these layers will be hydrostatic. Such an assumption is commonly employed in the description of groundwater flow in unconfined aquifers and is termed the Dupuit approximation (33). In this situation, an immobile water saturation is assumed to exist throughout the oil and air zones. Combining the idealized capillary pressure curves with this assumption of hydrostatic equilibrium yields the following two expressions for the interface heights: Here Darcy's law (Eq. 10) has also been incorporated in its potential form. Qa. are source/sink terms for the layers and Sy. represent the specific yield. This specific yield is defined as the mobile portion of the fluid saturation within the pore space of each layer. Eqs. (28a) and (28b) have the same form as the Boussinesq equation, commonly employed to model unsteady, horizontal flow in unconfined aquifers. The system of Eq. (27) through Eq. (28) represent a coupled set of equations in four unknowns: two fluid potentials h0,hw and two interface elevations (Zl,z2). In the Hochmuth and Sunada model, these equations are solved for each time level via an iterative technique. The flow Eqs. (28a) and (28b) are discretized and solved by an application of the Galerkin finite element method, using an implicit finite difference approximation for the time derivatives. Interface elevations are updated at the end of each time step. Constant head and impermeable boundary conditions are permitted. The resultant numerical model was verified by comparing a simulation to the analytical solution for the zero capillary pressure case. Some small scale laboratory experiments were also conducted which are presented in the fourth section of this paper.
Use of the Hochmuth and Sunada model is subject to some ofthe same restrictions inherent in the other sharp interface approximations discussed above. The specific gravity of the organic must be less than one, and although the formulation allows for heterogeneous hydraulic conductivities, the model is restricted to application for which the capillary pressure-saturation curves are uniform in space. As with any of the sharp interface models, the piston flow assumption restricts the application of this model to matrix materials for which the idealized capillary curve is a reasonable representation of the actual system behavior. For this numerical model, the investigators report computational difficulties (numerical oscillations) in the areas where the organic phase disappears (z1 = Z2). It should also be noted that no effort was made to incorporate a mass balance as a check on the numerical computations. As will be seen in the following sections, these observations are common to a number of computer models.
Immiscible Flow with Capillarity. Recognizing some of the limitations of the piston flow approximation, a number of investigators have developed computer models which can account for a transition zone where a number of phases are mobile (86)(87)(88)(89). These models solve immiscible phase flow equations of the form (11), coupled with constitutive relations for saturation and relative permeability. Such an approach is analogous to that employed in the so-called "black-oil" reservoir simulators discussed previously. Models of this type differ among themselves primarily in their dimensionality, in the number of mobile phases they can accommodate, and in their numerical solution techniques.
Two-phase (organic/water) models are presented by Little (86) and Osborne and Sykes (88). The Little model is a one-dimensional simulator that employs a finite difference technique to solve two-fluid phase flow equations in the form of Eq. (11). The Osborne and Sykes code extends this approach to two dimensions. Their model solves a system of flow equations for a vertical cross-section (x-z domain). Equations in the form of Eq. (11) are written for each phase (organic, water) and expanded in terms of capillary pressure as in Eq. (14). Fluid as well as matrix compressibilities are neglected, so that the first term in Eq. (14) is eliminated.
The resulting set of coupled second-order partial differential equations is solved in the Osborne and Sykes model by applying the method of weighted residuals. Linear, isoparametric, quadrilateral elements are employed. Time derivatives are approximated by finite difference discretization, and the model accommodates variable time-weighting. As in many of the oil industry models, upstream weighting is applied to the convective term. To achieve this, use is made of asymmetric weighting functions (90). The coupled nonlinear set of algebraic equations which results from the application of these numerical methods is solved in the computer model by a Newton-Raphson iterative technique. Both first-type (prescribed pressure) and second-type (prescribed flux) boundary conditions are incorporated. Nonhysteretic relative permeability and saturation constitutive relations are included in the model in tabular form.
Osborne and Sykes verified their code by a comparison with a one-dimensional simulation of the Little model. Reasonable agreement was found between the two simulations, considering model variations in equation solution techniques and methods of parameter evaluation. In their presentation, no mention is made of attempting to check the global mass balance performance of the model. Following model verification, the two-phase simulator was applied to a field site contamination event (the Hyde Park Landfill). An effort was made to qualitatively compare model predictions to the extent of nonaqueous phase organic migration at the site and to assess transport sensitivity to model parameters. No true validation could be attempted in this case because saturation and relative permeability relations for the site were unavailable. The situation was also complicated by the fact that the formation was highly fractured, and insufficient information regarding the width, frequency, and orientation of these fractures was available. Model simulations (within a reasonable range of uncertainty) failed to predict the depth or extent of nonaqueous phase migration actually encountered at the site. For this scenario, the predicted migration extent was shown to be quite sensitive to viscosity or intrinsic permeability and degree of anisotropy. Some limited numerical experimentation with capillary relations also indicated a strong sensitivity to this parameter. In addition, the authors conclude that the presence of anomalous high (or low) permeability zones could greatly influence the distribution of organic at the site.
Three-phase immiscible flow models are presented by Faust (87) and Kuppusamy et al. (89). Here the assumption is made that gas phase pressure is atmospheric so that, again, only two flow equations need to be solved. In the Faust model, equations in the form of Eq. (11) that also incorporate a source/sink term are discretized by the finite difference method. Matrix compressibility is included, but the fluids are assumed incompressible. A block-centered scheme is employed for a two-dimensional vertical cross-sectional domain. As with the Osborne and Sykes model, a Newton-Raphson scheme is used to compute the solution to the non-linear algebraic equations. Upstream weighting is also incorporated for the relative permeability terms, and a fully implicit approximation is used to ensure stability. Instead of solving for two pressures, the model computes water saturation and organic phase pressure. Remaining saturations and pressures are computed from the capillary relations. Faust's model was verified against two sub-problems: a one-dimensional linear waterflood problem, and a twodimensional scenario of water flow in the unsaturated zone. In the first problem, model simulations were compared to the Buckley-Leverett sdlution (no capillary pressure). Good agreement was found between solutions, with some smearing of the front in the numerical solution due to use of upstream weighting and a fully implicit formulation. In the second test problem, model solutions for water flow in the unsaturated zone were compared to solutions obtained with two other numerical models written specifically for unsaturated flow simulation.
The Kuppusamy et al. simulator is a finite element model analogous to that of Faust. Saturation expressions are expanded as in Eq. (14) to solve for pressures. Fluid and matrix compressibilities are neglected. As in the Osborne and Sykes model, quadrilateral elements and bilinear shape functions are employed as well as a variable-weighted finite difference discretization of the time term. No mention is made of upstream weighting. Solution of the resulting set of nonlinear algebraic equations is achieved by a relaxed Picard iteration scheme.
In the Kuppusamy et al. paper, model verification is not discussed, but an effort is made to validate the model with one-dimensional experimental data. (This validation is discussed in the next section.) Some heuristic evidence for convergence is given for this one-dimensional application of the model by considering the effect of time and spatial discretization on the solution.
In addition to the immiscible flow models described above, two other similar models are identified in a draft report prepared for the U.S. Environmental Protection Agency (26): a one-dimensional, three-phase flow model developed at Colorado State University (91) and a proprietary three-dimensional, two-phase flow code developed by Intera Technologies, Inc. This draft report also describes the un-documented three-dimensional extension of the Faust model (SWANFLOW). It should also be noted here that any unsaturated flow model could be adapted to treat organic immiscible phase migration in the unsaturated zone if only residual water is present in this zone. Oil reservoir models might also be used in such applications. Christensen et al. (92) report the adaptation of a three-dimensional, commercial, black-oil reservoir simulator (ECLIPSE) to describe subsurface migration of an oil spill. All of the models described in this section employ similar assumptions and have similar limitations. By incorporating capillarity, these models are capable of representing saturation distributions more precisely than the sharp interface models. Heterogeneous aquifer properties can be easily accommodated and attempts to simulate subsequent leakage or rainfall events can be made using this approach. None of the models incorporates hysteresis. There is no reason, which relates to model structure, that hysteresis submodels could not be incorporated into such codes if they were available. Recent publications propose such a submodel for the evaluation of hysteresis effects (93,94). To date, model dimensionality (of nonproprietary codes) is restricted to one and two dimensions. All models assume a passive air phase and neglect interphase mass transfer. An absence of global mass balance considerations is also apparent.

Interphase Mass Transfer
Parallel to the refinement of immiscible phase flow models, a number of models have been developed that consider the interphase partitioning of organics between the nonaqueous and water or vapor phases. All but one of the models which have been documented in this category neglect the migration of NAPL. Thus, one can think of these as treating a subproblem of the full multiphase flow scenario. Early investigations focused on describing the aqueous phase transport of soluble oil fractions emanating from an oil lens. Such models are the counterparts of the sharp interface models that treated lateral spreading of oil along the water table. Later efforts broadened the scope of these models to treat other organics of interest and to include vapor phase transport as well. More recently, abiotic and biotic transformations of the contaminants have been incorporated.
In early modeling efforts to describe solute transport away from a floating oil lens, analytical solutions to the transport equation [Eq. (18)] were employed. The chief problem then lay in the accurate characterization of the source term (specification of the boundary conditions). Hoffmann (96) views the nonaqueous organic phase as an oil lens floating upon the water table. This lens is assumed immobile for his computations and is treated as a point source. The strength of this point source (S) is proportional to what he terms the face of contact (F): The contact face is assumed in his analysis to be approximately equal to the undersurface area of the oil body (assuming a sharp oil-water boundary). Realistically, this contact face will depend on oil saturation within the lens. The constant of proportionality, E, in Eq. (29) is termed the "exchange constant." It is assumed to depend upon the type of organic fluid and the degree of contact. Degree of contact here implies some dependence of mass exchange on velocity. Hoffmann considers solute transport for the case of constant, unidirectional groundwater velocity and a homogeneous soil matrix. He suggests that point source solutions of the transport equation (18) for constant source strength can provide reasonable estimates for concentration evolution in time at some distance from the source.
A similar analytical approach is employed by Fried et al. (97) [also (98,99)]. Here, source strength is characterized by the expression: where C0 is the equilibrium concentration of the soluble oil fraction a is a cross-sectional area perpendicular to the flow direction q is the Darcy velocity within the contaminated zone As in Hoffmann's analysis, the source is assumed to be a point source of constant strength. The steady-state solution to the transport equation is examined to determine the maximum areal extent of a plume generated by contact with the oil lens, given some threshold concentration of environmental concern. Related investigations which model the spread of soluble hydrocarbon fractions under constant source strength conditions are given in (100,101).
There are some subtle differences between the specification of the source terms in the investigations described above. Hoffmann views the contact between the two phases as occurring along an abrupt and essentially horizontal interface. He acknowledges that the mass transfer between these two phases may be velocity dependent. Thus, as groundwater flow velocity varies, the source term may also vary. In the approach of Fried et al., however, the water is viewed as flowing through a residual oil zone within the capillary fringe. Again the source strength is a function of contact area, but here that area is perpendicular to the flow field. In addition, instantaneous equilibrium between the fluid phases is assumed. Thus, the source strength will be velocity independent. Neither approach allows for soluble component depletion as dissolution proceeds. Also, the nonaqueous zone must have reached its final distribution in the soil before this method of analysis is appropriate. It is clear that these approaches cannot be used to determine concentration profiles near the source, since both are predicated on the assumption of a point source. If the size and final configuration of the spill can be measured, however, estimation of the extent of contaminant migration via these methods appears feasible.
All of the transport models presented above focus on the aqueous transport of a soluble oil fraction. Organic depletion due to water phase partitioning is neglected in these models and source terms are assumed constant or simple functions of velocity. A more sophisticated model is presented in (102). Here the objective is still the modeling of the organic solute transport, but an attempt is made to account for the organic source depletion in the unsaturated zone above the water table. A one-dimensional water balance model (SESOIL), that was developed for the U.S. EPA, is coupled with a USGS two-dimensional groundwater flow code (SUTRA) to produce a new code called ULTRA. With this approach, monthly moisture content fluctuations in the vadose zone are modeled as a function of soil characteristics and precipitation, evaporation, and runoff data. In order to use this model, the initial distribution of the organic above the water table must be known or estimated, and this organic must be at its residual value (no organic flow). Equilibrium partitioning between the phases is assumed and biodegradation of organic in the vadose and saturated zones is incorporated. Zero-order kinetics for the biodegradation reactions is used at high concentrations, while a first-order reaction term is employed at low concentrations. The model accounts for dissolution and transport to the saturated zone by infiltrating rainfall and for contamination due to a rising water table. Volatilization from the water table is also modeled.
The authors employed their model to study two hydrocarbon spill sites in Florida. Data from the first site was used to calibrate the model, and then the calibrated model was run to simulate conditions at the second site. Limited data were available, but concentration trends similar to those observed in the field were predicted. The model severely overestimated the plume area. However, based on their study, the authors suggest that calibrated model concentration predictions for similar spill sites in Florida would be accurate to within a factor of five for actual concentration values.
Baehr and Corapcioglu (103,104) present a model which treats interphase partitioning in a more rigorous manner. Their model describes the one-dimensional transport of organic vapor as well as organic solutes within the unsaturated zone. Biotic and abiotic reactions are incorporated in the general model and a multicomponent composition of the organic phase is permitted.
Corapcioglu and Baehr (103) base the formulation on species mass balance in the form of Eq. (6). A Fickian nonadvective flux and modified Darcy's law expression are incorporated to obtain equations in the form of Eq. (15). Equilibrium partitioning of the organic is assumed with partition coefficients [Eq. (16)] as functions of constituent activity coefficients. Biological degradation is accounted for by assuming that it is oxygen-limited. Thus, a separate mass balance equation in the form of Eq. (15) is also written for oxygen. Adsorption is included but other abiotic reactions are neglected.
A one-dimensional numerical model is presented by Baehr and Corapcioglu (104) that solves the transport problem for the situation where the organic phase is at residual level (immobile) within the unsaturated zone. Constant water flux, constant water phase density, and uniform moisture content are assumed as well as a homogeneous, isothermal soil. Organic vapor transport is by diffusion only and diffusive transport within the organic phase is neglected. A uniform initial saturation and concentration distribution within the contaminated zone is assumed. To further simplify the analysis, the water phase is treated as an ideal solution and a linear adsorption isotherm is employed. Boundary conditions specified for the organic species include a zero concentration of contaminant at the soil surface and a zero diffusive flux at the water table. (Contaminants are assumed to enter the saturated groundwater zone only by advective transport of the water phase.) The resulting system of (weakly) nonlinear, parabolic equations is discretized within the model using a time-centered finite difference scheme. Nonlinear coefficients are evaluated by a method of forward projection and no iterations are performed. Solutions for each individual chemical component are computed sequentially. To evaluate the performance of the code, a mass balance check is incorporated in the model.
A related model is presented by Baehr (105). Here model applications are restricted to nonreactive, nondegradable contaminant phase components and the transport equation is rewritten to encompass a twodimensional (rz) system. As with the one-dimensional approach, the organic phase is assumed immobile, and the water flux is fixed at a constant rate in the vertical direction. Linear partition relations are also assumed. The resulting transport equation is very much like that presented in Eq. (18). It is a linear partial differential equation which accounts for vapor transport via diffusion and solute transport via dispersion, diffusion, and convection. One such equation can be written for each chemical species of interest. Silka (106) presents a similar cross-sectional vapor diffusion model that neglects water flux and considers only one organic component. The model has been employed to determine optimal grid spacing for soil-gas surveys in the field, when the objective is to detect surface contaminant sources at a distance.
Allan (107) also explores the movement of organic vapors within the unsaturated zone. His model is fundamentally an adapted finite-element contaminant transport code that solves an equation in the form of Eq. (18). Here the carrier fluid is gas rather than liquid water. Unlike the Baehr model discussed above, this model considers convective movement of the organic vapor as well as its diffusive transport. The constant pressure gradient driving force is input by the user, and the organic contaminant is introduced to the domain via the boundary conditions of the problem. As with the earlier analytical modeling efforts of Fried et al. (97), the source term is thus viewed as a constant. There is no depletion of contaminants and local equilibrium is assumed. The model accommodates neither organic phase nor water phase mobility and no reactions are considered.
One model has appeared in the literature that treats both the organic species transport problem and the immiscible phase flow of fluids. Many of the models discussed above can be viewed as a subset of this modeling approach. The equation development for this model is given by Abriola and Pinder (32). A one-dimensional numerical simulator is presented by Abriola and Pinder (108) and a two-dimensional (cross-sectional) extension of this approach is documented in (104,110,111). The model considers the twoand three-phase flow of a twocomponent organic liquid mixture in the unsaturated and saturated groundwater zones and permits the partitioning of one of these components into the gas and water phases. The other organic component is assumed nonreactive. Adsorption onto the soil matrix, chemical reaction, and biological degradation are neglected. As with the immiscible flow models discussed above, the gas phase is assumed to remain at atmospheric pressure. The model, however, is somewhat more general than many of these simulators, incorporating matrix and fluid phase compressibilities as well as composition-dependent organic phase density and viscosity. As in the organic transport models presented above, local equilibrium is assumed.
In the simulator, three species mass balance equations in the form of Eq. (17) are discretized by an implicit finite difference technique. The resulting set of nonlinear algebraic equations is solved simultaneously by a Newton-Raphson iteration procedure. Upstream weighting is incorporated as in the model of Faust (87). The model solves for two unknown capillary pressures and one mass fraction. Phase saturations and other mass fractions are obtained from saturation-pressure constitutive relations and equilibrium partition equations of the form (16). As in the immiscible flow models discussed above, hysteresis is neglected.
Convergence of the Abriola and Pinder model is demonstrated heuristically for a problem involving the onedimensional infiltration of a heavy oil and propane mixture into the unsaturated zone. Global mass balance results are also examined. Stability and convergence restrictions are reported by Pinder and Abriola (111) that limited the size of nodal spacing and time steps for the examined scenarios, making the simulations CPU intensive. Although such restrictions are not described in the papers of other investigators, it is anticipated that given the same input data and parameter relations, the immiscible flow models discussed above would also exhibit such behavior since they are based on similar equations.
The models discussed in this section on interphase transfer are quite varied. They range from the analytical description of solute transport in the aqueous phase to the numerical modeling of multi-species vapor and solute transport in the unsaturated zone. With the exception of the Abriola and Pinder model, all focus on some aspect of the species transport problem, neglecting the migration of the organic liquid phase and the nonsteady movement of water in the unsaturated zone. The weakness of such an approach is that it cannot quantify the complex evolution of the source term (the nonaqueous phase configuration) within the saturated or unsaturated zones nor can it treat such eventualities as a rising water table or nonsteady infiltration. Although the Abriola and Pinder model is capable of dealing with these scenarios, there is a large cost in computational time and model complexity. Using the Abriola and Pinder compositional approach, consideration of an additional chemical component involves the simultaneous solution of another mass balance equation, making large multispecies problems computationally intractable.

Model Validation
To date, few experiments have been documented in U.S. literature that are specifically designed to verify the mathematical models discussed in the previous section. A larger number of such experiments, however, were conducted in Europe during the late 1960s and early to mid-1970s and are, as yet, relatively unavailable to the American researcher in translation. These European experiments focused almost exclusively on the behavior of petroleum hydrocarbons in the subsurface and most centered on the validation of sharp interface models. Of the experiments conducted more recently in the United States on a wider variety of organic compounds, many have been unsuccessful or inconclusive in their attempt to validate models. In this section both groups of experiments are presented and discussed.
It should be noted that an additional and much larger group of experiments may also be assembled. This group would include all those experimental observations that are qualitative in nature and also those investigations that focus on the validity of some single assumption that is employed in one or more of the models. Although a thorough review of all such peripheral experimental work is beyond the scope of this work, many experiments having some bearing on the model assumptions in each category are presented below. Some observations related to field studies are highlighted in this regard. A number of experimental studies that focus on the measurement of the parameters for these models are also enumerated.

Sharp Interface Models
The ability of the piston-flow infiltration models to describe contaminant-infiltration episodes have been examined by a number of researchers. The most careful and complete measurements have apparently been conducted in Europe at the Technical University of Hannover (75,112,113); at the Eidgenossischen Technichen Hochschule (ETH), Zurich (114,115); and at the Bundesanstalt fur Gewasserkunde, Koblenz (22,(116)(117)(118)(119). One of the ETH reports has recently been translated into English for the U.S. Department of Energy (120). The other reports, however, are available only in their original German. Each of the experimental investigations cited above is described briefly below.
Mull (75) conducted a series of column experiments using a variety of sand mixtures for the porous matrix and two organic fluids: a fuel oil and benzene. Some of these investigations are summarized (in English) by Mull (76). Both horizontal and vertical infiltration studies were conducted for the constant head as well as the organic redistribution cases, as outlined in the previous chapter. Examination of Eqs. (24) and (25) reveals that the following parameters are required for front velocity predictions: hc,kro(so), 1, i0p',n,k. In Mull's experiments, porosities, intrinsic permeabilities, and organic phase viscosities and densities were measured directly. For each of the soils, moisture retention curves for a water-gas system were also determined. Based upon these retention curves, a capillary tube network model was employed to estimate h, and relative permeability functions.
Mull compared predicted and measured infiltration depths for water infiltration in unsaturated columns and organic phase infiltration. Frontal positions were determined by the use of monochromatic X-rays. In the water infiltration experiments, reasonable agreement was found between measured and computed frontal depths. For the experiments involving organic infiltration, frontal location predictions were made for a range of effective permeability values. These predicted infiltration curves were then shown to bracket the measured values. Mull concludes that predictions are reasonably accurate as long as infiltration distances exceed the mean capillary height hi. In Mull's experiments, fluid saturation distributions as well as frontal positions were determined within the columns. Saturation measurements during the organic redistribution experiments appear to validate the use of the assumed rectangular profiles shown in Figure 5.
Kleeberg (112) continued the work of Mull by examining the three-dimensional infiltration and lateral spreading along the water table of a fuel oil and benzene in sand tanks of various sizes. Capillary pressure saturation curves were measured for the organic-gas as well as the water-gas systems for one sand type. Organic-gas curves were also determined in the presence of water. It appears that relative permeabilities were neither predicted nor measured. Instead, representative curves were taken from the literature.
Kleeberg's report includes numerous plots of measured versus calculated frontal positions for vertical infiltration as well as lateral spreading. Data includes some saturation distributions as well as oil lens thickness measurements. Additional two-dimensional oil infiltration profiles are also given by Sprenger (113). The data presented in these reports reveal good agreement between measured and calculated values in most cases. Nearly all experiments were performed on a medium sand with mean grain size approximately 0.25 mm. For this material, then, it appears that the sharp front approximation works reasonably well. In the description of the experimental investigations, it is not clear, however, whether or not the relative permeability was adjusted to obtain the fits presented in the plots.
Schiegg (114,115) presents another series of one-and two-dimensional oil infiltration experiments. These are perhaps the most documented and detailed of the available laboratory investigations. Again, a relatively uniform medium quartz sand was used as the porous matrix. The selected organic phase was a mixture of boiling range benzene and a mercury-containing hydrocarbon compound. This second component was added to provide greater resolution for the gamma ray attenuation device which was employed to determine fluid saturations. In addition to saturation measurements, fluid pressures were determined using specially designed pressure transducer probes. Capillary pressure-saturation curves were measured for all two-phase combinations (organic-gas, water-gas, and water-organic). Apparently no attempt was made to estimate or measure relative permeabilities.
Although detailed measurements are presented for a number of different infiltration experiments along with precise information on boundary conditions, Schiegg does not compare his sharp interface model predictions to these measurements in his reports. Some comparisons are made, however, in a second series of experiments. Here, minimum oil lens thicknesses on the water  table were examined for different water table slopes. Schiegg did not find good agreement between his predictions based on a sharp interface model and these experimental measurements.
Results of one-and two-dimensional laboratory studies with a fuel oil (Heizol EL) as the infiltrating fluid are reported by Schwille (116)(117)(118). In comparison to the work of Schiegg discussed above, the observations are more qualitative in nature. Frontal positions were observed by the use of ultraviolet light and recorded with time-lapse photography. Layered soil systems were also considered. Intrinsic permeabilities and pore size distributions of the porous matrix materials were presented along with water-gas capillary curves. No comparisons with sharp interface modeling efforts were made in the paper. Some recent and limited experimental work in the United States has been directed towards verifying the sharp interface models. Illangasekare and Reible (80,81) report results of one-dimensional column studies. Organic infiltration and redistribution episodes were examined. Two different sand mixes (of mean diameters 0.3 mm and 0.4 mm) were employed in these experiments, as well as a range of organic compounds. Experimental results are presented for two immiscible organics: carbon tetrachloride and transmission oil, and two miscible compounds: chloroform and methanol. No attempt was made to measure either moisture retention or relative permeability curves. In these studies, the position of the wetting front was determined by visual observation. Columns contained a residual water saturation at the beginning of infiltration. The investigators compared model predictions of the frontal positions to experimental observations in a series of plots. Here, best fit predictions were obtained by calibration of three parameters: the effective conductivity and the effective capillary suction at the leading and trailing edges of the organic body. The authors argued that their sharp front model provides a good representation of the infiltration event, as is evidenced by the excellent agreement with the data. Many of their best-fit suction parameters, however, are negative. The authors do not comment on this puzzling result. The agreement achieved between model and measurements in this study is far superior to any reported in the European studies presented previously. Such a result is not surprising, however, since the data was fit with three adjustable parameters. It is clear that validation of the sharp interface model can only be achieved by independent measurement or es-timation of h, and kro as was done in some of the European studies. Illangasekare and Reible report a best fit value of kro = 1 in most of their experiments, another rather surprising result that would require further investigation.
Hochmuth and Sunada (85) report an experiment designed to verify their two-dimensional sharp-interface spreading model. The selected porous medium was comprised of 2.5-mm glass beads, and an 80% soltrol C-20% automatic transmission fluid mixture was employed as the organic fluid. Organic phase viscqsity and density were measured as well as saturated conductivity. As with the work of Illangasekare and Reible, average suction heads and corresponding values of relative permeability were not determined. Here, however, a range of reasonable values was examined. Residual saturations of organic were also varied. Visual observations were used to determine the shape and size of the oil lens and interface elevations were recorded as a function of time and distance. Difficulties were reported in estimating the position of the oil-air interface. No sharp interface was visible. Reasonable agreement, however, was found between model predictions and the observed position of the water-organic interface. A best-fit effective conductivity value was approximately one-eighth of the saturated value. The authors attribute this finding to capillary effects at the leading edge of the lens. A similar experiment was attempted with a finer matrix material (Ottawa sand) with no success. Capillary effects proved too significant in this case for the sharp front approximation to be reasonable.
From consideration of the experimental investigations described previously, it appears that the sharp interface approach to modeling organic infiltration in homogeneous media has been at least partially validated for a range of fuel oil products and a limited range of porous matrix materials. Its application to other types of organic fluids has yet to be demonstrated convincingly. It also appears critical that a larger range of porous matrix types be examined to determine the reasonable application limits of this approach.
Evidence for the validity of lateral spreading models is less conclusive and somewhat contradictory. Only the Kleeberg (112) study reports good agreement between model predictions and experimental observations. In a discussion of the multiphase flow problem, Schwille (121) notes that the size of the oil lens computed using sharp interface models has not been confirmed by field observations. He suggests that this may be due to water table fluctuations, which entrap more of the oil within the unsaturated zone. A recent field case study supports this view (122). It describes the importance of recharge and water table fluctuations to the distribution of hydrocarbons from a gasoline spill. Even though the number of experimental studies pertaining to lateral spreading is rather limited, it appears that this available data base has not yet been fully exploited. For example, the data presented in (112,(115)(116)(117)(118) could be compared to model predictions of Hochmuth and Sunada. The language barrier has apparently been an impediment to the careful analysis of the available data. It is also possible that quite a bit of the data from the experimental studies described above could be used to help in the validation of other types of immiscible flow models.

Immiscible Phase Flow Models
To date, few experimental studies have appeared in the literature that attempt to validate the use of blackoil type models in describing nonaqueous phase organic contamination events. This is primarily due to the difficulties inherent in the determination of the parameters required by such models and the measurement of model outputs such as fluid phase pressures and saturations. Because the models describe precise pressure and saturation distributions, visual observations of frontal positions are insufficient to validate model predictions. Laboratory techniques to determine pressure and saturation distributions are currently being refined (123)(124)(125)(126).
The two model validation experiments that have appeared in the literature are closely tied to validation of parametric models. In Kuppusamy et al. (89), results of a one-dimensional, two-phase flow study are compared to model predictions. Here, p-cymene (a derivative of benzene) was infiltrated into a water-saturated column. The porous matrix was a well-graded sand of known porosity and intrinsic permeability. For this study, capillary retention curves were measured for the two-phase system. Relative permeabilities were estimated using an extension of a capillary tube network model first proposed by Mualem (127) and presented by Parker et al. (128). In the column experiment, saturation distributions were not measured. Instead, water outflows were compared to model predictions. Visual inspection ofthe measured and computed outflow curves reveals a rather large discrepancy in curve shape as well as in outflow values. However, the authors suggest that confidence limits based on parameter uncertainty estimates will encompass the experimental data at moderate probability levels. Based on the data and measured parameters, it is not possible to determine the source of the apparent discrepancies. They could easily be due to an inappropriate model for relative permeability or to some problem inherent in the immiscible phase flow model mathematics, as well as to parameter uncertainty.
In a second study, results of a three-phase transient flow column experiment are presented by Lenhard et al. (123). Here, the porous matrix was again a wellgraded sandy material. The infiltrating liquid was an organic composed of a 1:9 mixture of 1-iodoheptane and soltrol-170. Two-phase capillary retention curves were measured for each fluid pair. As in the previous study, relative permeabilities were estimated from the capillary pressure data. A dual gamma apparatus was employed in conjunction with tensiometers to measure fluid phase saturations and pressures along the column. NaBr was added to the water phase to enhance the differences in the gamma attenuation coefficients of water and organic. A transient flow experiment was devised preserving a monotonically decreasing pathway for water and total liquid saturations. This was done to avoid complications due to hysteresis. The column was initially saturated with water, and a fixed volume of oil was permitted to infiltrate as the water table position was lowered. Approximately 41 min into the experiment, the water table was again lowered to permit greater air infiltration from above. This same scenario was modeled using the numerical simulator presented by Kuppusamy et al. (89). The authors report numerical difficulties in modeling the water table drop. To validate the model, observed saturation and pressure distributions at various depths were compared with computed values. Plots of these values reveal fairly good agreement, with the computed curve shapes mirroring the data. Large discrepancies in values were found at only one depth. These discrepancies were attributed by the authors to heterogeneities in column packing.
Although not directly tied to model validation, a number of papers have recently appeared in the literature involving the determination of model parameters. Much of the work focuses on representing capillary pressuresaturation relations and estimating relative permeabilities. The only actual measurements of relative permeabilities (outside of the oil industry literature) appear to be those reported for a water-TCE-gas system (111,129). Similar to the relations employed in many of the European experiments discussed in the previous section, Parker et al. (128) proposed a three-parameter model for capillary pressure curves based on a scaling procedure. Use of this procedure permits the separation of matrix from fluid properties. Scaling parameters may be estimated from the surface tension between the fluids. Relative permeabilities are then predicted from these capillary curves by employing a capillary tube model. Further documentation and verification of this scaling approach can be found in (92,130,131). Extensions of the model to incorporate hysteresis have appeared in (93,94).
The extension of two-phase capillary curves to a three-phase fluid system is based on the assumption that fluid wettabilities follow the order water > organic > air. Thus, the total liquid saturation is assumed to be a function of the gas-organic capillary pressure (Pg9), and water saturation is assumed to depend solely on organicwater pressure (Pow) In the formulation of Parker et al. (128), the residual saturation of any wetting fluid in a two-phase system is assumed to be a constant value.
There is some experimental evidence suggesting that the above assumptions may be over-simplified. Eckberg and Sunada (125) investigated the three-phase static distribution of water, air, and an organic phase (consisting of soltrol and automatic transmission fluid) in a porous medium composed of glass beads. Calculations of the total liquid content based on oil-air capillary pressures were found by these investigators to be inadequate. Lenhard et al. suggested that this finding of Eckberg and Sunada was due to lack of consideration of hysteresis effects (123). In their own study, however, they found large discrepancies between observed and predicted total liquid saturation relations.
Schiegg also reported large variability in the position and form of capillary curves encountered in his column studies (115). Here, however, discrepancies were found with the water-organic curves. This variation of curves was present despite the use of a very homogeneous porous medium. Schiegg suggests that some of this variation may be due to the occurrence of reversed wettability (organic > water > air). The possibility and importance of reversed wettability in groundwater systems has been considered by very few. Alharthi et al. explored the importance of wettability to organic phase retention and effective conductivity (126). They examined organic fluid spills in water-wet and completely dry soils. Large variations were found in organic retention and mobility based on this wetting history.
Observations of Schwille in reference to his experiments with organic fluids also have direct relevance to immiscible flow modeling (20,121). Based on his experience, Schwille noted that the concept of a single residual for a given pore-size distribution may be inappropriate. This observation is confirmed in a paper by Wilson and Conrad (132). In reviewing the oil industry literature pertaining to the measurement of residual saturations, they note that entrapped organic in twophase (water-organic) systems may greatly exceed residual levels in the vadose (three-phase) zone.
Of even greater significance are the experiments summarized by Schwille (22,119). Here, a series of one-and two-dimensional two-and three-phase organic chemical migration studies are reported. Schwille also investigated the behavior of a number of chlorinated solvents. Although his observations are primarily qualitative, these represent one of the few studies using dense organic compounds. In the displacement of water by organics, Schwille observed the frequent occurrence of fingering. This behavior occurred despite careful attempts to assure a homogeneous porous medium. In this situation, fingering was apparently caused by the displacement of a wetting fluid by a less viscous, nonwetting fluid. The implications of this observation are immense. The migration rate of an organic beneath the water table may be controlled by the propagation of such finger. There is currently no provision in the immiscible flow models to account for such behavior. Some references to this effect and attempts to model it have appeared in the oil industry literature (133)(134)(135). Application of such approaches to the propagation of organics in groundwater systems needs to be investigated.

Interphase Mass Transfer
Because of its complexity and the large number of parameters required by the compositional approach to nonaqueous phase pollutant modeling, there have been no documented attempts to validate such models experimentally. The investigations discussed in the previous section, however, will obviously have direct im-plications for the applicability of the immiscible flow portion of these codes. Studies relating to organic vapor transport as well as organic liquid-water and vapor partitioning have also been undertaken that have direct bearing on the appropriateness of the assumptions employed in the compositional models. These studies are summarized below. It is also apparent that any advances in our knowledge of solute transport or the biotic and abiotic transformation of organics in groundwater systems will have impact on developing compositional models. A review of the numerous investigations in this area, however, is beyond the scope of this paper.
One of the key assumptions inherent in the compositional modeling approach of Corapcioglu and Baehr (103) or Abriola and Pinder (32) is the assumption of local equilibrium between the fluid phases. This is also true for most of the simpler transport models discussed previously (79,105). Indeed, only Hoffmann considers the possibility of some mass transfer process that governs the strength of the source term (95,96).
Although the local equilibrium assumption implies that equilibrium concentrations should be achieved in the field in the presence of a nonaqueous phase, this review uncovered little field evidence to indicate that this actually occurs. In one study involving a field-scale experiment, concentrations were found to be close to equilibrium values beneath an infiltrated oil lens (136). It is reported, however, that at many nonaqueous, organic-phase contaminated field sites, maximum contaminant concentrations in groundwater are not observed to exceed a few percent of equilibrium levels (137). Schwille reports that the extent of oil contamination plumes are generally less than would be predicted, based on the local equilibrium assumption (20). Although field concentration data are available at many contaminated sites, it appears that no systematic study of this question has been made. It must be noted that there are many possible reasons (such as the effects of fingering, heterogeneities, biodegradation, or sample dilution) for the lack of field confirmation of the local equilibrium assumption.
At the laboratory scale, a few studies have been performed which support the validity of the local equilibrium assumption for organic liquid-water and vapor partitioning. They are relatively limited in scope, focusing on a very few contaminants and soil types. Van der Waarden et al. (138) presented the results of a series of column experiments using fine glass beads as the porous matrix. A mixture of hydrocarbons was injected into the unsaturated zone and permitted to redistribute. Water was then trickled into the top of the column and the concentration in the effluent was measured. The investigators concluded that the extraction of hydrocarbons could be described by an equilibrium partitioning process. Interpretation of data, however, is made difficult by the large variation in fluid saturations along the column as well as the presence of a gas phase.
In an organic liquid-vapor partitioning study, the recovery of petroleum hydrocarbons by forced air venting of soil columns was examined by Marley and Hoag (139).
Total hydrocarbon recovery rates in the outflow were measured and results for three different air flow rates were reported. The authors found good correlation between the recovery rates and a simple theoretical model, based on local equilibrium considerations. Analysis of the data is complicated by the presence of many different hydrocarbon components.
The local equilibrium assumption has also been examined in larger scale laboratory studies. In one investigation, crude oil was spilled in a large sand tank (100). Only sketchy details of the experiment are given, but it appears that the tank was irrigated with water and the leachate was collected for analysis. It is reported that concentrations in the effluent were "similar" to those obtained from measuring concentrations beneath a floating crude oil lens in a large beaker. In a second study, a residual zone of perchloroethylene was developed in a large sand tank within the saturated groundwater zone (140). Maximum concentration levels in the contaminant plume were examined for dependence on groundwater velocity. Measured concentration peaks corresponded to equilibrium levels.
Another series of column experiments is described by Fried et al. (97). These studies were conducted for two different hydrocarbons: a mixture of toluene and isooctane, and a gas oil. The soil type and characteristics are not specified. Hydrocarbons were injected into a watersaturated soil column to form a residual organic saturation. Water was subsequently forced at various rates through the column and organic concentrations in the effluent were measured. Although the authors report the attainment of equilibrium concentrations, their data seem to suggest the dependence of mass transfer on velocity (141). Similar experiments conducted by Hoffmann show no such velocity dependence (95).
Based on mass transfer considerations, one would anticipate that the validity of the local equilibrium assumption would depend on fluid velocity as well as the configuration of the entrapped residual. Pfannkuch attempts to make sense of the experimental data of Hoffmann, Fried et al., and others by examining it in this context (141). He suggests that one can interpret the available data by acknowledging that at low Peclet numbers the mass exchange process will be governed by diffusion, while for larger velocities, it will be solubility controlled. Similar conclusions have been reached by other investigators (142).
There are some other experimental studies which support the existence of such a velocity effect. In a preliminary study, column experiments to examine the partitioning behavior of pure alkanes are described (143). Although experimental details are not given, a measured correlation between average velocity and inferred mass transfer coefficient is reported for velocity ranges consistent with natural systems. A similar velocity dependence has been observed in soil-venting experiments (144). Here, however, partitioning is occurring between the organic and gas phases.
In light of the contradictory laboratory and field evidence summarized above, it appears that the validity of the local equilibrium assumption warrants further investigation. Exploration of the configuration of organic residuals in the pore spaces will also be important in this regard (132). Nonequilibrium effects could conceivably become of great significance to the modeling of remedial cleanup schemes such as aquifer flushing or soil venting. The nonequilibrium effects would also have direct implications for the natural spread, distribution, and persistence of organics in the subsurface environment.
The importance of vapor transport to the subsurface migration of organics has only recently been appreciated. Wilson et al. (145) presented the results of a laboratory investigation of the transport of organic pollutants in soil columns. By infiltrating organic solutions in unsaturated columns, they found that 30 to 90% of the organic mass was lost to volatilization during percolation. In a field evaluation of fuel oil spill, it was reported that hydrocarbons were present hydraulically up gradient from the source (146). This spread was attributed to the volatilization of organics. Other documented case studies (147,148) also support the importance of vapor transport to the subsurface distribution of contaminants.
Shallow soil gas sampling has become a means for detecting the presence of a contaminant plume in the deeper subsurface (149,151). Marrin and Thompson (152) detailed the relation of soil gas organic concentrations to solute concentrations in a contaminant plume at a field site. They suggested that water-gas partitioning and subsequent gaseous diffusion are the dominant transport mechanisms. Some contradictory observations, however, are given by Thomson (153). In a combination field/laboratory study, Thomson found anomalously high concentrations of organics in the vapor phase than could be explained by equilibrium partitioning considerations. He suggests that this anomaly could be due to lateral vapor movement from other areas of the subsurface or organic trapping by less permeable layers.
In his works, Schwille stresses the significance of vapor transport to the distribution of contaminants (20,21,121). He reported an experiment that demonstrates the impact of vapor phase density as well as diffusion on the migration of organics (22). The gas venting experiments of Thornton and Wootan (144) also support the existence of such a density effect.
With the exception of the vapor transport model presented by Allan (107), none ofthe existing compositional models incorporate convection or density effects. From considering the above studies, it appears that such effects could play an important role in the distribution of contaminants. The significance of these transport mechanisms warrants further investigation at the laboratory level as well as the field scale. Numerical modeling studies which incorporate gas convection would be useful to isolate the important driving mechanisms in various hypothetical contamination scenarios.
Because diffusion is a significant driving force for vapor transport (the sole force in most ofthe mathematical models) , the estimation of vapor diffusion coefficients becomes of great importance. While the free air diffusion coefficients of many organics of interest have been measured or estimated (154), the effect of the porous medium and fluid saturations on this coefficient must be taken into account. Several corrections to the free air diffusion coefficient have been proposed in the literature (155)(156)(157)(158). These corrections are functions of gas phase saturation and matrix porosity. Such corrections have been used to model the one-dimensional transport of pesticides in the unsaturated zone (159,160) and are employed in at least one vapor transport code (105).
Laboratory studies to confirm the use of such correction factors for the diffusive transport of volatile organics are reported in Thomson (153), Bruell and Hoag (161) and Thibodeux et al. (162). In the Thomson study, fair agreement was found between measured and estimated effective diffusion coefficients for several volatile organics. Here, a Penman correction was used. Bruell and Hoag conducted steady-state diffusion experiments in dry and wet sands using a variety of analytical grade gasoline hydrocarbons. They found a good fit between their data and that of the Millington-Quirk model. It was also demonstrated that the correction factor was independent of the temperature or hydrocarbon component. Sorption effects were eliminated from their analysis by considering only steady-state diffusion rates.
In their investigation of organic vapor diffusion in clay-sand mixtures, however, Thibodeaux et al. found that measured effective diffusion coefficients were significantly higher than would be predicted by correction expressions. They suggested that surface diffusion may account for this diffusion enhancement. Substantial sorption of vapors onto the soil matrix was also observed. Most existing compositional models neglect such sorption effects. Thus, it appears that additional laboratory studies are warranted to isolate the various factors contributing to diffusive vapor transport and to confirm the application of Millington-Quirk type correction factors to free air diffusivities.

Summary and Conclusions
This paper examined the fundamental equations employed to describe multiphase flow in porous media and presented existing MPFT models which have been developed to simulate the behavior of such systems. These models have been classified into three groups: (1) sharp interface models, (2) immiscible flow models with capillarity, and (3) interphase mass transfer models.
The major advantages of the sharp interface models are the relative simplicity of their mathematical formulation and their reduced parameter requirements over the other model groups. These factors make this class of models especially attractive for planning or screening purposes.
A number of laboratory studies were reviewed that attempted to validate the sharp interface models. It appears that the one-dimensional infiltration models can reproduce the infiltration behavior observed in the laboratory for a range of organics and a limited number of soil types (glass beads and medium sands). Although it is acknowledged that the piston-flow assumption will be inappropriate in finer grained and well-graded materials, limits for its application have yet to be determined. Whether or not such models could be used in a predictive capacity depends on our ability to independently measure or estimate the required model parameters. Unfortunately, many of the studies have been inconclusive in that the parameters were determined by calibration with the data. Use of these calibrated parameters for the description of infiltration behavior in similar systems or for the simulation of infiltration episodes subject to different boundary conditions has not been adequately examined. Validation efforts for lateral spreading models have met with comparatively less success. A number of investigators reported that model predictions did not mirror field or laboratory observations. Discrepancies appear to relate to the dominance of capillary effects in the lateral spreading process as well as to natural fluctuations of the water table. Such effects are not incorporated into the sharp interface models.
Based on the development of the sharp interface models, there are a number of limitations which govern the models' usefulness. The first is that these models are only appropriate for the characterization of spill episodes in homogeneous systems. Although in theory one could extend the sharp interface approach to analyze infiltration in layered soil systems, it is anticipated that gravity forces may not dominate in such systems. Thus, the sharp interface assumption might not be appropriate. A second limitation of these models is their inability to describe subsequent infiltration episodes or displacement events such as infiltrating rainfall or a rising water table. Similarly, the models cannot be used to explore the success of aquifer restoration schemes. Finally, the lateral spreading models are clearly restricted in application to organics having a specific gravity less than one.
The second category of models is undoubtedly more versatile than the sharp interface models. Models within this group differ very little from one another. Variations lie chiefly in their method of equation solution or their dimensionality. In these models, capillary curves are not idealized, and more than one phase can be mobile in the pore space at any time. The approach is not restricted to any particular class of immiscible organic compounds or matrix materials. Complicating factors such as hysteresis, although not presently incorporated in any of the existing models, can also be handled using this approach.
Certain conceptual assumptions in this modeling approach, however, appear to warrant further investigation. One assumption is that of a static air phase. The immiscible flow models assume that variations in gasphase pressures do not play a role in the migration of contaminants. This premise has yet to be adequately verified. A second assumption of perhaps even greater significance is that the infiltration of organic is a stable displacement process. Some experimental evidence has been reviewed that demonstrates the presence of fingering at the displacement front. The development and growth of such fingers could be the controlling mechanism for the infiltration of organics in many contamination scenarios. Existing immiscible phase flow models have no capability to describe this fingering phenomenon.
Because the immiscible phase flow models can only be as accurate as their parameter inputs, true model validation is necessarily linked to the verification of parameter estimation techniques. Certain assumptions inherent in existing parameter estimation models appear to require further study. These include the assumed wettability order of the fluids in the porous medium, the uniqueness of the wetting fluid residual saturation value, and the validity of employing capillary tube network models to predict relative permeability relations (especially two-phase, organic-water permeabilities and extensions to three-phase systems).
One obvious drawback to these models is their large data requirements. In addition, the immiscible flow models also have large computational requirements in terms of CPU time and storage. Saturation front resolution requires the refinement of nodal spacing, while the highly nonlinear character of the governing equations restricts time step size. Existing models employ simultaneous solution methods throughout the domain. Thus, the investigation and development of alternative solution schemes would make the use of these models more feasible in a field setting.
Although not discussed frequently in the literature, it appears that computational difficulties often arise in the simulation of organic contamination scenarios with immiscible flow models. Few simulation examples have actually been documented in the literature, and global mass balances are not reported for those examples that are presented. It is this reviewer's experience that most of these computational difficulties arise from the improper or inconsistent specification of boundary or initial conditions or the use of inappropriate or inconsistent parameter models for capillary pressure and relative permeability relations. Thus, research needs to be directed towards the formulation of guidelines for applying these models.
A final limitation of the immiscible flow modeling approach is that it neglects the chemical interaction of the immiscible organic phase with the other pore fluids and the matrix material. From the environmentalist's perspective, these interactions are of primary importance since they ultimately control the persistence and distribution of contaminants in the subsurface. It is only through an understanding of these interactions that one can hope to predict dose levels to the biosphere or to design effective remedial strategies.
The third group of models attempts to address these chemical interactions of the immiscible organic phase with the other phases. Models in this category are the most varied of all the model groups. Some focus on the unsaturated soil zone and seek to describe organic vapor diffusion and dissolution. Others focus on the organicwater interactions in the saturated zone and the transport of organic solutes. With the exception of the Abriola and Pinder model (32,108), none of the reviewed models incorporates nonaqueous organic phase transport. Instead, the organic is treated as an entrapped residual fluid of known configuration. This residual phase serves as a source term for organic component partitioning into the other phases.
Most of the models in this category are based on the solution of convection-dispersion equations of the type customarily employed to model solute transport in ground water systems. Model equations are frequently linear or weakly nonlinear. Thus, the computational burden of the models in this group is generally less than that of the other models. This observation would change, however, if more complex chemical reactions or nonequilibrium effects were incorporated into the codes.
Because of the difficulties inherent in quantifying the configuration of the organic residual source, there are few documented attempts to validate these models in the field. Some laboratory efforts, however, have been directed towards verifying the assumptions on which these models are based. Two important assumptions are those of equilibrium partitioning between the phases and the dominance of diffusion as the primary vapor transport mechanism. From an examination ofthe available data, it appears likely that there are instances in which each of these assumptions is not justified. Thus, additional experimental work needs to be directed towards the delineation of the limits of the applicability of these assumptions and to the development of alternative formulations.
As with the immiscible flow models, parameter estimation is of great importance to the validation of the interphase transfer models. Models dealing with the saturated and unsaturated groundwater zones require the specification of dispersion coefficients. The difficulties inherent in predicting these coefficients are well documented in the solute transport literature, and a great deal of research is presently being conducted in this area. Those models dealing with unsaturated zone transport also require the specification of vapor phase diffusion coefficients. To date, relatively less attention has been focused on this topic, especially in relation to determining such coefficients for volatile organic compounds.
Although much progress has obviously been made in recent years in the modeling of multiphase organic chemical flow and transport, it is clear that this field is still in its infancy. A rather large number of models have been presented in the literature, but none have been convincingly validated with laboratory data. Confident applications of such models to field sites for organic transport prediction, thus, appears to be reasonably far off. This statement is not to imply, however, that these models do not have value. It is through the use of such models in conjunction with laboratory data that we will be able to isolate the important areas for future study.
Even though there is a very limited data base for multiphase flow models, it appears that this data base is under-exploited. The careful experiments conducted in Europe in the 1970s have not yet been fully utilized for model validation. Similarly, available field measurements could be reviewed to help support or contradict model assumptions.
Sharp interface models may be of more immediate applicability to field investigations because of their relative simplicity. It must be cautioned, however, that these models are severely limited in their capabilities. In aquifer restoration operations, we are ultimately interested in questions that can only be answered by the application of compositional-type models. The longer term research challenge, then, is to validate and refine these models to make them more flexible, efficient, and accessible.