VZJ Journal of Natural Resources and Life Sciences Education
HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS
QUICK SEARCH:   [advanced]


     

This Article
Right arrow Abstract Freely available
Right arrow Figures Only
Right arrow Full Text (PDF) Free
Right arrow Alert me when this article is cited
Right arrow Alert me if a correction is posted
Services
Right arrow Similar articles in this journal
Right arrow Similar articles in ISI Web of Science
Right arrow Alert me to new issues of the journal
Right arrow Download to citation manager
Right arrow reprints & permissions
Citing Articles
Right arrow Citing Articles via HighWire
Right arrow Citing Articles via ISI Web of Science (3)
Right arrow Citing Articles via Google Scholar
Google Scholar
Right arrow Articles by Oostrom, M.
Right arrow Articles by Wietsma, T. W.
Right arrow Search for Related Content
PubMed
Right arrow Articles by Oostrom, M.
Right arrow Articles by Wietsma, T. W.
GeoRef
Right arrow GeoRef Citation
Agricola
Right arrow Articles by Oostrom, M.
Right arrow Articles by Wietsma, T. W.
Related Collections
Right arrow Vadose Zone Processes and Chemical Transport
Right arrow Multiphase Fluid Flow
Published in Vadose Zone Journal 4:163-174 (2005)
© 2005 Soil Science Society of America
677 S. Segoe Rd., Madison, WI 53711 USA

ORIGINAL RESEARCH

A Comparison of Models Describing Residual NAPL Formation in the Vadose Zone

M. Oostroma,*, M. D. Whitea, R. J. Lenhardb, P. J. Van Geelc and T. W. Wietsmad

a Hydrology Group, Pacific Northwest National Laboratory, P.O. Box 999 MS K9-33, Richland, WA 99352
b Subsurface Science Initiative, Idaho National Environmental and Engineering Laboratory, P.O. Box 1625, Idaho Falls, ID 83415-2025
c Department of Civil and Environmental Engineering, Carleton University, 1125 Colonel By Drive, Ottawa, ON, Canada K1S 5B6
d EMSL, Pacific Northwest National Laboratory, P.O. Box 999 MS K8-96, Richland, WA 99352
* Corresponding author (mart.oostrom{at}pnl.gov)

Received 5 February 2004.



   ABSTRACT
TOP
 ABSTRACT
INTRODUCTION
 THE THEORETICAL MODELS
 COMPARISON WITH STATIC PRESSURE...
 COMPARISON WITH TRANSIENT COLUMN...
 SUMMARY
 REFERENCES
 
A major shortcoming of multifluid flow simulators is the inability to predict the retention of nonaqueous phase liquid (NAPL) in the vadose zone after long drainage periods. Recently, three theoretical models—the Wipfler and van der Zee (WVDZ) model; the Van Geel and Roy (VGR) model; and the Lenhard, Oostrom, and Dane (LOD) model—have been proposed for describing residual NAPL formation. The WVDZ model assumes a critical total liquid saturation below which all NAPL becomes residual. The VGR and LOD models are extensions of an existing hysteretic relative permeability–saturation–capillary pressure model and assume formation of residual NAPL during NAPL drainage and imbibition, respectively. In this study, we compared model predictions against results of a series of static pressure cell experiments. We found no experimental evidence supporting the WVDZ concept of a critical total liquid saturation. The other two models yielded reasonable predictions. The VGR and LOD models were then incorporated into a multifluid flow simulator, and simulations of two transient column experiments were conducted. Both models performed considerably better than simulations without considering the formation of residual NAPL, underwriting the importance of incorporating this process in simulators. Although the VGR and LOD models are based on different conceptual models, no clear performance differences could be observed when simulation results were compared against the transient experimental data.

Abbreviations: LOD, Lenhard, Oostrom, and Dane • NAPL, nonaqueous phase liquid • NH, nonhysteretic Brooks and Corey–Burdine constitutive model • PCE, perchloroethylene • VGR, Van Geel and Roy • WVDZ, Wipfler and van der Zee


   INTRODUCTION
TOP
 ABSTRACT
 INTRODUCTION
THE THEORETICAL MODELS
 COMPARISON WITH STATIC PRESSURE...
 COMPARISON WITH TRANSIENT COLUMN...
 SUMMARY
 REFERENCES
 
DESCRIPTIONS OF MOBILE AND ENTRAPPED NAPL behavior are now standard features of some multifluid flow simulators. However, despite recent contributions describing NAPL residual saturation formation (Lenhard et al., 2004; Van Geel and Roy, 2002; Wipfler and van der Zee, 2001), these models do not generally incorporate retention of NAPL in the vadose zone following NAPL imbibition events (e.g., surface spills, tank leaks). Commonly used constitutive relations assume a nonzero NAPL relative permeability when the nonentrapped NAPL saturation is greater than zero. As a result, this NAPL is allowed to drain from the pore spaces. Neglecting residual NAPL saturation formation might cause an overestimation of the volume of NAPL that reaches the groundwater (Van Geel and Roy, 2002). Because NAPL retained in the vadose zone can serve as a long-term source for groundwater contamination via transport through the gas or aqueous phase, understanding and predicting the processes of residual NAPL formation and remobilization are critical when considering restoration or management of a contaminated subsurface site.

Residual saturation is a term used in the hydrologic and petroleum sciences to describe different processes and phenomena. The use of this term may be confusing unless clear definitions are provided and supporting information is given, such as the wettability of the porous medium. When applied to water in strongly water-wet porous media, the residual water saturation is the water saturation at which the water relative permeability decreases to zero. Theoretically, it refers to the water contained in the smallest pore spaces and as films on solid surfaces. The residual water saturation is the minimum water saturation that can be obtained by water drainage when there are no restrictions to drainage. In the petroleum industry, the residual water saturation is often referred to as the "irreducible water saturation."

When applied to NAPL in strongly water-wet porous media, the term residual NAPL has been used to describe the NAPL saturation that becomes occluded by water as water displaces NAPL into larger pore spaces during water imbibition (e.g., Steffy et al., 1997). Additionally, residual NAPL has been applied to describe the NAPL saturation in the vadose zone that does not drain from the pore spaces when no drainage restrictions exist (e.g., Oostrom and Lenhard, 2003). In strongly water-wet porous media, the NAPL is a nonwetting fluid. Typically, the vadose zone is assumed to be strongly water wet. The term residual NAPL has also been applied to mixed- and oil-wet porous media in the petroleum industry (e.g, Salathiel, 1973; Willhite, 1986). In those systems, the residual NAPL saturation is similar in concept to the residual water saturation in strongly water-wet porous media because it is the wetting fluid in the oil-wet pore spaces, and NAPL is contained in the smallest pore spaces available.

Given the descriptions in the previous paragraph, it is obvious that the term residual NAPL is used by various disciplines (e.g., petroleum vs. vadose zone hydrologists) to describe different processes. It has been used to describe either a nonwetting or a wetting fluid and is conceptualized as being continuous, discontinuous, or both throughout the pore spaces. In all cases, however, the residual NAPL is considered to be immobile. To avoid confusion with the terminology, we denote NAPL that is occluded by water as entrapped NAPL and NAPL that does not drain from the pore spaces as residual NAPL.

Laboratory experiments have shown that formation of residual NAPL might be an important process under certain conditions and needs to be considered in modeling NAPL migration, redistribution, dissolution, and volatilization. Jarsjö et al. (1994) and Hoag and Marley (1986) demonstrated that residual saturation was inversely correlated to aqueous saturation. Eckberg and Sunada (1984) reported that fine sands have larger residual NAPL saturations than coarser sands. Hofstee et al. (1997) observed relatively large amounts of the nonspreading NAPL perchloroethylene (PCE) retained in the vadose zone. Two recent intermediate-scale experimental studies (Oostrom and Lenhard, 2003; Oostrom et al., 2003) clearly demonstrate the formation of residual carbon tetrachloride saturation in heterogeneous porous media. Just like the PCE in Hofstee et al. (1997), the carbon tetrachloride in these experiments showed behavior consistent with a negative spreading coefficient, which can be defined for a NAPL spreading over the aqueous surface as (Adamson, 1982)

[1]
where {sigma}gl is the gas–aqueous, {sigma}nl the NAPL–aqueous, and {sigma}gn the gas–NAPL surface tension. If the spreading coefficient is zero or positive, then the NAPL spreads across the aqueous surface and the contact angles tend toward zero. If the spreading coefficient is negative, then the NAPL forms a lens on the aqueous–gas surface.

The NAPLs used in the experimental studies by Hofstee et al. (1997), Oostrom and Lenhard (2003), and Oostrom et al. (2003) are laboratory-grade organic liquids with relatively high oil–water interfacial tensions. Dwarakanath et al. (2002) found that interfacial tensions of several field NAPLs, including weathered gasoline and degreasing solvents, had values between 8 and 18 mN m–1, resulting in positive spreading coefficients. The low interfacial tensions are likely caused by surface-active agents. Dwarakanath et al. (2002) indicated that under vadose zone conditions, these NAPLs are not likely to exist in the form of discrete ganglia but rather will spread across the water–oil interface, resulting in low residual saturations.

Neglecting residual NAPL in multifluid subsurface flow and transport simulations has been noted as the primary reason for discrepancies between predicted and observed NAPL saturations. In modeling NAPL infiltration experiments, neglecting residual NAPL yields under estimations of the NAPL retained in the vadose zone and overestimation of the NAPL volumes reaching the capillary fringe. Van Geel and Sykes (1994a)(1994b) conducted two-dimensional laboratory experiments and numerical simulations, considering nonwetting fluid entrapment and saturation hysteresis of NAPL movement and redistribution for a controlled spill of n-heptane. Comparison of the numerical simulations against the laboratory experiment at later times showed that the modeled NAPL front had advanced deeper and NAPL saturation declined more rapidly in the vadose zone, indicating the need to include residual NAPL formation. Simulations of the experiments described by Hofstee et al. (1997), using a multifluid constitutive theory that did not include residual NAPL saturation formation, failed to predict NAPL retained in the vadose zone. Van Geel and Sykes (1997) modeled a two-dimensional heptane spill under fluctuating water table conditions and demonstrated the importance of accounting for residual NAPL formation. In review of work by Oostrom and Lenhard (1998), who compared predictions using different constitutive theories against one-dimensional NAPL infiltration experiments, Van Geel and Roy (2002) suggested that inclusion of a NAPL residual in the unsaturated zone may have improved the results for all of the tested retention relations. As shown in Oostrom and Lenhard (2003) and Oostrom et al. (2003), their experimental results could not be accurately predicted with a constitutive theory that did not account for residual saturation formation.

Several models allow the inclusion of a constant residual NAPL saturation value for a certain porous medium. For instance, the screening model by Weaver et al. (1994) for NAPL transport in the vadose zone uses a modified relative permeability term, allowing NAPL to flow only when its saturation is larger than an a priori chosen residual saturation value. Although such an approach might be acceptable for screening purposes, advanced models that allow the residual NAPL formation to be dependent on fluid saturation history are needed to predict NAPL movement more accurately. Recently, three theoretical models—Wipfler and van der Zee (2001), Van Geel and Roy (2002), and Lenhard et al. (2004)—have been published describing modifications to current constitutive relations to include the formation of residual NAPL. The three models consider the effects of fluid saturation history on the formation processes. Our main objective is to compare these advanced models against measured fluid saturations in both static pressure cell and transient column experimental results. Using the initials of the last names of the authors, we refer to these models as the WVDZ, VGR, and LOD models.


   THE THEORETICAL MODELS
TOP
 ABSTRACT
 INTRODUCTION
 THE THEORETICAL MODELS
COMPARISON WITH STATIC PRESSURE...
 COMPARISON WITH TRANSIENT COLUMN...
 SUMMARY
 REFERENCES
 
The three models all apply to water-wet porous media only. For more detailed descriptions, the reader is referred to the original manuscripts: Wipfler and van der Zee (2001), Van Geel and Roy (2002), and Lenhard et al. (2004). The symbol notations for fluid saturations vary considerably among the three models. To allow comparison between the models, the notations are converted to follow the definitions listed in this section. The apparent water, l, and total liquid saturation, t, are defined as

[2]

[3]
where Sl is the water saturation; Sn the NAPL saturation; Slr the residual aqueous phase saturation, which is assumed to be immobile; Sgel the entrapped air saturation in water; Sne the entrapped NAPL saturation in water; and Sge the total entrapped air saturation. Effective fluid saturations are defined as follows:

[4a]

[4b]

The total-entrapped air saturation is the sum of entrapped air occluded by water (Sgel) and by the NAPL (Sgen):

[5]
The relations between the various forms of NAPL discussed in this paper (free, entrapped, residual, mobile, and immobile) are given in Eq. [6], where the subscripts m, i, r, e, and f refer to mobile, immobile, residual, entrapped and free, respectively:

[6a]

[6b]

[6c]
The NAPL saturation is the sum of mobile and immobile NAPL (Eq. [6a]). Immobile NAPL is the sum of residual and entrapped NAPL (Eq. [6b]). Free NAPL, as opposed to entrapped NAPL, comprises mobile and residual NAPL and occurs within pore-space radii larger than those occupied by the apparent aqueous saturation, which is defined as the sum of the effective aqueous and entrapped-NAPL saturations (Eq. [6c]). Residual NAPL is the immobile component of free NAPL. It can exist as continuous or discontinuous liquid in the pore space and appears distinct from mobile NAPL when free NAPL is replaced by a fluid of lower wettability (e.g., gas phase). Mobile NAPL is the mobile component of free NAPL that only occurs as continuous liquid.

WVDZ Model
The saturation–pressure relationships in the WVDZ model are based on pore-scale fluid behavior of a nonspreading NAPL. The model does not include relative permeability–saturation relations. The model assumes a porous medium consisting of uniform spheres, with the centers situated at a 60° angle to each other in the horizontal and vertical directions. Each sphere has eight contact points. A drawing of the configuration is shown in Fig. 1 . The formation of residual NAPL in this model is illustrated by Wipfler and van der Zee (2001) by using the following scenario: Water is initially at residual saturation, in the form of thin layers of water spreading over the surfaces and capillary rings around the sphere contact points. An excess amount of a nonspreading NAPL is added and subsequently allowed to drain. During drainage, NAPL rings are formed around the water rings. The NAPL continues to drain until the rings become disconnected. The NAPL that remains in place is considered to be nondrainable residual NAPL. Figure 1 shows the fluid distribution just before the NAPL rings disconnect. The combined saturations of water and NAPL retained in the pendular rings when the NAPL becomes residual are defined as the critical total liquid saturation Stc. According to the WVDZ model, each fluid–porous medium combination has a unique value determined by fluid properties and porous medium geometry. When Sl = Slr, the residual NAPL saturation at Stc is equal to a maximum value, which can be expressed as Snrmax = Stc – Slr, or in terms of effective saturations, nrmax = tc. Adding water results in growth of the rings around the contact points, and some of the NAPL will become continuous and subsequently drain. For the situations where Slr < l ≤ tc, the residual NAPL saturation nr = tc – l. If, after drainage of the mobile NAPL until t = tc, the water saturation is reduced, the residual NAPL saturation will remain the same and t < tc. Wipfler and van der Zee showed through a calculation of fluid volumes that the value of the spreading coefficient has a large impact on the value of tc: The more negative the spreading coefficient, the larger the value of the critical liquid saturation.



View larger version (21K):
[in this window]
[in a new window]
  		Fig. 1. Pore-scale distribution of spheres and fluids according to the WVDZ model. After Fig. 2 in Wipfler and van der Zee (2001).

 


View larger version (15K):
[in this window]
[in a new window]
  		Fig. 2. The Pl–Pn field for unsaturated conditions with relations for Pnrc (dashed line) and Pnc (dotted line). Adapted from Fig. 3 in Wipfler and van der Zee (2001).

 
The pore-scale definition of tc is upscaled to the continuum scale by tying this saturation into the van Genuchten (1980) saturation–capillary pressure relations. The van Genuchten (1980) relations for water and total liquid saturations in three-phase systems are

[7a]

[7b]
where {alpha}, n, and m = 1 – 1/n are porous medium parameters, ßnl is the NAPL water scaling factor, ßgn is a air–NAPL scaling factor (Parker and Lenhard, 1987), Pnl is the NAPL–water capillary pressure, and Pgn is the gas–water capillary pressure. The capillary pressure Pij is defined as Pij = Pi – Pj. Inserting tc into Eq. [7b] and solving for the associated critical residual NAPL pressure, Pnrc, gives

[8]
For a certain porous medium and NAPL combination, Pnrc is a constant for Pl < 0. This relation is shown as a horizontal line in Fig. 2 , where a Pl–Pn plane for unsaturated conditions is shown. The figure, adapted from Wipfler and van der Zee (2001), also shows a sloped line for the critical NAPL pressure, Pnc. This relation is obtained by equating Eq. [6] and [7] and solving for the NAPL pressure. The critical NAPL pressure denotes the transition from a two-phase water–air to a three-phase water–NAPL–air system. The wedge between Pnrc and Pnc indicates the domain where the NAPL is residual.



View larger version (23K):
[in this window]
[in a new window]
  		Fig. 3. Predicted residual NAPL saturation as a function of tmax for a hypothetical case with (a) l = 0.0, (b) l = 0.3. For the VGR and LOD model predictions, nrmax = 0.2. For the WVDZ model, tc = 0.2.

 
VGR Models
The VGR models provide residual NAPL saturation formation terms that can be combined with existing saturation–capillary pressure models, such as the three-phase relations based on either the van Genuchten (1980) or the Brooks and Corey (1964) model. From the data obtained during a series of static three-phase experiments with residual heptane, Van Geel and Roy (2002) proposed a total of four conceptual models for describing residual NAPL: a linear model, a modified Land (1968) model, a constrained exponential model, and an unconstrained exponential model. The exponential models require the use of additional parameters that do not have a clear physical meaning and that have to be obtained through curve fitting of experimental data. For these reasons, the exponential models are not considered in this analysis. The linear model and the modified Land model will be referred to as the VGR Linear model and the VGR Land model, respectively.

For the VGR Linear model, the residual NAPL saturation, nr, is the product of four terms:

[9]
The A term is the absolute maximum residual NAPL saturation, nrmax, that can be obtained in a porous medium. This maximum value occurs physically after complete NAPL drainage from t = 1.0 to l = 0.0. Experimental techniques to determine this value are discussed in Van Geel and Roy (2002) and White et al. (2004). The B term is introduced to estimate a maximum residual NAPL saturation on the basis of the reversal point, or the historic maximum total liquid saturation, tmax, from the primary wetting to a drainage scanning curve. For the VGR Linear model, B = tmax. The C term is a scaling factor to account for the effect of water saturations larger than the residual water saturation on residual NAPL formation. In this model, C = /tmax. If the three-phase system is at residual water saturation, l = 0 and the C term is equal to 1. As the water saturation increases, this term decreases, thus reducing nr. The D term was applied to account for the position of the actual apparent total liquid saturation from the reversal point tmax. Van Geel and Roy proposed using D = /. The denominator of this term denotes the change in t that may occur in a three-phase system. Combining the four terms yields

[10]
In the VGR Land model, the residual NAPL saturation is a function of three terms:

[11]
The C and D terms are identical to the terms used in the VGR Linear model. The E term combines the information of the absolute maximum residual NAPL saturation and the total saturation reversal point according to a Land (1968)–type formulation: E = tmax/, where RL = – 1. Combining the E, C, and D terms yields:

[12]

In both VGR models, residual NAPL is formed during NAPL drainage from tmax. The total liquid saturation range for this formation to occur is between tmax and l. This concept is not compatible with the total liquid formation in Eq. [3]. Using that formulation, the minimum t on a NAPL drainage path cannot be less than the sum of the water saturation plus an amount of residual NAPL formed between tmax and t. The minimum t will never be able to reduce to l to form the maximum possible residual saturation between tmax and t. To allow residual NAPL formation in that range, the t definition for the VGR models is modified to

[13]
Van Geel and Roy (2002) did not specify explicit relative permeability–saturation relations but rather provided general recommendations on how to adjust these relations through the use of mobile NAPL saturation terms.

LOD Model
The LOD model also provides a modification to the existing three-phase saturation–capillary pressure models to account for residual NAPL formation. Besides saturation–capillary pressure relations, the LOD model also includes a comprehensive set of relative permeability–saturation equations. As opposed to the VGR models, where residual NAPL is formed during NAPL drainage, the LOD model assumes that residual NAPL is created in the form of immobile NAPL in small pores, pore wedges, films, or lenses during NAPL imbibition. Several figures illustrating the formation process are included in Lenhard et al. (2004). The LOD model for residual NAPL formation is a power function of the form

[14]
The F term is identical to the A term in the VGR models and is equal to nrmax. Lenhard et al. (2004) reasoned that the different behavior of spreading and nonspreading NAPL would be captured in this calibration term.

The term G is a measure of the volume fraction of pore space occupied by NAPL. Lenhard et al. (2004) proposed the relationship G = /. The term H is a measure of the size (radii) of the pores containing the NAPL; it is used to estimate the volume of NAPL that is immobile in thin films or lenses on water surfaces. According to Lenhard et al. (2004), H = 1 – l. Both terms G and H ensure that they are 1 when the residual NAPL saturation is at its maximum for a given water saturation path history and 0 when Snr = 0. Substituting the expressions for F, G, and H in Eq. [14] gives

[15]
To determine the power law exponent {kappa} in Eq. [15], Lenhard et al. (2004) hypothesized that immobile NAPL in pore wedges and in small pores will increase as the volume of pore space that contains NAPL increases. As the NAPL invades larger pores, small pores can become accessible to the NAPL that otherwise would not be accessible. However, the rate of increase in residual NAPL decreases as the volume of pore space that contains NAPL becomes larger. Lenhard et al. (2004) proposed using {kappa} = 0.5 as a first approximation. To determine the power law exponent {eta}, they hypothesized that the volume of immobile NAPL in thin films or lenses on water surfaces will be greater when NAPL is in smaller radii pore spaces relative to larger radii pore spaces. As smaller radii pores become accessible to the NAPL because of water drainage, the volume of immobile NAPL will increase at an increasing rate. Accordingly, Lenhard et al. (2004) proposed using {eta} = 2 on the basis of the assumption that the increase in the volume of NAPL present in thin films or as lenses on water surfaces as the water saturation decreases follows a pattern similar to the tortuosity variation as the water saturation decreases, as proposed by Burdine (1953). Substituting these exponents in Eq. [15] results in

[16]

Essentially, in the LOD model, the residual NAPL saturation is scaled to vary between the maximum residual NAPL saturation and 0, depending on where in the pore space NAPL was present. In contrast with the VGR models, the LOD model assumes that a residual NAPL is formed as soon as NAPL enters the porous medium.

The relative permeability relations in the LOD model were developed using the same approach used in Parker and Lenhard (1987) and Lenhard and Parker (1987). With respect to this model, the water and air permeabilities are identical and only the NAPL relative permeability relations are modified. The key assumption is that all of the residual NAPL is contained in pore spaces indexed between the apparent water saturation and the sum of the apparent water saturation and the effective residual NAPL saturation. This assumption was necessary to obtain an expedient method to predict the NAPL relative permeability since a well-defined theory describing the distribution of residual NAPL is not available. The reader is referred to Lenhard et al. (2004) for a detailed description of the relative permeability–saturation model. The fact that the LOD and VGR models assume that residual NAPL is formed when NAPL imbibes and when mobile NAPL is drained, respectively, has an effect on the NAPL relative permeability. For instance, assuming the same NAPL saturation, the NAPL relative permeability during NAPL imbibition is larger in the LOD model than in the VGR model.

Example Model Predictions
To demonstrate the residual NAPL saturation formation behavior of the theoretical models, two example calculations of static fluid distributions are presented. It is important to note that the results of both examples are independent of a saturation–capillary pressure model. The WVDZ model predicts a constant residual NAPL saturation equal to tc – l for conditions where tmax ≥ tc and a linearly increasing residual saturation from 0 to tc – l when 0 < tmax < tc. For the VGR Linear, VGR Land, and LOD models, Eq. [10], [12], and [16] are used to compute the residual NAPL saturations, respectively.

In the first example, residual NAPL saturations are computed as a function of the maximum total liquid saturation for the case where the water saturation is kept at residual. To evaluate comparable scenarios for all models, the maximum residual NAPL saturation, nrmax, for the VGR and LOD models, and critical total liquid saturation, tc, for the WVDZ model were assumed to be 0.2. No fluid entrapment is considered; hence, effective and apparent saturations are equal. The modeled scenario for the first example is as follows: Initially, the water saturation is at residual. Nonaqueous phase liquid is subsequently added until tmax = 1. Finally, mobile NAPL is drained from the system until only residual NAPL is left. The results of the first example are shown in Fig. 3a . For all models, the range in the residual NAPL saturation is 0 to 0.2. The plot shows that the WVDZ model predicts a constant residual NAPL saturation of 0.2 for cases where tmax ≥ tc and a linearly increasing residual saturation from 0 to 0.2 when 0 < tmax < tc. The VGR Linear model predicts a linear increase from 0 to 0.2 for the whole tmax range. The VGR Land and the LOD models predict decreasing rates of NAPL formation with increasing tmax. This result indicates that these models predict relatively larger amounts of residual NAPL in the smaller pores than in the larger pores.

For the second example, a similar scenario was followed except that the effective water saturation, l, was kept constant at 0.3. The results of the computations are shown in Fig. 3b. An important result is that the predicted residual NAPL saturation for the WVDZ model is 0 for all tmax values. For this scenario this means that the WVDZ model predicts that up to 70% of the pore space can be drained without leaving any residual NAPL behind. Since the water saturation is larger, the predicted residual NAPL saturations using the VGR and LOD models are smaller for the second example than the first example. The predicted maximum value at tmax = 1 for the LOD model is less than for the VGR models. The maximum values for the VGR Linear and Land model are the same. This can be seen by inserting tmax = 1 into Eq. [10] and [12]. For this case, both equations reduce to nr = nrmax.


   COMPARISON WITH STATIC PRESSURE CELL EXPERIMENTS
TOP
 ABSTRACT
 INTRODUCTION
 THE THEORETICAL MODELS
 COMPARISON WITH STATIC PRESSURE...
COMPARISON WITH TRANSIENT COLUMN...
 SUMMARY
 REFERENCES
 
As an initial test of the predictive equations for residual NAPL, predictions were compared with experimental results published by van Geel and Roy (2002). In this analysis, it is assumed that no air entrapment occurred. In the pressure cell experiments, a spreading NAPL (heptane mixed with Sudan III dye) was used. The residual NAPL saturation was measured at three initial different water saturations, (Sl = Slr, 0.2, and 0.4) on the main drainage saturation–capillary pressure branch with varying initial NAPL saturations. During NAPL drainage, some of the water also inadvertently moved out of the pressure cells, resulting in reduced final water saturations. A listing of the measured and computed values using the WVDZ, VGR, and LOD models is provided in Table 1 for the two sets of experiments (Parts I and II) where initial and final water saturations were reported. The residual saturations for the VGR Linear, VGR Land, and LOD models were computed using Eq. [10], [12], and [16], respectively. The nrmax, needed to compute the residual NAPL saturations according to these equations, was experimentally determined by Van Geel and Roy (2002) to be 0.20. The residual NAPL saturations for the WVDZ model were computed using nr = tc – l for tmax ≥ tc and nr = tmax – l for tmax < tc. Since no measured values of tc are available, residual saturations were computed using two arbitrary values: tc = 0.20 and 0.40. Because of the static nature of the experiments, the calculations are independent of the selected saturation–capillary pressure model.


View this table:
[in this window]
[in a new window]
  		Table 1. Measured and computed residual NAPL saturations for Van Geel and Roy's (2002) static pressure cell experiments.{dagger}

 
Measured and computed residual NAPL saturations for Parts I and II are shown in Table 1 and increase with increased tmax values. This trend is honored by the VGR and LOD models, but not by the WVDZ model. The WVDZ model predicts constant residual NAPL saturations for all measurements except for cases when tmax < tc. The computed sums of squared differences for the WVDZ model are therefore considerably larger than for the other methods (Table 2). Based on the information provided in Tables 1 and 2 it is obvious that the Van Geel and Roy (2002) experimental data do not support the concept of a critical total liquid saturation below which the NAPL becomes residual. Wipfler (2003) came to the same conclusion after completing a number of residual saturation formation experiments using the nonspreading NAPL dodecane. Wipfler stated that, "a critical total liquid saturation that indicates the transition to nondrainable NAPL saturation is not observed." The inability to find experimental evidence for the tc concept forwarded in the WVDZ model is casting doubt on attempts by Wipfler and van der Zee (2001) to model experimental three-phase fluid retention results by Hofstee et al. (1997). To predict observed trends in NAPL and water saturation behavior, Wipfler and van der Zee (2001) had to assume tc = 0.7, which is an unrealistically large value.


View this table:
[in this window]
[in a new window]
  		Table 2. Sum of squared differences between measured and model predictions.

 
The measured saturations indicate that most of the NAPL becomes residual in the smaller pores. This behavior is predicted more accurately by the VGR Land and the LOD models. The VGR Linear model underestimates the residual saturation in the smaller pores and consequently overestimates the nondrainable NAPL saturation in the larger pores. As a result, the sum of squared differences for the VGR Linear model is larger than this sum for the VGR Land and LOD models.


   COMPARISON WITH TRANSIENT COLUMN EXPERIMENTS
TOP
 ABSTRACT
 INTRODUCTION
 THE THEORETICAL MODELS
 COMPARISON WITH STATIC PRESSURE...
 COMPARISON WITH TRANSIENT COLUMN...
SUMMARY
 REFERENCES
 
To further test the theoretical models, model results were compared with measured fluid saturation from two transient column experiments. In the first experiment, Soltrol 220 (Chevron Phillips Chemical, The Woodlands, TX) was injected into a variably water saturated column. The second experiments involved infiltration and redistribution of dodecane into a draining column where the total liquid saturation was kept on a draining path. The fluid saturations were measured with a dual-energy {gamma} radiation system using procedures described by Oostrom and Dane (1990) and Oostrom et al. (1998). To enable this comparison, the VGR and LOD models were incorporated into the STOMP simulator (White and Oostrom, 2000, 2003). The models were incorporated into both the three-phase extension of the Brooks and Corey (1964)–Burdine (1953) relative permeability–saturation–capillary pressure model and the van Genuchten (1980)–Mualem (1976) model. The simulations presented in this paper use the VGR and LOD models in conjunction with the Brooks and Corey–Burdine relative permeability–saturation–capillary pressure relations. The WVDZ model was not implemented into the STOMP simulator because its major premise, the existence of a total liquid saturation below which NAPL becomes residual, appeared to be invalid for the static experiments reported by Van Geel and Roy (2002) and Wipfler (2003). In addition, the description of the WVDZ model by Wipfler and van der Zee (2001) is not sufficiently detailed to allow for all necessary phase condition changes to occur in the simulator.

The STOMP simulator solves systems of coupled nonlinear equations, which are solved using linearization techniques (e.g., Newton–Raphson schemes). Models for residual NAPL would preferably have smooth and continuous partial derivatives of the residual NAPL with respect to the primary system unknowns. A practical model for predicting the formation and behavior of residual NAPL saturation must show agreement with experimental observations and be amenable to incorporation into a coupled nonlinear equation solution scheme. Complicating the numerical implementation is that the VGR and LOD models do not provide for the disappearance or reduction of residual or entrapped NAPL via aqueous dissolution or evaporation. Moreover, the complexities that develop in the hysteretic constitutive relations for regions that experience multiple wetting and draining paths in the Lenhard et al. (2004) and the van Geel and Roy (2002) theories do not lend themselves to being easily resolved as constitutive equations within a set of nonlinear conservation equations. Therefore, the VGR and LOD models have been implemented with practicality in mind. The model applies to multifluid flow and transport in water-wet porous media and accounts for entrapped and residual NAPL saturation formation, allows for dissolution, but avoids the numerical complexities of pore-geometry hysteresis and wettability alterations. For details of the LOD model implementation into the STOMP simulator, the reader is referred to White et al. (2004). The VGR model incorporation was completed using the same numerical methods. To contrast the results obtained using the residual saturation formation models with a nonhysteretic model, both transient experiments were also simulated using the nonhysteretic Brooks and Corey–Burdine constitutive model, indicated in this paper as NH simulations.

Column Experiment 1
Experiment 1 was conducted in a 1-m-long column with a 5.0-cm inside diameter. The column was packed under water-saturated conditions with 90 cm of a medium-grained Hanford sand. Denoting the bottom as z = 0 cm, the column was calibrated for dual-energy {gamma} radiation measured at nine locations: z = 5, 15, 25, 35, 45, 55, 65, 75, and 85 cm. The water table was subsequently lowered at a rate of 10 cm h–1 from z = 90 to 20 cm, after which the porous medium system was allowed to reach quasi-static equilibrium. The column was scanned daily and after 7 d no temporal changes in aqueous saturation were observed. After this draining period, 120 mL of Soltrol 220 was injected at a rate of 1 mL min–1 for 2 h. After the injection, the NAPL was allowed to redistribute. The column was scanned daily for the duration of the experiment (50 d) to determine NAPL and water saturation. Measured hydraulic properties of the porous medium are listed in Table 3, while fluid properties are listed in Table 4. For details of the measurement techniques, the reader is referred to White et al. (2004).


View this table:
[in this window]
[in a new window]
  		Table 3. Hydraulic properties of porous media used in the two transient column experiments.

 

View this table:
[in this window]
[in a new window]
  		Table 4. Fluid properties of the NAPLs used in the transient column experiments.

 
The Soltrol 220 imbibition experiment was modeled starting with a drained column (i.e., hydrostatic aqueous phase) with the water table location at z = 20.0 cm. The domain was discretized into 90 1-cm grid cells. A total of 120 mL of NAPL was then injected during a 2-h period in the top grid cell. The infiltration and redistribution of the NAPL was simulated for 50 d using VGR models, the LOD model, and a standard Brooks Corey–Burdine nonhysteretic relative permeability–saturation–capillary pressure model (NH model).

Comparisons between the numerical simulations and laboratory experimental results are shown at z = 85 in Fig. 4 and z = 25 cm in Fig. 5 . Figure 6 shows the final measured and computed distribution of the NAPL in the column after 50 d. The sums of the squared differences between predicted and all measured NAPL and water saturations at all locations are shown in Table 5.



View larger version (30K):
[in this window]
[in a new window]
  		Fig. 4. Measured and simulated fluid saturations at z = 85 cm for Exp. 1.

 


View larger version (31K):
[in this window]
[in a new window]
  		Fig. 5. Measured and simulated fluid saturations at z = 25 cm for Exp. 1.

 


View larger version (17K):
[in this window]
[in a new window]
  		Fig. 6. Measured and simulated (a) NAPL and (b) water saturations at the conclusion of Exp. 1.

 

View this table:
[in this window]
[in a new window]
  		Table 5. Sum of squared differences between measured and simulated fluid saturations using the NH, LOD, VGR Linear and VGR Land methods for Exp. 1.{dagger}

 
At locations with unsaturated conditions, a rapid increase in NAPL saturation was observed, followed by a more gradual decrease. For instance, at z = 85 cm, the NAPL saturation increases rapidly to a value of approximately 0.7 before it slowly decreases to a value close to 0.15 at t = 50 d. For the locations with unsaturated conditions, the NH simulation predicts earlier NAPL arrival times than the simulations that account for residual saturation formation. Although the simulations accounting for residual saturation formation generally perform better than the NH simulation, the discrepancies with experimentally obtained NAPL saturations are the largest at intermediate levels. These differences can be explained in part by the observation that at these locations, the predicted maximum saturation is less than the experimental maximum. Another reason explaining the differences might be the inherent heterogeneities in the column. At the saturated locations in the lower part of the column, the NH simulation predicts a premature arrival of the NAPL and water drainage (Fig. 5). Since no residual saturation is formed in the NH simulation, more NAPL is available for downward movement, resulting in larger relative permeability values and faster movement. The formation of residual NAPL reduces the NAPL relative permeability, thus slowing downward migration and resulting in retention of NAPL in the upper region of the column as it is unable to migrate to the lower elevations.

At all measurement levels, the simulation results considering residual NAPL showed better agreement with the experimental results than the NH model (Fig. 6). This visual indication is supported by the values in Table 5. The sums of squared differences for both NAPL and water are consistently larger for the NH model than for the other methods. The differences between the VGR and LOD models are relatively small. The VGR Land models recorded the lowest sum of squared differences for both fluids, followed by the LOD and the VGR Linear models.

Column Experiment 2
Experiment 2 was conducted in a 95-cm-long glass column with a 7.62-cm inside diameter. The column was filled with 85 cm of a medium-grained sand. The sand was a mixture of three Unimin sands: 40% (by weight) F30, 30% F40, and 30% F50-60. Denoting the bottom of the sand as z = 0 cm and the top as z = 85 cm, the column was calibrated for dual-energy {gamma} radiation measured at eight locations: z = 5, 15, 25, 35, 45, 55, 65, and 75 cm. After the water table was set level with the top of the sand, 150 mL of dodecane was placed on top of the saturated sand. The initial dodecane thickness on top of the saturated sand was 3.3 cm. The water table was then lowered from z = 85 to 15 cm in seven 10-cm steps imposed on seven consecutive days. On each of the 7 d, the water level was lowered at a rate of 0.5 cm min–1. The column was scanned daily for the duration of the experiment (35 d) to determine fluid saturations. The experiment was designed to avoid gas entrapment by either the NAPL or the aqueous phase.

Porous medium and fluid properties are listed in Tables 3 and 4, respectively. The methods to determine these properties are described in White et al. (2004). The experiment was simulated for 35 d with VGR, LOD, and NH models. For this experiment, the VGR Linear and VGR Land models produce identical results because at each locations that encountered NAPL, tmax = 1. Equations [10] and [12] become identical when this scenario occurs.

In contrast to Soltrol 220, dodecane is a nonspreading NAPL with a spreading coefficient of –0.006 N m–1. Contacts angles, computed with an algorithm discussed by White et al. (2004), yielded a NAPL–aqueous phase contact angle of 12.5° and a gas–NAPL contact angle of 32.3° (Table 4). The nonspreading nature of this NAPL is also evidenced by the experimentally obtained nrmax value of 0.232, which is more than twice as large as for Soltrol 220 in Exp. 1.

Comparisons between the numerical simulations and laboratory experimental results are shown at z = 75 cm in Fig. 7 and z = 25 cm in Fig. 8 . Figure 9 shows the final measured and computed distribution of the NAPL in the column after 35 d. The sums of the squared differences between predicted and measured fluid saturations are shown in Table 6.



View larger version (25K):
[in this window]
[in a new window]
  		Fig. 7. Measured and simulated fluid saturations at z = 75 cm for Exp. 2.

 


View larger version (24K):
[in this window]
[in a new window]
  		Fig. 8. Measured and simulated fluid saturations at z = 25 cm for Exp. 2.

 


View larger version (16K):
[in this window]
[in a new window]
  		Fig. 9. Measured and simulated (a) NAPL and (b) water saturations at the conclusion of Exp. 2.

 

View this table:
[in this window]
[in a new window]
  		Table 6. Sum of squared differences between measured and simulated fluid saturations using the NH, LOD, and the VGR methods for Exp. 2.{dagger}

 
At all observation levels, the simulation results considering residual NAPL show better agreement with the experimental results. The figures clearly show the stepwise lowering of the lower aqueous-phase boundary condition. In agreement with the experimental observations, no NAPL was predicted to enter the system after the first water pressure head lowering from z = 85 to 75 cm. Nonaqueous phase liquid started to infiltrate the domain only after the second 10-cm reduction in pressure head from z = 75 to 65 cm. At z = 75 cm (Fig. 7), the NAPL saturation increased until the end of Day 2 and dropped after the pressure head was lowered from z = 65 to 55 cm. The predicted saturations using the VGR and LOD models and measured NAPL saturations at this location are virtually time-invariant after about 5 d until the end of the experiment. This behavior is consistent with the nonspreading characteristics of the dodecane. The NH simulation consistently underpredicts the NAPL saturation when the NAPL is draining from the column. In this simulation, NAPL is always mobile and more NAPL is allowed to move downwards in the column. Figure 8 (z = 25 cm) represents a location that remains fully liquid saturated during the experiment. The NAPL arrival time and saturation values are predicted well with the VGR and LOD simulations, but the NH simulation underestimates the arrival time and allows more NAPL to move downwards into the liquid-saturated zones. These results were also observed in the other column experiments and indicate that an erroneous assumption of no residual saturation formation in the vadose might have implications not only for the vadose zone but also for the underlying saturated regions. The saturations at the end of the experiment (Fig. 9) illustrate this clearly. The NH model predicts almost complete drainage of NAPL from the vadose zone and collection of the NAPL at the top of the saturated region. The experiment results and the VGR and LOD model simulations predict considerable residual NAPL saturation in the vadose zone and less mobile NAPL in the saturated region of the column. The sums of squared differences for NAPL saturations (Table 6) clearly show that the NH model produces much larger sums than the VGR and LOD simulations. As was observed for Exp. 1, the sums of squared differences for the LOD and VGR were similar. In contrast with the results for Exp. 1, the LOD model performed slightly better than the VGR model at the majority of the locations. However, the results of the two transient experiments do not clearly favor one of the residual NAPL saturation models.


   SUMMARY
TOP
 ABSTRACT
 INTRODUCTION
 THE THEORETICAL MODELS
 COMPARISON WITH STATIC PRESSURE...
 COMPARISON WITH TRANSIENT COLUMN...
 SUMMARY
REFERENCES
 
Recently, the WVDZ, VGR Linear, VGR Land, and LOD models have been proposed, providing constitutive relations for residual NAPL formation in the vadose zone. The WVDZ model assumes a critical total liquid saturation below which all NAPL becomes residual. The VGR and LOD models are extensions of an existing hysteretic relative permeability–saturation–capillary pressure model and assume formation of residual NAPL during NAPL drainage and imbibition, respectively.

Following a description of each model, model predictions were compared against results of a series of static pressure cell experiments. No experimental evidence was found supporting the WVDZ concept of a critical total liquid saturation. The other two models produced acceptable predictions.

The VGR and LOD models were incorporated into a multifluid flow simulator and model predictions were compared against results from two transient column experiments. The WVDZ model was not implemented because of its lack of agreement with the static experiments and incomplete description of phase conditions changes required for numerical simulation. Both the VGR and LOD models performed considerably better than the nonhysteretic simulations. An additional important observation was that neglecting a residual NAPL term in the simulator tends to overestimate NAPL volumes reaching the water table. The improvements obtained with the VGR and LOD models underwrite the importance of incorporating the residual NAPL saturation formation process in multifluid flow simulators. Although the VGR and LOD models are conceptually different, no clear performance differences could be observed.


   ACKNOWLEDGMENTS
 
Pacific Northwest National Laboratory (PNNL) is operated by the Battelle Memorial Institute for the Department of Energy (DOE) under Contract DE-AC06-76RLO 1830. This research is part of the Groundwater/Vadose Zone Integration Project funded through the DOE's Richland Operations Office. A portion of the research described in this paper was performed in the Environmental Molecular Sciences Laboratory, a national scientific user facility sponsored by the DOE's Office of Biological and Environmental Research and located at PNNL. R.J. Lenhard's contribution was supported by the Laboratory Directed Research & Development Program of the Idaho National Engineering & Environmental Laboratory (INEEL) under the Subsurface Science Initiative.


   REFERENCES
TOP
 ABSTRACT
 INTRODUCTION
 THE THEORETICAL MODELS
 COMPARISON WITH STATIC PRESSURE...
 COMPARISON WITH TRANSIENT COLUMN...
 SUMMARY
 REFERENCES