PDF995, Job 6 ______________ *Corresponding author. 1 Agricultural Sciences, 7(2):1-7 (2002) © 2002 Sultan Qaboos University Modeling Residual NAPL in Water-Wet Porous Media R.J. Lenhard1*, A.R. Kacimov2, A.M. Tartakovsky1, and H. AbdelRahman2 1Idaho National Engineering and Environmental Laboratory (INEEL), Idaho Falls, ID, USA, 83415-2025 2Department of Soil and Water Sciences, College of Agricultural and Marine Sciences, Sultan Qaboos University, PO Box 34, Al-Khod 123, Muscat, Sultanate of Oman ABSTRACT: A model is outlined that predicts NAPL which is held in pore wedges and as films or lenses on solid and water surfaces and contributes negligibly to NAPL advection. This is conceptually referred to as residual NAPL. Since residual NAPL is immobile, it remains in the vadose zone after all free NAPL has drained. Residual NAPL is very important because it is a long-term source for groundwater contamination. Recent laboratory experiments have demonstrated that current models for predicting subsurface NAPL behavior are inadequate because they do not correctly predict residual NAPL. The main reason for the failure is a deficiency in the current constitutive theories for multiphase flow that are used in numerical simulators. Multiphase constitutive theory governs the relations among relative permeability, saturation, and pressure for fluid systems (i.e., air, NAPL, water). In this paper, we outline a model describing relations between fluid saturations and pressures that can be combined with existing multiphase constitutive theory to predict residual NAPL. We test the revised constitutive theory by applying it to a scenario involving NAPL imbibition and drainage, as well as water imbibition and drainage. The results suggest that the revised constitutive theory is able to predict the distribution of residual NAPL in the vadose zone as a function of saturation-path history. The revised model describing relations between fluid saturation and pressures will help toward developing or improving numerical multiphase flow simulators. Keywords: Porous media, advection, groundwater contamination, multiphase flow, constitutive theory, vadose zone. !"#$%&!'()*+!,(-!.)&#/%&!0+!'(1234%&!'(45%&!'67489:;<=>!? #%&!!"! 4%&!!+*/! !"#$%&'()*+',#-*./012/3+43*+'540/%&-*62",748*+25&-*9+4:;%*<$=-+** ! !'##$%&*>********************+2?%&* ,!4"*$,#*$!&20%&*!"* ,%&('"%&* ,%"/%&*()2'%%*4%*+2"5* %)+2%&*,-%.*/%*23012345***<4"40"%4=* ?(3$"%&* ***************-.*6* =$'%&-*74"%&*8?79*:;)* ,<#9*-9* &,)+*<4&(%=* >%,.*:%;)*?2%/3*.%'%&-* ,;,@%7A&*B.403*C* ,&('"%&* ,"/%&*, ******** ,"40"%&*D=4;"%&*.1*$!&20%&*E$F3* ,;")*.1*G***********.1*:&(3*4:5H1*6*E$F'3*C*$!&20%&*,-.*!"* ,&('"%&* ,"/%&*?9*4"=- **************:%I*J#/5*')*E$F';%* ;=4&%&*<4,"/%&*?2/3*?9*'K=* ,&('"%&* ,"/%&*,-.*?2/'%3-*$%L2=*!%"/%* =$'%;%* ,"4%0"%&*D=4;%"%& *****D%")9*D=4;%"*G********* ,%&('"%&* ,%"/%&*,-%.*$(%'K30***4.$M%N*B%#+*1************,4,"%&*O2;'%*&+'P"*?2/3*?9*!/"L*4%:5Q*R$,%(8* ,%".9*<&+ *****',K(%&*S'"%&*:;)* ,12T%&*G*** ,148*J0,%*$!&20%&*,-.*E2;7*()2'%* ,%4F%&*U+4";%&*?9*R$,VQ&* ,;"K"%&*W+4T'%&*<$:X9 ****Y(%;'3*C-=*8,FP%&*$/<%4=*E2;0%&*&-:*G0,!$%&*Z(0%&*?I* L$.2T%&*<4L$[;%&*$<1*.1*!"/L*U+4";%&*,-.*$<1*7&+-*. ************* L9'K%&*<484F"%&*.1*$"K'0L*\-%&-*+&2=Q&*9'K'%"%&*D1''%;%* ,%%4F%&*G******.'%&*..*+&2=Q&*R9'K'"%&* L$.2T%&* L$[;%&*?I ****$!&20%&*]M^-*((<'%&-* ,(0;%&* L+4@;%&*!,=*<4%)_K%&*B%/F30,!4"*$,#*$!&20%&-*6*74"*6*7&2.*`*\9+2?%&* *G*1*,-.*ab;3 *********** ,%4F%&*+&2=Q&*R9'K'"* L$.2T%&* L$[;%&*("*c*4"9&*!/"L*.'%&-*$!&20%&*]M^-*((<'%&*!,=*<4)_K%&*dPL*4*+2"5* )+2%& ************$!&20%&*,-.*!%"*.%&('"%&*e4%P"9&-*dL$%P3-*()2'%%*6**G********()23*:;)*R+94)* ,;L2/'%&* L$.2T%&* L$[;%&*?Y=*f!4';%&*ab;3- ***,-.*e4P"9&*(%Lg23**********((<'%&*+40"*:;)* %&'8*74"%&*:%I* 14^A4=*$!&20%&**G******\-%&-*(*&$"%&*U+2";%&*(T<3* 7&+'%&*,-. **************************** L9'K%&*<484F"%&*$L2?3*:;)*')4%0,7*,+-'%=*\-%%&-*$!&2%0%&*]M%^-*((%<'%&* &?;%"*.%1*W$'%%&*!%,=*<4%)_K%&*d%PL *+&2=Q&*9'K'"*D1';%* ;0F"%&G* LENHARD et al. 2 he contamination of groundwater by fuel and industrial chemicals occurs throughout the world. A major class of contaminants is nonaqueous phase liquids (NAPLs). The risk that they pose to public health and the environment can be severe because many of these compounds are known or suspected carcinogens. Several cases of NAPL migration to groundwater have been reported in the Sultanate of Oman. A leak at a British Petroleum gasoline station in Al Kamil contaminated the vadose zone as well as the groundwater. Discovery of petrol compounds in groundwater wells located several tens of meters from the leak source initiated a series of actions, including construction of an air-sparging system (MWR, 2000). Contamination of groundwater close to major roadways is another environmental concern in the Sultanate of Oman. During dry periods, NAPLs (principally petroleum products) become emplaced along roadways because of improper disposal and leaks from automobiles. The NAPL will occupy soil-pore wedges and small pores, because the water content in sandy soils is very low. The NAPL will be held in the vadose zone against gravity because it is immobile. During brief rainy periods, which include torrential rain, compounds from the NAPL will dissolve and move downward with the water to the groundwater, effectively delivering a pulse of contaminants to the groundwater. The migration of NAPLs toward groundwater is a complex process. In addition to the physical movement of NAPL, volatile and soluble components of the NAPLs can move toward groundwater in the gas and aqueous phases. Once in contact with groundwater, these components can partition into groundwater. NAPLs less dense than water (LNAPLs) will accumulate above the water table and may depress the water-saturated region. NAPLs more dense than water (DNAPLs) can penetrate below the water table. Fuel compounds, such as hydrocarbons, are LNAPLs and industrial solvents, such as trichloroethylene and carbon tetrachloride, are DNAPLs (Fetter, 1993). As NAPL migrates downward through the unsaturated zone (vadose zone) toward groundwater, some of the NAPL remains in the vadose zone. That portion which remains likely resides in pore wedges and as thin films or lenses on water surfaces. As such, this NAPL is relatively immobile (or negligibly mobile). The negligibly mobile NAPL that does not drain from the pore spaces in the vadose zone is referred to as residual NAPL and can be continuous or discontinuous throughout the pore spaces. Residual NAPL can serve as a long-term source for groundwater contamination. Because of its persistent nature, residual NAPL needs to be considered when developing plans for site characterization and cleanup strategies and technologies. Modeling has been recognized as a scientific method to help cost-effectively manage subsurface contamination. To predict multiphase flow, a set of nonlinear partial differential equations needs to be solved. A solution can be obtained for given boundary and initial conditions only if the constitutive relations are known. For multiphase flow, constitutive relations are the relationships among relative permeability (k), saturation (S), and pressure (P), i.e., k-S-P relations. A deficiency in current constitutive theory is its inability to predict residual NAPL (Oostrom et al., 2003). The commonly used constitutive theory predicts that NAPL will completely drain from the vadose zone over time, which is contrary to laboratory and field observations. In a controlled experiment by Hofstee et al. (1998), a residual NAPL saturation ranging to 20% was measured, but computer simulations failed to accurately predict the behavior of the residual NAPL. In another investigation by Oostrom and Lenhard (2002) that included sand and porous media high in calcium carbonate, numerical simulations using common constitutive theory were unable to satisfactorily predict the measured residual NAPL saturation in the experiment. Improvements in constitutive theory are needed to accurately predict the subsurface behavior of NAPLs in the vadose zone. In this paper, we will describe improvements to commonly used constitutive theory that allow residual NAPL to be predicted using numerical simulators. The improvements are based on the work of Lenhard et al. (2003). Also in this paper, we will apply the revised constitutive theory to a scenario involving NAPL imbibition and drainage and water imbibtion and drainage. Unfortunately, we will not be able to present the detailed revised constitutive theory proposed by Lenhard et al. (2003) or discuss other recently proposed modifications to current constitutive theory (Wipfler and van der Zee, 2001; van Geel and Roy, 2002) due to space constraints. Model Description Probably, the most commonly used constitutive theory for air-NAPL-water flow is that published by Parker et al. (1987), Parker and Lenhard (1987), and Lenhard and Parker (1987). The work by Parker et al. (1987) pertains to nonhysteretic k-S-P relations. The works by Parker and Lenhard (1987) and Lenhard and Parker (1987) pertain to hysteretic k-S-P relations. In the latter work, the concept of an apparent saturation is introduced and entrapment of NAPL and air are modeled. Both models are for strongly water-wet T MODELING RESIDUAL NAPL IN WATER-WET POROUS MEDIA 3 porous media where water is the wetting fluid and NAPL and air are nonwetting fluids, in which the NAPL is more wetting than air. The improvements presented in this paper are modifications to the hysteretic k-S-P relations. The key modifications are 1) redefining the total-NAPL saturation, 2) proposing an equation to predict the residual-NAPL saturation as a function of saturation-path history, and 3) calculating entrapped air in the residual NAPL. Major aspects of these modifications will be presented below as well as major components of the hysteretic k-S-P relations. For a more detailed description of the modifications, the readers are referred to Lenhard et al. (2003). Before presenting the model, we want to inform the readers about the terminology that we will employ to describe fluid saturations. We define trapped fluid as that which becomes discontinuous and occluded by more wetting fluids. Since trapped fluids are discontinuous, they are immobile in sandy or finer porous media at normal fluid flow conditions. We define residual fluid as that which is either immobile or negligibly mobile and not occluded by another fluid. We believe residual fluid can be either discontinuous or continuous. Residual fluid that is contained in pore wedges or in lenses on water surfaces (nonspreading liquids) may be discontinuous or continuous, depending on the direction of fluid movement and the fluid saturation. Residual fluid that is a thin film on water surfaces is likely to be continuous (spreading liquids). We believe that discontinuous residual fluid is immobile and continuous residual fluid is negligibly mobile. This definition differs than that is typically used in the petroleum industry with regard to residual NAPL. What we call trapped NAPL, petroleum engineers call residual NAPL. To avoid confusion between residual water and residual NAPL, we have adapted the above definitions. We define free fluid as that which is continuous and mobile. S-P Model The apparent water saturation is defined as: )1( rw rwatwotw w S1 SSSS S − −++ = where Sw is the water saturation, Sot is the trapped-NAPL saturation occluded by the aqueous phase, Satw is the trapped-air saturation occluded by aqueous phase, and Srw is the residual-water saturation, which is assumed to be negligibly mobile and the value of which is independent of the water saturation-path history. Satw also includes the trapped air in the trapped NAPL that is occluded by the aqueous phase. The apparent total-liquid saturation is defined as: rw rwatow t S1 SSSS S − −++ = )2( where So is the total-NAPL saturation and Sat is the total trapped-air saturation. Note that Sat is the sum of trapped-air saturations occluded by the aqueous phase and the NAPL: atoatwat SSS += (3) where Sato is the trapped-air saturation occluded by the NAPL. The terms Sot, Satw, and Sat in equations (1) and (2) are functions of the water and total-liquid saturation-path histories [see Lenhard (1992) for descriptions]. The first modification to the constitutive theory involves redefining the total-NAPL saturation. In the hysteretic S-P model described by Parker and Lenhard (1987) and Lenhard (1992), the total-NAPL saturation consisted of two components: free NAPL and trapped NAPL. The modification is defining the total-NAPL saturation to consist of three components: free NAPL, trapped NAPL, and residual NAPL. Correspondingly, the total-NAPL saturation (So) is defined as: (4)rootofo SSSS ++= where Sof is the free-NAPL saturation and Sro is the residual-NAPL saturation. To determine So, each component needs to be calculated based on the total- liquid and water saturation-path histories and the NAPL-water and air-NAPL capillary pressures. Note that So does not approach zero as Sof approaches zero. The second modification involves developing an equation to predict the residual-NAPL saturation as a function of saturation-path history. Based on arguments given in Lenhard et al. (2003), the predictive equation is: )5(5.1w 5.0 w max tromrot )S1()SS(SS −−= )7( )6( rw rom rom rw atrro rot S1 S S S1 SS Swhere − = − + = in which rotS is the apparent residual-NAPL saturation, romS is the effective maximum residual-NAPL saturation, maxtS is the historic maximum apparent total-liquid saturation, and Satr is the trapped-air saturation in the residual NAPL. rotS is composed of LENHARD et al. 4 Sro and Satr. An effective nonwetting fluid saturation, such as the effective maximum residual-NAPL saturation, is determined accordingly: )8( rw x x S1 S S − = where Sx is the nonwetting fluid saturation. The third modification involves modeling entrapped air in the residual NAPL. Because the trapped air in the residual NAPL is occluded only by the NAPL (i.e., it will never also be occluded by water such as trapped air in trapped NAPL), it is a component of Sato [see equation (3)]. Therefore, Sato consists of two components: trapped air in free NAPL (Satof) and trapped air in residual NAPL (Satr): )9(atratofato SSS += Furthermore, the trapped air in the residual NAPL consists of two components, depending on the saturation-path history, which are defined as: )10(atroatrwatr SSS += where Satrw is the trapped-air saturation in the residual NAPL that resulted from advancing air-water interfaces and Satro is the trapped-air saturation in the residual NAPL that resulted from advancing air-NAPL interfaces. Both Satrw and Satro and their corresponding scaled values (effective saturations) are determined from the saturation-path history. To calculate Satrw and Satro, the smallest pores where air-water and air-NAPL interfaces occupied relative to those indexed by rotw SS + need to be considered. The smallest pores that possessed air-water interfaces are indexed by ∆ S w aw , and the smallest pores that possessed air-NAPL interfaces are indexed by minS . For the saturation paths in which rotw SS + is greater than aw wS∆ , air trapped in residual NAPL may have resulted from either advancing air-water or air-NAPL interfaces. The condition where mintS is larger than or equal to rotw SS + means that all trapped air in the residual NAPL resulted only from advancing air-water interfaces. )12(0S )11( S1 S,SmaxSS SS SSSS SS atro aw w aw wwrotw atrwatrw rotw min t aw wrotw andfor = − −+ = +≥>+ ! " ! # $ ! % ! & ' ( ) *+ , - ∆ ∆ ∆ The condition where mintS is smaller than rotw SS + means that the air trapped in the residual NAPL may have resulted from either advancing air- water or air-NAPL interfaces. The amount of air trapped in the residual NAPL by advancing air-NAPL interfaces can be determined from )13( S1 SSmaxSS SS SSS andS SSfor min t min t,wrotw aroatro rotw min t aw wrotw ! " ! # $ ! % ! & ' − ( ) * + , -−+ = +<>+ ∆ To determine the amount of air trapped in the residual NAPL by advancing air-water interfaces for this saturation path, mintS is compared to .S aw w ∆ )14( S1 S,SmaxS,Smax SS SS aw w aw www min t arwatrw aw w min t ! ! " !! # $ ! ! % !! & ' ( ) *+ , - ∆ ∆ ∆ − (( ) * ++ , - − = > )15(0S SS atrw aw w min t = ≤ ∆ For the saturation paths where rotw SS + is less than or equal to awwS∆ , any air trapped in residual NAPL resulted only from advancing air-NAPL interfaces. Depending on the value of mintS , there may be no air trapped in the residual NAPL. )17( )16( 0S 0S SSSandSSSfor atro atrw rotw min t aw wrotw = = +≥≤+ ∆ MODELING RESIDUAL NAPL IN WATER-WET POROUS MEDIA 5 )19( S1 S,SmaxSS SS )18(0S SSSandSSSfor min t w min trotw aroatro atrw rotw min t aw wrotw ! " ! # $ ! % ! & ' − ( ) *+ , -−+ = = +<≤+ ∆ In Equations (11)-(19), aroS and arwS are the maximum effective trapped-air saturations for the given saturation-path history, awwS∆ is the effective water saturation at the drying-to-wetting reversal from the main drainage S-P path, and mintS is the historic minimum apparent total-liquid saturation. aroS is the maximum effective trapped-air saturation that results from advancing air-NAPL interfaces, and arwS is the maximum effective trapped-air saturation that results from advancing air-water interfaces. All air trapped in a three-phase fluid system is a function of aroS , and all air trapped prior to initiation of a three-phase fluid system is a function of arwS . Readers are referred to Parker and Lenhard (1987) and Lenhard (1992) for descriptions of how aroS and arwS are determined. Discussion A computer program was written to calculate apparent, actual, trapped, and residual saturations for inputted capillary pressures, incorporating the above S- P modifications. A saturation-path scenario involving NAPL and water imbibition and drainage was used to test the revised S-P constitutive model. Table 1 lists the parameters that we used to generate the fluid saturations as functions of the capillary pressures. The parameters, which are all defined in Parker and Lenhard (1987) and Lenhard and Parker (1987) except Srom, reflect a sandy porous medium. The hypothetical saturation-path scenario begins with water drainage, yielding a two-phase air-water fluid system. Water drains until an air-water capillary pressure of 80 cm and a water saturation of 0.356. All fluid saturations unless otherwise noted are actual saturations (not scaled). The next condition, which we call Path-2, involves NAPL imbibition. The total-liquid saturation increased to 0.751. We simulated a slight decrease in the water saturation (i.e., 0.005 saturation units) as the NAPL front passed. During NAPL imbibition, air will become trapped in NAPL by advancing air-NAPL interfaces. The residual-NAPL saturation will increase as maxtS increases, and there is no entrapped NAPL because water is on the main drainage path. Path- 3 involves NAPL drainage, while the water saturation TABLE 1 Parameters used to calculate saturations corresponding to the hypothetical saturation-path scenario. Parameter Value Parameter Value Parameter Value dα 0.05 cm-1 aoβ 1.80 aro i S 0.20 iα 0.10 cm-1 owβ 2.25 arw i S 0.25 n 2.0 romS 0.20 or i S 0.25 rwS 0.15 increased only slightly. A rebound was simulated in the NAPL-water capillary pressure as the NAPL front passed (NAPL drained), which caused the water saturation to increase from 0.351 to 0.353 – thus causing very minor NAPL entrapment in water. The residual-NAPL saturation remains constant until the water saturation increases. Inspection of Equation (5) shows that as long as maxtS and wS are constant, rotS will not change. maxtS obtained its highest value during NAPL imbibition (i.e., maxtS = 0.751). Air that was trapped during NAPL imbibition is now being released as NAPL drains from the larger pore spaces, except the trapped air which remains in the residual NAPL. Path- 4 involves water imbibition (e.g., water table elevation increased). The NAPL saturation was held near 0.136, while the water saturation increased from 0.353 to 0.610. During water imbibition, free and residual NAPL will become occluded by water (trapped NAPL). The entrapped air saturation will also increase as the air-NAPL interfaces move into larger pore spaces. In the last saturation path (Path-5), the total- liquid saturation continually decreases. Early in the saturation path, water drains at a constant NAPL saturation of 0.136, and later in the saturation path, NAPL drains at a constant water saturation of 0.315. In the final stage of Path-5, all of the NAPL became residual NAPL. Results are shown in Figures 1 and 2. Figure 1 shows the saturations for Path-2 and Path-3. Figure 2 shows the saturations for Path-4 and Path-5. Note that Figure 1B is a continuation of Figure 1A, Figure 2A is a continuation of Figure 1B, and Figure 2B is a continuation of Figure 2A. The changes in saturations are shown as a function of the air-NAPL capillary pressure. By inspecting Figure 1 for Path-2, it can be seen that as the total-liquid saturation increases, both residual and free NAPL increase. When the difference between the total-liquid and water saturations is small, the relative amount of residual NAPL to the total NAPL is greater than for free NAPL; however, there is some free NAPL even at low total-NAPL saturations. LENHARD et al. 6 Figure 1. Total-liquid, water, residual-NAPL, free-NAPL, and trapped-NAPL saturations as a function of the air – NAPL capillary pressure for A) Path-2 and B) Path-3. Free NAPL first appeared when the total-NAPL saturation was only 0.015. The total-NAPL saturation at which free NAPL appears will be a function of the apparent water saturation. As the total-NAPL saturation increases, the relative amount of free NAPL increases. For Path-3, the total-NAPL saturation decreases while the water saturation is relatively constant. The decrease in the free-NAPL saturation corresponds to the decrease in the total- NAPL saturation, because the predicted residual-NAPL saturation remains constant. Further investigations are warranted to determine if the residual-NAPL saturation should remain constant for the saturation-path conditions of Path-3. Since it is very difficult to measure residual- and free-NAPL saturations simultaneously, the best approach to investigate this issue Figure 2. Total-liquid, water, residual-NAPL, free-NAPL, and trapped-NAPL saturations as a function of the air-NAPL capillary pressure for A) Path-4 and B) Path-5. is to compare numerical simulations to transient air- NAPL-water flow experiments. In Figure 2 for Path-4, the water saturation increased while the total-NAPL saturation was constant. This caused the total-liquid saturation to increase. As the water moved into larger and larger pores, some of the residual and free NAPL became trapped NAPL in water. The free-NAPL saturation decreased as the water saturation increased. The residual NAPL was not mobilized; it became trapped in water along with some free NAPL. Eventually, all of the NAPL became entrapped (the residual- and free- NAPL saturations are zero). When the apparent total- liquid saturation was equal to the apparent water saturation (all NAPL is entrapped), water drainage from the porous medium was simulated. As water MODELING RESIDUAL NAPL IN WATER-WET POROUS MEDIA 7 drained, trapped NAPL was released and became residual and free NAPL. Initially, most of the trapped NAPL that was released became residual NAPL. Shortly thereafter, most of the released NAPL became free NAPL. The amount of residual NAPL relative to free NAPL is different depending on whether water is imbibing (Figure 2A) or draining (Figure 2B). The reason for this is because the saturation-path history is different for Path-4 conditions and Path-5 conditions. The water saturation decreased until all of the trapped NAPL was released. At this stage, the residual-NAPL saturation is predicted to be greater than the free-NAPL saturation. During the later periods of Path-5, the water saturation was held constant while the free NAPL was allowed to drain. Eventually, the free- NAPL saturation became zero (all of the NAPL was residual NAPL). When air-NAPL capillary pressure was further increased, the model would not allow further drainage of NAPL and the total-liquid and water saturations remained constant. Comparing Figures 1 and 2, one can observe that the calculations of residual-, trapped-, and free-NAPL saturations for Path-4 and Path-5 yield very different results compared to Path-2 and Path-3. The predicted trapped-, residual-, and free-NAPL saturations can have complex patterns even though the total-liquid and water saturation increase or decrease in a smooth, monotonic manner. Conclusions The revised S-P constitutive relations proposed by Lenhard et al. (2003) were tested. A hypothetical saturation-path scenario was developed from which apparent, actual, trapped, free, and residual saturations were predicted. Analyses of the results suggest that the revised constitutive theory is able to predict the distribution of residual NAPL in the vadose zone as a function of saturation-path history. The revised model describing relations between fluid saturations and pressures will help toward developing or improving numerical multiphase flow simulators. However, further investigations are warranted to fully evaluate the constitutive relations. In particular, the revised constitutive relations need to be incorporated in numerical simulators and tested against transient air- NAPL-water flow experiments. Acknowledgment Appreciation is extended to the Laboratory Directed Research and Development program at the Idaho National Engineering and Environmental Laboratory (INEEL) for funding this work. The INEEL is operated for the US Department of Energy (DOE) by Bechtel BWXT Idaho, LLC under DOE Idaho Operations Office Contract DE-AC07- 99ID13727. References Fetter, C.W. 1993. Contaminant Hydrogeology. Macmillan Publishing Company, New York. 458pp. Hofstee, C., R.C. Walker, and J.H. Dane. 1998. Infiltration and redistribution of perchloroethylene in partially saturated, stratified porous media. J. Contam. Hydrol. 34:293-313. Lenhard, R.J. 1992. Measurement and modeling of three-phase saturation-pressure hysteresis. J. Contam. Hydrol. 9:243-269. Lenhard, R.J., M. Oostrom, and J.H. Dane. 2003. A constitutive model for air-NAPL-water flow in the vadose zone accounting for residual NAPL in strongly water-wet porous media. Submitted for publication J. Contam. Hydrol. Lenhard, R.J. and J.C. Parker. 1987. A model for hysteretic constitutive relations governing multiphase flow, 2. Permeability-saturation relations. Water Resour. Res. 23:2197- 2206. Ministry of Water Resources (MWR). 2000. Report on the Hydrocarbon Contamination in Willayat Al-Kamil. Report number N756. Ministry of Water Resources, Sultanate of Oman. Oostrom, M. and R.J. Lenhard. 2002. Flow behavior in unsaturated Hanford caliche material: An investigation of residual NAPL. Submitted for publication Vadose Zone J. Oostrom, M., C. Hofstee, R.J. Lenhard, and T.W. Wietsma. 2003. Flow behavior and residual saturation formation of liquid carbon tetrachloride in unsaturated heterogeneous porous media. J. Contam. Hydrol. (in press). Parker, J.C. and R.J. Lenhard. 1987. A model for hysteretic constitutive relations governing multiphase flow, 1. Saturation- pressure relations. Water Resour. Res. 23:2187-2196. Parker, J.C., R.J. Lenhard, and T. Kuppusamy. 1987. A parametric model for constitutive properties governing multiphase flow in porous media. Water Resour. Res. 23:618-624. Van Geel, P.J. and S.D. Roy. 2002. A proposed model to include a residual NAPL saturation in a hysteretic capillary pressure- saturation relationship. J. Contam. Hydrol. 58:79-110. Wipfler, E.L. and S.E.A.T.M. van der Zee. 2001. A set of constitutive relationships accounting for residual NAPL in the unsaturated zone. J. Contam. Hydrol. 50:53-77. ________________________________________ Received September 2002. Accepted November 2002.