Effect of soil moisture dynamics on dense nonaqueous phase liquid (DNAPL) spill zone architecture in heterogeneous porous media

Journal of Contaminant Hydrology 90 (2007) 159 – 183 www.elsevier.com/locate/jconhyd

Effect of soil moisture dynamics on dense nonaqueous phase liquid (DNAPL) spill zone architecture in heterogeneous porous media Hongkyu Yoon ⁎, Albert J. Valocchi, Charles J. Werth Department of Civil and Environmental Engineering, Environmental Engineering and Science, University of Illinois at Urbana-Champaign, 205 N. Mathews Avenue, Urbana IL 61801, USA Received 11 June 2005; received in revised form 24 May 2006; accepted 4 September 2006 Available online 20 December 2006

Abstract The amount, location, and form of NAPL in contaminated vadose zones are controlled by the spatial distribution of water saturation and soil permeability, the NAPL spill scenario, water infiltration events, and vapor transport. To evaluate the effects of these processes, we used the three-phase flow simulator STOMP, which includes a new permeability–liquid saturation–capillary pressure (k–S–P) constitutive model. This new constitutive model considers three NAPL forms: free, residual, and trapped. A 2-D vertical crosssection with five stratigraphic layers was assumed, and simulations were performed for seven cases. The conceptual model of the soil heterogeneity was based upon the stratigraphy at the Hanford carbon tetrachloride (CT) spill site. Some cases considered co-disposal of NAPL with large volumes of wastewater, as also occurred at the Hanford CT site. In these cases, the form and location of NAPL were most strongly influenced by high water discharge rates and NAPL evaporation to the atmosphere. In order to investigate the impact of heterogeneity, the hydraulic conductivity within the lower permeability layer was modeled as a realization of a random field having three different classes. For six extreme cases of 100 realizations, the CT mass that reached the water table varied by a factor of two, and was primarily controlled by the degree of lateral connectivity of the low conductivity class within the lowest permeability layer. The grid size at the top boundary had a dramatic impact on NAPL diffusive flux just after the spill event when the NAPL was present near the ground surface. NAPL evaporation with a fine grid spacing at the top boundary decreased CT mass that reached the water table by 74%, compared to the case with a coarse grid spacing, while barometric pumping had a marginal effect for the case of a continuous NAPL spill scenario considered in this work. For low water infiltration rate scenarios, the distribution of water content prior to a NAPL spill event decreased CT mass that reached the water table by 98% and had a significant impact on the formation of trapped NAPL. For all cases simulated, use of the new constitutive model that allows the formation of ⁎ Corresponding author. Fax: +1 217 333 0687. E-mail address: [email protected] (H. Yoon). 0169-7722/$ - see front matter © 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.jconhyd.2006.09.007

160

H. Yoon et al. / Journal of Contaminant Hydrology 90 (2007) 159–183

residual NAPL increased the amount of NAPL retained in the vadose zone. Density-driven advective gas flow from the ground surface controlled vapor migration in strongly anisotropic layers, causing NAPL mass flux to the lower layer to be reduced. These simulations indicate that consideration of the formation of residual and trapped NAPLs and dynamic boundary conditions (e.g., areas, rates, and periods of different NAPL and water discharge and fluctuations of atmospheric pressure) in the context of full three-phase flow are needed, especially for NAPL spill events at the ground surface. In addition, NAPL evaporation, densitydriven gas advection, and NAPL vertical movement enhanced by water flow must be considered in order to predict NAPL distribution and migration in the vadose zone. © 2006 Elsevier B.V. All rights reserved. Keywords: Nonaqueous phase liquid (NAPL); Vapor transport; NAPL evaporation; STOMP; Hanford site

1. Introduction Dense nonaqueous phase liquids (DNAPLs), like chlorinated solvents, are ubiquitous soil and groundwater contaminants at many hazardous waste sites. DNAPLs are difficult to remediate because they are hard to locate, persistent over long time periods, and sorb strongly to soils and sediments (Pankow and Cherry, 1996). DNAPLs are present in both the vadose and saturated zones and act as a continuing source of groundwater contamination for very long periods. Our focus in this study is on the vadose zone, where uncertainties in predicting DNAPL distribution and migration create challenges for the design of proper remediation strategies (DOE, 2001). These uncertainties have been attributed to the heterogeneity and spatial variability of hydraulic properties (Poulsen and Kueper, 1992; Essaid et al., 1993), hysteresis with NAPL entrapment (Essaid et al., 1993; Van Geel and Sykes, 1997; Oostrom and Lenhard, 1998), and deficiencies in understanding of NAPL migration processes (DOE, 2001). Controlled laboratory experiments of NAPL infiltration and redistribution in layered porous media conducted by several investigators (Hofstee et al., 1997, 1998; Oostrom et al., 2003a) showed that existing multiphase flow models under-predicted NAPL saturation, primarily because residual NAPL was not considered. Several recent works have improved understanding of residual NAPL formation (Wipfler and van der Zee, 2001; Van Geel and Roy, 2002; Keller and Chen, 2003; Lenhard et al., 2004). Lenhard et al. (2004) proposed a new relative permeability–saturation–capillary pressure (k–S–P) relation that accounts for the formation of residual and trapped NAPLs in a threephase system. This new constitutive relation considers three NAPL forms: free, residual, and trapped. Both free and residual NAPL have direct access to the pore gas; free NAPL is mobile, but residual NAPL is immobile because it is held by capillary forces. Trapped NAPL is also immobile, but it is surrounded by water. The new hysteretic k–S–P model has been incorporated into a continuum-based multiphase flow simulator (STOMP) and verified against column experimental results (White et al., 2004; Oostrom et al., 2005). Additional difficulties in understanding physical processes that affect NAPL distribution in the vadose zone arise when sites have been contaminated with complex mixtures along with high volumes of discharged wastewater. A notable example is the Hanford Site in Washington where carbon tetrachloride (CT) was discharged to the subsurface with large volumes of wastewater. Based on a rough estimation before remediation began in 1992, approximately 64% of the original inventory was assumed to remain as NAPL phase in the vadose zone (Last and Rohay, 1993). Nevertheless, only 10% of the original CT discharged was removed from the vadose zone using soil vapor extraction between 1992 and 1999 (Rohay, 2000) and only 1% of the CT discharged was removed from groundwater using a pump and treat system between 1996 and

H. Yoon et al. / Journal of Contaminant Hydrology 90 (2007) 159–183

161

2003 (Poston et al., 2004). Since NAPL phase CT has not been directly observed in the vadose and saturated zones, a new conceptual model is needed to account for the physical processes that hinder the recovery of the remaining CT and to better describe mass balance of the total CT discharged (DOE, 2003). An experimental study in a two-dimensional intermediate scale flow cell by Oostrom and Lenhard (2003) showed that a large portion of the CT discharged can be trapped if high water discharge follows the initial NAPL discharge. Since the trapped NAPL can persist due to diffusion through pore water during soil vapor extraction, it is important to predict the formation of trapped NAPL for scenarios with large volumes of discharged wastewater. Recently, three-dimensional modeling using a three-phase simulator, STOMP, was conducted to simulate the large-scale features of CT distribution at the Hanford Site (Oostrom et al., 2003b). Sensitivity simulations suggest that several parameters including permeability, hydraulic properties, and the value of maximum residual NAPL saturation affect the NAPL distribution. Also, a large portion of the initial CT mass was expected to exist as NAPL phase in the vadose zone, but CT mass predicted by simulations has not been verified by the current inventory model (DOE, 2003) and the formation of the trapped NAPL has not been evaluated. Further, the mechanisms that affect the migration and distribution of NAPL under large volumes of discharged wastewater were not thoroughly evaluated. Although NAPL flow and transport in the vadose zone occur in the presence of three mobile phases (water, gas, and NAPL), most multiphase flow simulations assumed that gas was at atmospheric pressure. Several studies have shown that transport in the gas phase can significantly impact the contaminant mass distribution in the presence of NAPL. Transport processes considered were density-driven advective gas flow (Sleep and Sykes, 1989; Mendoza and Frind, 1990a; Lenhard et al., 1995), gas advection enhanced by water infiltration (Sleep and Sykes, 1993; Zhu and Sykes, 2000), and gas diffusion to the atmosphere (Thomson et al., 1997; Zhu and Sykes, 2000, Jellali et al., 2003, Broholm et al., 2005). A recent field experiment in the vadose zone investigated the fate of a fuel-NAPL mixture emplaced 0.8–1.3 m below ground surface (Broholm et al., 2005). Analysis within the NAPL source zone revealed that mole fractions of several volatile compounds were depleted within about 3 months due to gas diffusion to the atmosphere. Experimental and modeling studies also indicate that atmospheric pressure fluctuations affect the advective vapor flux through the soil surface (Choi et al., 2002; Tillman et al., 2003) and gas flow to a well in the vadose zone (Ellerd et al., 1999; Neeper, 2002). However, the effects of these fluctuations on DNAPL migration and distribution in the vadose zone have not been investigated. Although volatile organic chemicals (e.g., chlorinated methanes, ethanes, and ethenes) can evaporate into the atmosphere at the ground surface during DNAPL and water spill events, the importance of DNAPL and water spills at the ground surface on vapor movement and mass balance calculations have been neither particularly emphasized nor thoroughly evaluated in previous investigations. The goal of this work is to identify what mechanisms control the migration and distribution of NAPL in a strongly layered heterogeneous vadose zone under dynamic soil–moisture conditions. Our work was motivated by an assessment of CT in the vadose zone at the Hanford site, so many of our simulation parameters are representative of this site. However, we do not intend to simulate actual site conditions at Hanford. Seven different scenarios were studied by numerical simulation to illustrate the effect of heterogeneity, gas advection, NAPL evaporation, barometric pumping, water discharge rates, and water saturation on the distribution and migration of NAPL in heterogeneous porous media for different NAPL and water loading histories (i.e., dynamic soil– moisture conditions). To illustrate important mechanisms related to the Hanford CT site, we emphasize (1) the importance of vapor movement in the upper layers on NAPL migration, (2) the

162

H. Yoon et al. / Journal of Contaminant Hydrology 90 (2007) 159–183

effect of heterogeneity on NAPL migration in a deep vadose zone, and (3) the formation of trapped and residual NAPL. We used the three-phase flow simulator STOMP, which includes the new k–S–P constitutive model developed by Lenhard et al. (2004). An enhanced understanding of the mechanisms that control NAPL behavior under dynamic soil–moisture conditions, like those at the Hanford site, will give insights into developing conceptual models for contaminant distributions at strongly layered heterogeneous sites. 2. Mathematical model Numerical simulations were conducted using the STOMP (Subsurface Transport Over Multiple Phases) simulator (White and Oostrom, 2000, 2004). Since a volatile organic compound, CT is considered in this study, the Water–Air–Oil mode of the STOMP simulator was used to simulate NAPL transport in the gas phase and changes in gas pressure. For the Water–Air–Oil mode, the STOMP simulator solves mass conservation equations for components simultaneously (water, air, and oil) that describe multiphase flow and transport through variably saturated porous media. Phase partitioning of each component is computed assuming equilibrium conditions. The general conservation of mass equations for water, air and oil components in a three-phase system can be written as: " #  X X A X Mw w w w w ð/Sa xa qa Þ þ ½jd xa qa qa − jd /Sa qa ðsa Da þ Dha Þjxa At a¼1;g Ma a¼1;g a¼g;1 X ¼ xw a qa Qa

ð1Þ

" #  X X A X Ma a a a a ð/Sa xa qa Þ þ ½jd xa qa qa − jd /Sa qa ðsa Da þ Dha Þjxa At a¼1;g Ma a¼1;g a¼1;g X ¼ xaa qa Qa

ð2Þ

a¼1;g

a¼1;g

A At

"

X

a¼1;g;n

−

# ð/Sa xoa qa Þ

X

a¼1;g

þ

ð1−/Þxos qs

þ

X

½jd xoa qa qa 

ð3Þ

a¼1;g;n

 X Mo o o jd /Sa qa ðsa Da þ Dha Þjxa ¼ xoa qa Qa Ma a¼1;g;n

where α denotes the phase (l = aqueous, g = gas, n = NAPL, s = solid), ϕ is porosity, Sα is the saturation of phase α, ρα is the density of phase α, ρs is porous media grain density, xαc is the mole fraction of component c (w = water, a = air, o = oil) in phase α, qα is the Darcy velocity of phase α, τα is the tortuosity of phase α (Millington and Quirk, 1961), Dαc is the molecular diffusion coefficient for component c in phase α, Dhα is the hydraulic mechanical dispersion coefficient of phase α, M is the molecular weight of component c or phase α, and Qα is the source or sink term

H. Yoon et al. / Journal of Contaminant Hydrology 90 (2007) 159–183

163

for phase α per unit of bulk volume. The Darcy velocity of phase α is computed from Darcy's law, kra k qa ¼ − ðjPa þ qa gzÞ ð4Þ la where k is the intrinsic permeability tensor, krα is the relative permeability of phase α, μα is the viscosity of phase α, Pα is the pressure of phase α, g is the gravitational constant, z is a unit vector in the direction of gravitation. Constitutive relations are needed to close the above set of equations. In particular, the constitutive k–S–P relation is a critical component of predictive multiphase flow models (White et al., 1995). A complete description of constitutive relations employed in the STOMP simulator can be found in White and Oostrom (2000, 2004) and White et al. (2004). In this section, the k–S–P relations including the formation of residual NAPL are described briefly. In the k–S–P model updated by Lenhard et al. (2004), the total-NAPL saturation (Sn) consists of three components: free (Snf), residual (Snr), and trapped NAPL (Snt): Sn ¼ Snf þ Snr þ Snt

ð5Þ

The free NAPL consists of mobile NAPL that can exist at water–gas interfaces. The trapped NAPL is discontinuous and immobile, and occluded by water, the most wetting phase. The residual NAPL is not trapped by the aqueous phase, but it is immobile. Trapped NAPL saturation is computed using the parametric model of Kaluarachchi and Parker (1992). In the residual NAPL formation theory developed by Lenhard et al. (2004), the residual NAPL saturation is not constant, but estimated as a function of saturation-path history. It is nonlinearly scaled between zero and the maximum residual NAPL saturation. The effective residual NAPL saturation (S¯nr) is computed as max ¼ max

Snr ¼ S nr ðS t

¼

¼

− S 1 Þ0:5 ð1− S 1 Þ1:5

ð6Þ ¼ max

max where S¯nr is the maximum effective residual NAPL saturation, S t is the historic maximum ¼ apparent total-liquid saturation, and S 1 is the apparent aqueous saturation. Effective saturations ¼ (S¯i) and apparent saturations (S i ) ignoring air entrapment are expressed as

S¯1 ¼ S1 −S1r ; 1−S1r ¼

S1 ¼

S¯n ¼ Sn ; 1−S1r

S1 þ Snt −S1r ; 1−S1r

¼

St ¼

S¯g ¼ Sg 1−S1r

ð7Þ

S1 þ Sn −S1r 1−S1r

ð8Þ ¼

where S1r is the residual aqueous saturation and S t is the apparent total-liquid saturation. Since preliminary tests showed that air entrapment had no significant impact on water and NAPL migration under the conditions considered in this study, air entrapment was not considered. The maximum residual NAPL saturation can be observed after complete drainage of NAPL from a column initially saturated with NAPL at the residual aqueous saturation (i.e. water wetting system). It was observed by Van Geel and Roy (2002) that the residual NAPL saturation increases as the maximum NAPL saturation prior to NAPL drainage increases. The Brooks–Corey S–P model was used in this study because this model considers a distinct nonwetting fluid entry pressure (Brooks and Corey, 1964). The Brooks–Corey S–P model with the Burdine relative permeability model including fluid entrapment predicted

164

H. Yoon et al. / Journal of Contaminant Hydrology 90 (2007) 159–183

experimental data in the unsaturated and saturated regions better than both the nonhysteretic and hysteretic van-Genuchten S–P models with the Mualem relative permeability model (Oostrom and Lenhard, 1998). Since the hysteretic Brooks–Corey model is not available in the current STOMP, the Brooks and Corey S–P model with the Burdine relative permeability model including fluid entrapment was used. In the new k–S–P relation updated by Lenhard et al. (2004), the relative permeability for NAPL becomes zero as the free NAPL saturation becomes zero. Since residual NAPL is immobile, residual NAPL saturation does not contribute to NAPL relative permeability. The aqueous, NAPL, and gas relative permeabilities are expressed as 2þ3k kr1 ¼ ðS¯1 Þð k Þ

ð9Þ

  2  ¼  2þk  ¼ ð2þk k Þ ðkÞ ¯ ¯ S S krn ¼ nf − S 1 þ nr St

ð10Þ

 2   ¼  2þk  ð Þ 1− S t k krg ¼ S¯g

ð11Þ

where krl, krn, and krg are the aqueous, NAPL, and gas relative permeabilities and λ is the pore size distribution index. Table 1 Porous media properties used in simulations(1) Stratigraphic k (3) (m2) units(2)

Depth Variance Range (m) (m)(4)

HUC Mean Low HFS Mean High Low PPlz Mean High Low PPlc Mean High RG Mean

12

–

–

15.5

0.12

45, 1.5

3

0.12

15, 1.2

4

0.25

9, 0.9

–

–

(1)

6.8E-11 6.6E-13 6.6E-12 6.6E-11 1.7E-14 1.7E-13 1.7E-12 7.8E-13 7.8E-12 7.8E-11 1.3E-11

22

Composition (%)(5)

ϕ(6)

λ(7)

ψ d (8) (cm)

100 20 60 20 15 70 15 20 60 20 100

0.166 0.382 0.382 0.382 0.424 0.52 0.424 0.337 0.337 0.337 0.226

1.0 8.8 0.503 25.4 0.503 15.4 0.503 10.4 0.855 143.7 0.855 123.7 0.855 103.7 0.720 60.4 0.720 40.4 0.720 20.4 0.538 50.4

Slr

(9)

0.12 0.20 0.15 0.10 0.40 0.25 0.20 0.25 0.25 0.25 0.10

Snt,max (10) K d (11) (L/kg) 0.10 0.25 0.20 0.15 0.25 0.25 0.18 0.25 0.20 0.15 0.10

0.1 0.1

0.5

0.1 0.1

Maximum residual NAPL saturation was set at 0.05. (2) Indicator classes are based on the three soil classes and no nugget effect. Low and high values represent permeability an order of magnitude lower and higher than a mean value in heterogeneous layer, respectively. (3) Permeability, based on Khaleel and Freeman (1995) and Oostrom et al. (2003b). (4) Horizontal and vertical ranges. (5) Composition of the soil classes for each layer. (6) Porosity, based on Khaleel and Freeman (1995). (7) Brooks and Corey pore size distribution index, based on Khaleel and Freeman (1995) and Yoon et al. (2003). (8) Air entry pressure, based on Khaleel and Freeman (1995) and Yoon et al. (2003). (9) Residual water sat., based on Khaleel and Freeman (1995). (10) Maximum trapped NAPL sat., based on Essaid et al. (1993), Oostrom and Lenhard (2003), and White et al. (2004). (11) Linear sorption coefficient, based on Yonge et al. (1996) and Truex et al. (2001). (7) and (8) The k–S–P parameters were converted to equivalent to Brooks–Corey–Burdine parameters using Lenhard et al. (1989).

H. Yoon et al. / Journal of Contaminant Hydrology 90 (2007) 159–183

165

3. Numerical simulations 3.1. Hanford carbon tetrachloride site Modeling efforts focus on capturing the effects of important stratigraphic features and operational parameters at the Hanford Site, but we do not attempt to mimic precise site-specific conditions. The Hanford Site near the carbon tetrachloride (CT) discharge areas (e.g., the 216-Z-9 Trench) consists of five major stratified layers in the vadose zone that span a depth from 60 to 70 m. These include, from top to bottom, the Hanford upper coarse (HUC), the Hanford fine sand (HFS), the Plio-Pleistocene Z (PPlz), the Plio-Pleistocene C (PPlc) (often both layers are referred to as the caliche layer), and the Ringold Unit E (RG). Hydraulic properties of the porous media used in this study are listed in Table 1. For the k–S–P parameters, the equivalent Brooks–Corey parameters were obtained using data reported in Khaleel and Freeman (1995) by adopting the conversion equation proposed by Lenhard et al. (1989). Detailed descriptions of these hydrogeologic units can be found in several publications (Rohay and McMahon, 1996; Swanson et al., 1999; Rohay, 2000; Serne et al., 2002; William et al., 2002). Wastewater and NAPL mixtures were discharged to the soil at three primary disposal facilities near the Z Plant (called the Plutonium Finishing Plant) at the Hanford Site. NAPL mixtures consist primarily of CT. Our efforts focus on discharge at the 216-Z-9 Trench. Here, approximately 4,090,000 L of wastewater, which included DNAPL CT, was discharged from July 1955 to June 1962 (Last and Rohay, 1993). Water and DNAPL discharge rates at the 216-Z-9 Trench are shown in Fig. 1. The area of the Trench was 18.3 m (length) by 9.1 m (width). Natural water infiltration rates from precipitation vary at the Hanford Site, depending upon soil texture and the degree of vegetation. A rate of 0.005 m/yr during the dry season from March to October and a rate of 0.075 m/yr during the wet season from November to February were previously estimated from natural infiltration rates at the Hanford CT Site (Kincaid et al., 1998; Hoitink et al., 2002). 3.2. Model simulation A simplified two-dimensional cross-section was selected to represent strongly-stratified layered features at the Hanford Site. The model domain is 662 m wide by 58.5 m deep. As

Fig. 1. Water and NAPL discharge rates for the 216-Z-9 Trench.

166

H. Yoon et al. / Journal of Contaminant Hydrology 90 (2007) 159–183

previously described, five stratigraphic layers are considered. However, unlike the Hanford Site, variations in layer thickness and slope are ignored. We focus on the importance of the NAPL and water spill scenario at the ground surface on vapor movement in the upper layers, and the effect of heterogeneity on the distribution and migration of DNAPL in a deep vadose zone. The two PlioPleistocene and the Hanford fine sand layers are fine-grained and have low permeabilities. They are important horizontal barriers for preventing vertical migration of water and NAPL. Heterogeneities in these layers are expected to significantly affect water and NAPL migration so they are considered. In contrast, the Hanford upper coarse and the Ringold Unit E layers are assumed homogeneous with an anisotropy ratio of 10 to capture the strongly-stratified feature. The model domain was discretized using 131 (x-direction) × 50 (z-direction) cells. The grid was refined vertically in both PPlc and PPlz layers, vertically near the top boundary, and horizontally near the NAPL and wastewater release area (i.e., source zone). Vertical grid spacings for the PPlc and PPlz layers were 0.5 m and 0.3 m, respectively, and vertical grid spacings near the top boundary ranged from 0.1 m to 0.5 m. The vertical grid spacing for the PPlz was chosen to be smaller than the corresponding air entry pressure in order to ensure reasonable capillary effects at the interface between the PPlz and the PPlc. The vertical grid spacing of the top boundary cells (0.1 m) was chosen in order to allow NAPL evaporation to occur by a diffusion process through the source zone. The horizontal spacing for the source zone was 1 m. Preliminary results suggested that the horizontal spreading of water and NAPL was limited to a short lateral distance surrounding the source zone. Hence, random fields for low permeability layers were generated

Fig. 2. Numerical domain, boundary conditions, and distribution of heterogeneous permeability in stratigraphic units (Case 1). White stands for mean permeability and black and gray colors stand for permeability an order of magnitude lower and higher than a mean value in each layer, respectively. The model domain is 662 m wide by 58.5 deep. A part of the domain is shown.

H. Yoon et al. / Journal of Contaminant Hydrology 90 (2007) 159–183

167

over only the distance of 86 m as shown in Fig. 2, while the rest of the domain was assumed to be homogeneous in each layer. The sequential indicator simulation model (SISIM) from the GSLIB package (Deutsch and Journel, 1998) was used to generate a two-dimensional random permeability field consisting of three discrete soil classes for the heterogeneous layers. The three indicator classes represent low, mean, and high permeabilities and the permeability of each class differs by an order of magnitude. The soil properties, correlation ranges, variances, and soil compositions for the soil classes are listed in Table 1. To describe continuous subunits in the heterogeneous layers, the random fields were generated by conditioning upon data that were arbitrarily selected, based upon Hanford data (Khaleel and Freeman, 1995; Rohay, 2000; Serne et al., 2002; William et al., 2002). A random field, which contains continuous lenses and composition of soil classes similar to input data, was chosen for each layer. Each layer is statistically independent with anisotropy ratios (kx /kz) of 5 and 10 for the heterogeneous and homogeneous layers, respectively. An anisotropy ratio of 10 was used for the Hanford Site in previous investigations (Piepho, 1996; Oostrom et al., 2003b) and recent field scale experiments have confirmed that is a reasonable anisotropy ratio (Yeh et al., 2005). The random fields used in this study are shown in Fig. 2. Boundary conditions used in the simulations are as follows (Fig. 2). For the top boundary, a constant atmospheric gas pressure (99,560 Pa) (except when barometric pumping was considered) and zero gas concentration were imposed. The value of the constant gas pressure was an average value of measured pressure fluctuations (Rohay, 1996). For the bottom boundary, no flux conditions (no gradient in pressure and concentration) were imposed. For the left- and right-boundaries, a time-varying Dirichlet condition (water pressure) for the water phase was derived from a one-dimensional simulation with natural water-infiltration event for 38 years, a free gradient boundary condition was imposed for gas phase, and a no flow condition was imposed for the NAPL phase. The gas pressure at the boundary was linearly extrapolated from two adjacent interior nodes by maintaining the pressure gradient (White and Oostrom, 2000). For NAPL components in the gas and aqueous phases, zero concentration was used. The domain was wide enough to avoid boundary effects on contaminant transport at the left- and right-boundaries. The NAPL and wastewater release area is 18 m long and is located at the center of the top boundary. A water table was placed approximately 8.0 m from the bottom boundary at the center of the domain and is inclined at a 0.005 slope from left to right. Initial conditions for the model domain were obtained by conducting a long-term simulation (i.e., 1000 years) using a constant water infiltration rate of 0.005 m/yr that is considered a Table 2 Chemical properties used in simulations Variable

CT −3 a

Density (kg m ) Viscosity (Pa·s)a Vapor pressure (Pa)a CT–contaminated water–air surface tension (N m− 1) CT–water equilibrium interfacial tension (N m− 1) CT–air equilibrium surface tension (N m− 1)b Henry's constant (Pa) Vapor diffusion coefficient (cm2 s− 1)

1426 1.11E− 3 10830 0.0608 0.0401 0.0207 1.3062 E + 8 0.083

a A value represent the property of mixture of 8.8% TBP, 14.7% DBBP, 2.9% lard oil, and 73.6% CT at 20 °C (Oostrom et al., 2003b). b A value was computed, assuming a zero spreading coefficient.

168

H. Yoon et al. / Journal of Contaminant Hydrology 90 (2007) 159–183

reasonable natural infiltration rate during the dry season at the Hanford Site (Oostrom et al., 2003b). Initial water saturations were very close to the residual water saturations for each soil (Table 1). Initially, there was no NAPL in the domain. Using the initial condition, simulations of NAPL migration were performed for 38 years, just prior to the start of SVE operations at the 216Z-9 Trench. Organic chemicals discharged at the Hanford Site were from 50% to 85% CT by volume. This mixture has different properties from pure phase CT. Measured properties of average mixture composition of 73.6% CT, 8.8% tributyl phosphate, 14.7% dibutyl butyl phosphonate, and 2.9% lard oil for density, vapor pressure, solubility, and viscosity (Oostrom et al., 2003b) were used for modeling. Other properties including interfacial tensions were assumed to be those for pure CT. Chemical and physical properties used in simulations are summarized in Table 2. The wastewater discharged was assumed to have water properties. No data are available regarding maximum residual NAPL saturation and maximum trapped NAPL saturation for the Hanford soils. Maximum residual NAPL saturation was reported to range from 5% to 15% for medium to fine sand materials (Wilkins et al., 1995; Van Geel and Roy, 2002; Oostrom et al., 2003a). Zalidis et al. (1998) observed that the residual NAPL saturation of a natural sandy–clay–loam soil was about 3% at residual water saturation. In this study, maximum residual NAPL saturation was set at 5% in order to simulate the effect of the formation of residual NAPL on the NAPL migration. Maximum trapped NAPL saturation has been reported in the literature, ranging from 10% to 25% for coarse to fine materials (Essaid et al., 1993; Oostrom and Lenhard, 2003; White et al., 2004); here the value was set between 10% and 25%, in inverse proportion to the permeabilities for each soil. Limited information is available regarding the CT sorption parameters at the Hanford Site. Linear sorption coefficients were chosen based upon CT sorption results with the Hanford sand and silty sand (Yonge et al., 1996) and estimated values using organic carbon contents at the Hanford Site (Truex et al., 2001). To simplify the model, the effect of temperature upon parameters used in simulations was not considered; however, in reality, temperature varies over time and would affect chemical properties such as vapor pressure and viscosity (Hoitink et al., 2002), especially near the ground surface. In

Table 3 Conditions used in simulation casesa Case Permeability fieldb 1 2 3 4 5 6 7

Vertical grid size of a cell adjacent Barometric to the top boundaryc pumping

Heterogeneous 10 cm Heterogeneous Uniform 5 cm

No

Artificial water discharge

Natural water infiltration (m/yr)

Initial condition

Yes

0.005–0.075 d

Dry

No No

0.30 e 0.30 e

Wet

Yes

a

Case 1 is the reference case. Other Cases use conditions defined in Case 1 except where conditions are described. Anisotropy ratios of 5 and 10 were used for heterogeneous and uniform fields, respectively. c This vertical grid size represents a relative thickness of boundary layer for NAPL evaporation to the atmosphere. d A rate of 0.005 m/yr from March to October and a rate of 0.075 m/yr from November to February were repeated. This cycle was repeated for the entire simulation period. e A water infiltration rate of 0.30 m/yr for 0.025 years was followed by no water infiltration for 0.075 year. This cycle was repeated for the entire simulation period. b

H. Yoon et al. / Journal of Contaminant Hydrology 90 (2007) 159–183

169

addition, thermal effects were also ignored due to the lack of information at the time of discharge. All phase partitioning was assumed to be at equilibrium, indicating that for conditions considered in this study the time scales for equilibrium partitioning are assumed to be smaller than those for advective transport. 3.3. Model case descriptions Table 3 summarizes the simulations that were performed for seven different cases to illustrate the effects of heterogeneity in the PPlz, NAPL evaporation to the atmosphere, barometric pumping, water discharge rates, and water saturation on the distribution and migration of NAPL. Case 1 is considered the base case and the permeability field is the one shown in Fig. 2. Water and NAPL discharge rates for Case 1 are divided into two stages. During Stage I, wastewater and NAPL discharges occurred over a 7-year period in the Trench as shown in Fig. 1. The NAPL and wastewater discharges were modeled as a continuous source (Piepho, 1996; Oostrom et al., 2003b). Total CT mass released during Stage I was 24,204 kg. No further NAPL or wastewater discharge occurred during Stage II, which lasted an additional 31 years. Natural water infiltration rates were maintained over the entire top boundary during Stage I and Stage II. Because the PPlz layer is believed to play a significant role in slowing NAPL migration to the water table, Case 2 is used to investigate the effects of heterogeneity in this layer on NAPL migration and distribution. For Case 2, six different permeability field realizations for the PPlz among 100 realizations were chosen, while all other layers were the same as in Case 1. Parameters used for the indicator simulations are listed in Table 1. The six realizations included two that have the most continuous horizontal distributions of the low permeability soil class of three indicator classes, and four that have the most discontinuous horizontal distributions of this soil class among 100 realizations. For Case 3, a uniform field for each layer with an anisotropy ratio of 10 was used. The geometric mean permeability for each layer in Case 1 was very similar to the average values in Case 3. Although full Monte Carlo analysis was not performed, we believe that the primary impacts of connectivity of the low permeability soil class on NAPL migration can be determined through these simulations because the six of 100 realizations chosen include the extreme cases of connectivity. Case 4 focuses on NAPL evaporation because this process can have an important impact upon NAPL mass loss to the atmosphere during NAPL release events at the ground surface. Large grid spacings (e.g., 1 m) at the top boundary without considering the NAPL evaporation to the atmosphere were used by previous investigations of the Hanford site (Piepho, 1996; Oostrom et al., 2003b). Without using an explicit model of NAPL evaporation to the atmosphere (e.g., stagnant boundary layer model or empirical mass transfer model from NAPL surface to the atmosphere), for most finite volume models the diffusive mass flux to the atmosphere at the top boundary is computed based on the concentration gradient between the center of the top boundary cell and the top boundary. Since the gas concentration was set to zero at the top boundary, large grid spacing at the top boundary is likely to under-predict the diffusive flux to the atmosphere because the length across which diffusion occurs in the porous media is artificially large. Therefore, the vertical grid size at the top boundary cells in Case 4 was refined to 5 cm, which is half of that in Case 1 in order to perform sensitivity analysis of the effect of the top grid spacing on NAPL evaporation to the atmosphere. Case 5 is used to investigate the effects of barometric pressure fluctuations on NAPL migration and distribution. At the Hanford Site, barometric pressure fluctuations averaged about 1.3% of the average value (Rohay, 1996), but no simulations reported in the literature have taken this into

170

H. Yoon et al. / Journal of Contaminant Hydrology 90 (2007) 159–183

account. These fluctuations were approximated with a square wave in which the cycle of high and low pressures within ± 0.65% of the average pressure was repeated every 0.02 years. Since the barometric fluctuation pattern causes the maximum time step for Case 5 to be much smaller than that for Case 1, and preliminary tests showed that the convergence rate was much slower in Case 5 than in Case 1, a reduced horizontal distance of 222 m including the center of the domain was used for Case 5 in order to save computational time. Comparison of results between the full horizontal domain of 662 m and a smaller domain of 222 m without barometric pressure fluctuations indicated that the effect of the domain size on NAPL migration was not significant.

Fig. 3. Snapshots of water saturation (color bar), total NAPL saturation, trapped NAPL saturation, and residual NAPL saturation (Case 1) at (a) 3.8 yr, (b) 6 yr, (c) 9 yr, and (d) 38 yr. Contour lines stand for NAPL saturation. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

H. Yoon et al. / Journal of Contaminant Hydrology 90 (2007) 159–183

171

The pressure fluctuation was handled by changing the top boundary condition. All other conditions are the same as those in Case 1. Cases 6 and 7 are used to investigate the effects of water infiltration rate on NAPL migration and distribution in the vadose zone with different initial water saturations. For these cases, artificial water discharge to the trench was not considered, and a higher water infiltration rate of 0.30 m/yr for 0.025 years was followed by no water infiltration for 0.075 years. This cycle was repeated for the entire simulation period. The DNAPL discharge scenario was the same as in Case 1. For Case 6 (dry case), the initial water distribution was the same as in Case 1. For Case 7 (wet case), the initial water distribution was obtained by conducting a long-term simulation using a continuous water infiltration rate of 0.075 m/yr. The wet condition had an initial water saturation 0.1 to 0.2 greater than the dry condition. 4. Results 4.1. Base case (Case 1) Contour profiles of water saturation, total NAPL saturation, residual NAPL saturation, and trapped NAPL saturation at different times for Case 1 are shown in Fig. 3. At 3.8 years, the water front already reached the groundwater table due to high water discharge rates, and the NAPL front reached the PPlz. At this time the bottom of the PPlz was almost completely saturated with water (up to 99.5%) because the top of the PPlc acts as a capillary barrier for the imbibing water phase. Since the top of the PPlz was unsaturated, NAPL infiltrated into the PPlz with slight spreading at the interface between the PPlz and the HFS. However, a short time later, the NAPL formed a pool and spread laterally until the pressure increased enough for the NAPL to displace the underlying water in the bottom part of the PPlz. After approximately 5.0 years, NAPL reached the bottom of the PPlz where the lowest permeability layer is located. It again formed a pool and spread laterally outward. This spreading resulted from higher horizontal permeability due to an anisotropy ratio of 5, high water saturation in the PPlz, and capillary action at the interface. Once the overlying pressure of the NAPL pool exceeded the entry pressure of the PPlc, the NAPL continued to migrate down. NAPL reached the groundwater table after 8.3 years. Comparison of total, residual, and trapped NAPL saturations in Fig. 3 indicates that at all times, most NAPL existed in the free form. Residual NAPL saturation was relatively small because water contents were high. Trapped NAPL saturation was also very small. It reached a maximum after approximately 6 years because large amounts of infiltrating water continued to migrate through the system after the main NAPL front advanced through the system. Trapped NAPL was primarily located in the low permeability layers (e.g., the PPlz and HFS). The highest trapped NAPL saturation at 6 years (= 0.040) was much lower than the theoretical maximum trapped NAPL saturation (= 0.25). After Stage I, water saturation started to decrease because only natural infiltration, which was very small, supplied water to the system. As a result, the residual NAPL saturation is greater at 9 years than at 6 years (Fig. 3 (b) and (c)). In contrast, a portion of the pore space that was occupied by trapped NAPL became exposed to air and the trapped NAPL saturation decreased. Eventually, trapped NAPL saturation decreased to near zero after 9 years. Cumulative mass fluxes with time are shown in Fig. 4. The total CT mass retained in the caliche layer (i.e., the difference between mass flux into the PPlz and mass flux out of the PPlc) reached quasi-steady state after approximately 15 years when the relative NAPL permeability was small due to the decreased water and NAPL saturation. The low permeability layer (PPlz) slowed the movement of NAPL downward. After 38 years, 62.4% of the total CT mass discharged was

172

H. Yoon et al. / Journal of Contaminant Hydrology 90 (2007) 159–183

retained in the vadose zone. Of the total CT mass discharged, 11.3% of this crossed the water table. Although the amount of water discharged was large and equilibrium partitioning between aqueous and NAPL phases is assumed, the amount of CT dissolved by advective water flow in the vadose zone that crossed the water table was 0.2% of the total CT mass discharged. Results in Fig. 4 also show that substantial mass was lost through volatilization. During Stage I, volatilized mass out of the domain increased linearly over time. This is mainly due to a diffusive flux to the atmosphere through the top boundary. Most of the diffusive flux occurred over the 30 m top boundary, which includes the source zone (18 m). The effect of the diffusive flux on NAPL migration will be discussed in the next section. The other mechanism of vapor movement is density-driven gas advection affected by the flow of water and NAPL. To clarify the development of advective gas flow during the early stage of water and NAPL discharge, contour profiles of vapor density and gas velocity vectors are shown in Fig. 5. As NAPL is discharged from the source zone, most CT mass up to 0.5 years is dissolved into water, sorbed to soil, and volatilized because of the low NAPL discharge rate and high CT vapor pressure. As a result, the vapor density beneath the source zone rapidly increased from 1.17 kg/m3 to above 1.6 kg/m3. This resulted in gradients in the total gas potential (i.e., the sum of the gas pressure and gravitational potentials), which induced advective gas flow downward. The induced vertical downward gas flow overcame the upward gas pressure gradient beneath the source zone, resulting in a downward advective gas flow into the subsurface through the top boundary. Similar advective gas flow from the top boundary was reported by Mendoza and Frind (1990b) where an immobile NAPL source zone was located 0.3 m below the ground surface. The induced gas flow changed the gas pressure distribution in the HUC layer, causing pressure gradients to develop in the lateral direction. Because the anisotropy ratio was 10 in the HUC layer, the water saturation increased at the bottom of the HUC layer, and the HCU layer was underlain by the lower permeability HFS layer, advective gas flow moved away from the center in the HUC as shown in Fig. 5. In addition, the upward gas flow toward the open top boundary created large circulation flows as previously observed (Mendoza and Frind, 1990b; Lenhard et al., 1995). The increased NAPL discharge rate after 0.5 year caused a NAPL front to form beneath the source zone and move downward. This enhanced the downward and lateral migration of high density CT vapor. The maximum vertical gas velocity at 0.7 year was 0.156 m/day, which is approximately

Fig. 4. Cumulative mass fluxes into layers, out of boundaries, and across water table (Case 1).

H. Yoon et al. / Journal of Contaminant Hydrology 90 (2007) 159–183

173

Fig. 5. Snapshots of vapor density (color bar) and gas velocity (vectors) at (a) 0.05 yr, (b) 0.1 yr, (c) 0.2 yr, and (d) 0.7 yr. Pure air density is 1.173 kg/m3. Maximum horizontal Darcy gas velocities (m/d) are (a) 0.076, (b) 0.072, (c) 0.077, and (d) 0.156. Arrow sizes represent the relative magnitude of the Darcy gas velocities. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

two times higher than that at 0.2 year. After Stage I, there are no further water and NAPL discharges and hence the top of the NAPL region was first removed due to density-driven gas flow. As shown in Fig. 6, advective gas flow was continuously induced through the center region of NAPL. This pattern of gas flow quickly removed the top of the NAPL region. This behavior continued until most of the NAPL mass in the top layer (HUC) was removed (Fig. 6 (d)). Although advective gas flux is the dominant mechanism for vapor movement, diffusion processes affected the gas concentration plume toward the top boundary. At 15 years, a 0.1 g/L contour of gas concentration was uniformly distributed approximately 2.5 m below the top boundary. This indicates that diffusion was fast enough to maintain a uniform gas concentration contour (0.1 g/L) a fixed distance from the zero gas concentration at the top boundary. After 38 years, about 11% of the total CT discharged moved away from the NAPL region via advective gas flow. About 25.7% of the total CT mass discharged moved out of the domain as vapor through the top boundary, while less than 0.01% left the domain through the left and right boundaries due

174

H. Yoon et al. / Journal of Contaminant Hydrology 90 (2007) 159–183

to the large domain size. The total mass remaining in the gas, aqueous, and sorbed phases after 38 years was 8.0% of the total CT mass discharged. 4.2. Effect of heterogeneity in the PPlz (Cases 1–3) Values of mass in the vadose zone, mass that crossed the water table, and mass that volatilized for Cases 1 through 7 after 38 years are shown in Fig. 7. Case 2 in Fig. 7 includes mass distributions for all six realizations. Only the results for Cases 1 through 3 are discussed in this

Fig. 6. Snapshots of total NAPL saturation (color bar), gas concentration (g/L, contour lines), and gas velocity (vectors) at (a) 1 yr, (b) 7 yr, (c) 8 yr, and (d) 15 yr. Saturated gas concentration is 0.65 g/L. The model domain is 662 m wide. Arrow sizes represent the relative magnitude of the Darcy gas velocities. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

H. Yoon et al. / Journal of Contaminant Hydrology 90 (2007) 159–183

175

section. The mass of CT retained in the vadose zone after 38 years in Case 2 realizations varied depending upon the presence of discontinuities in the low permeability soil class in the PPlz layer. The mass of CT retained in the vadose zone in Case 3 was similar to that in Case 1. Compared with Case 1 (Fig. 3 (d)), NAPL in the PPlc and PPlz layers in Case 2 had slightly different profiles of high NAPL saturation (results not shown). This is mainly due to different distributions of the low permeability soil class in the PPlz for the six realizations. For the six extreme cases of 100 realizations, the CT mass that reached the water table varied by a factor of two (Fig. 7). The two realizations with continuous lenses of the low permeability soil class had slower downward NAPL migration and retained higher CT mass, resulting in less CT mass crossing the water table. The four realizations with discontinuities in the low permeability soil class provided a preferential NAPL flow path, resulting in high CT mass crossing the water table. These four realizations had similar amounts of CT mass crossing the water table, as shown in Fig. 7. This result indicates that in strongly-stratified heterogeneous porous media, the connectivity of the low permeability soil class can play an important role in NAPL migration. For Case 3, each layer has a uniform permeability field with an anisotropy ratio of 10. For the PPlc and PPlz layers in Case 3, the vertical permeability was close to the effective vertical permeability of the three different soil classes in Case 1, where the anisotropy ratio was 5. NAPL mass distributions were similar between Cases 1 and 3 (Fig. 7). For Cases 1–3, volatilized CT mass exiting the domain was very similar (Fig. 7) because advective–dispersive transport of dense CT vapor occurred mainly in the top layer and in the upper region of the HFS. 4.3. Effect of NAPL evaporation (Cases 1 and 4) Cumulative mass fluxes in the gas phase with time through the top boundary for Cases 1 and 4 are shown in Fig. 8. The only difference between Cases 1 and 4 is the vertical grid spacing of the

Fig. 7. Comparison of CT mass in the vadose zone, CT mass across water table, volatilized CT mass out of the domain after 38 years. For Case 2, the box contains results from four realizations that have the most discontinuous horizontal distributions of the low permeability soil class of three indicator classes among 100 realizations and volatilized mass (triangle) includes results from the six realizations.

176

H. Yoon et al. / Journal of Contaminant Hydrology 90 (2007) 159–183

top boundary cells. The total CT mass retained in the vadose zone and total CT mass across the water table was significantly affected by the thickness of the boundary layer as shown in Fig. 7. After 38 years in Case 4, 54.3% of the total CT mass discharged was retained in the vadose zone; 42.7% flowed out of the model domain in the gas phase, and only 2.9% moved across the water table. During NAPL discharge (Stage I), the gas concentration at the center of the top boundary cell was near or at the saturation vapor pressure. Since the gas concentration was set to zero on the top boundary, the total mass flux of CT vapor to the atmosphere through the top boundary was approximately two times greater with the vertical grid size of 5 cm than with that of 10 cm at the top boundary (Fig. 8). As expected, most volatilized mass diffused out of the domain through the 30 m top boundary including the source zone during Stage I (Fig. 8); total mass flux to the atmosphere was almost the same as the diffusive flux during Stage I. Hereafter, the diffusive zone (30 m) refers to the 30 m top boundary including the source zone (18 m). As NAPL was removed from the upper region via volatilization during Stage II, gas concentrations adjacent to the top boundary decreased. As a result of decreased gas concentrations, the rate of mass loss to the atmosphere decreased sharply through the diffusive zone (30 m), while it decreased more gradually through the entire top boundary due to the length of the domain (662 m) (Fig. 8). The mass loss to the atmosphere decreased NAPL mass migrating downward through the PPlz and the PPlc layers and slowed vertical NAPL movement. Hence, a small amount of NAPL crossed the water table. 4.4. Effect of barometric pumping (Cases 1 and 5) Cumulative mass fluxes in the gas phase with time through the top boundary for Cases 1 and 5 are shown in Fig. 8. Barometric pumping slightly affected total mass flux through the top boundary. After 38 years, the CT mass distributions were similar for Cases 1 and 5 (Fig. 7).

Fig. 8. Comparison of cumulative volatilized mass fluxes out of the domain (Cases 1, 4–5). The mass fluxes through the source zone were calculated over the distance of 30 m surrounding the source zone.

H. Yoon et al. / Journal of Contaminant Hydrology 90 (2007) 159–183

177

During Stage I, CT vapor movement was primarily governed by density-driven gas advection. Vertical pressure gradients induced by barometric pumping were small due to the anisotropy ratio of 10 in the top layer; these gradients were overwhelmed by the gradients in the total gas potentials induced by dense CT vapor. As a result, NAPL distribution and gas concentrations in Case 5 were almost identical to those in Case 1 (results not shown). During Stage II, the total CT mass flux into the atmosphere for Case 5 was greater than that for Case 1, but after 38 years, the difference was subtle as shown in Fig. 8. After artificial water and NAPL discharges ceased, NAPL in the top of the domain was removed. As the gradients of the downward advective gas flow beneath the source zone gradually decreased, the barometric pumping had more impact on the vapor movement in the upper region, resulting in higher total CT mass flux into the atmosphere for Case 5 than for Case 1. However, because of the decreased gas concentrations and low vertical permeability in the upper region, the effect of barometric pumping on NAPL migration was very small in the continuous NAPL spill scenario considered in this study (Fig. 8). 4.5. Effect of water infiltration rate and initial water content (Cases 1, 6, and 7) Contour profiles of total NAPL saturation, residual NAPL saturation, and trapped NAPL saturation in Cases 6 and 7 after 38 years are shown in Fig. 9. Vertical NAPL movement was slower in Case 6. In Cases 6 and 7, more CT mass in the gas phase moved out of the domain through the top boundary, and in Case 7, less CT mass was retained in the vadose zone, compared to Case 1 (Fig. 7). This can be attributed to the low water saturation near the top boundary in Cases 6 and 7 because there was no artificial water discharge during Stage I. For a lower water saturation, a larger diffusive flux to the atmosphere near the diffusive zone (30 m) is expected due to the longer retention time of NAPL in the upper layer (HUC) and the smaller tortuosity in the gas phase.

Fig. 9. Snapshots of total NAPL saturation, residual NAPL saturation, and trapped NAPL saturation after 38 years: (a) Case 6 and (b) Case 7. Color bar stands for water saturation and contour lines stand for NAPL saturation. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

178

H. Yoon et al. / Journal of Contaminant Hydrology 90 (2007) 159–183

Residual NAPL mass in Cases 6 and 7 was smaller than that in Case 1 after 38 years, which can be attributed to large trapped NAPL mass for Case 6 and high water saturation for Case 7 (Fig. 9). In Case 6, only 0.2% of the total CT mass discharged moved across the water table (Fig. 7). This can be attributed to the NAPL relative permeability. In the Burdine relative permeability function (Eq. (10)), NAPL relative permeability increases with increasing water saturation at the same free NAPL saturation in unsaturated conditions because NAPL occupies large pore spaces where capillary forces are small. The average water saturation, NAPL saturation, and NAPL relative permeability over a 10 m horizontal distance in the PPlz (as indicated in Fig. 2) were computed for Cases 1 and 6 in Fig. 10. The NAPL front reached the top of the PPlz earlier in Case 1 than in Case 6 due to a higher water saturation in Case 1. In case 1, NAPL relative permeability was high until water saturation started to decrease at 8 years. After 8 years, NAPL relative permeability decreased with decreasing water saturation and relatively constant NAPL saturation.

Fig. 10. The change in (a) water saturation and NAPL saturation, and (b) NAPL relative permeability over time for Cases 1 and 6. The values are averaged over a 10 m horizontal distance in the PPlz as indicated in Fig. 2.

H. Yoon et al. / Journal of Contaminant Hydrology 90 (2007) 159–183

179

In contrast, NAPL relative permeability for Case 6 started to increase at 12 years when the NAPL front reached the PPlz, but the water front did not. In Case 6, the water saturation, NAPL saturation, and NAPL relative permeability were relatively constant over time between 17 years and 22 years due to the low water flow rate. After 22 years, NAPL relative permeability started to decrease with decreasing free NAPL saturation and increasing trapped NAPL saturation as the water front reached the PPlz. The formation of trapped NAPL was significantly affected by water infiltration rates for relatively dry initial conditions (Cases 1 and 6) and by initial water content for the same water infiltration scenario (Cases 6 and 7). As explained earlier, trapped NAPL mass in Case 1 was relatively small due to the high water discharge rate during Stage I. For Case 6, the water front reached the PPlz at about 22 years (Fig. 10 (a)) when water saturation was low and NAPL saturation was high, resulting in high trapped NAPL saturation as shown in Fig. 9. In contrast, for the wet condition (Case 7), NAPL occupied only the larger pore spaces at high water saturation when the water front swept downward, so the amount of trapped NAPL was low. 5. Discussion and summary Comparison of the migration and distribution of NAPL for Cases 1–3 demonstrated that the overall permeability and connectivity of low permeability soil in the PPlz and PPlc layers were important factors controlling the NAPL distribution. This conclusion is in agreement with results from a controlled field study in the vadose zone reported by Poulsen and Kueper (1992) that demonstrated that NAPL saturation is spatially variable and highly dependent upon small-scale spatial variations in hydraulic properties. In particular, simulation results using six realizations in Case 2 demonstrated that the connectivity of a low permeability soil in the lowest permeability layer (the PPlz) affected CT mass that crossed the water table by − 32% to 43%, compared to the base case (Case 1). However, variability of the NAPL distributions for six realizations was not significant, compared to when enhanced NAPL evaporation (Case 4) and water infiltration for initially dry conditions (Case 6) were considered. The CT mass that reached the water table in Cases 4 and 6 decreased by 74% and 98%, compared to Case 1, respectively. The formation of trapped NAPL in the low permeability layers was not important for all six realizations as with Case 1. In contrast to what would be expected in deep heterogeneous layers like the Hanford CT site, NAPL migration to the water table was governed by two factors: high water discharge rates during Stage I which enhanced NAPL migration downward to the PPlz and NAPL evaporation and gas advection which reduced downward NAPL mass migration into the lower layers. Comparison of Cases 1 and 4 showed that the vapor diffusive flux to the atmosphere had a dramatic impact upon NAPL migration during Stage I, when NAPL phase existed at the source zone, resulting in high CT mass flux to the atmosphere. This was due to the smaller vertical grid spacing of 5 cm or the smaller diffusion length in the top boundary in Case 4. This high diffusive flux to the atmosphere was also observed by a recent field experiment (Broholm et al., 2005) and several empirical models of NAPL evaporation rate from NAPL pools to the atmosphere (Huang, 1997). Since NAPL migration is very sensitive to the diffusive flux through the source zone which is proportional to the period of NAPL discharge and the size of the source zone, there is a strong need to develop release scenarios at the ground surface that include periodic NAPL discharges (spill rates and periods) and different source zone sizes. In addition, experimental work is needed to verify the rate of NAPL evaporation under various water and NAPL discharge scenarios. Comparison of Cases 1 and 5 showed that barometric pumping had a minimal impact upon NAPL migration (Case 5) during Stage I with a continuous water and NAPL spill. When the

180

H. Yoon et al. / Journal of Contaminant Hydrology 90 (2007) 159–183

NAPL phase existed at the source zone, high diffusive mass flux to the atmosphere occurred. Since high gradients in the total gas potential were induced inside the NAPL plume due to NAPL discharge at the ground surface, barometric pressure fluctuations were not enough to change the direction of the gas advective flow, resulting in minimal effect on the total mass flux to the atmosphere. Although barometric pumping only marginally affected the vapor movement under the continuous water and NAPL spill scenario considered in this study, barometric pumping may have greater impact on the advective mass flux to the atmosphere through the source zone when a periodic spill scenario is considered. Comparison of Cases 1 and 6 showed that NAPL flow was significantly affected by the high water discharge rates that resulted in high water saturation and high NAPL relative permeability during Stage I in Case 1. Although previous studies (Piepho, 1996; Oostrom et al., 2003b) addressed modeling DNAPL transport at Hanford, our study is unique because the effects of water saturation on NAPL movement are demonstrated. In contrast to Cases 1 and 7, results for Case 6 showed that water saturation had a great influence upon the trapped NAPL saturation. This was experimentally demonstrated by Oostrom and Lenhard (2003); CT was distributed at residual water saturation in a caliche layer (e.g., PPlc and PPlz layers) of a sand box, and this was followed by high water infiltration from a source zone on the top of the box. In contrast to what would be expected in semiarid regions like the Hanford Site, the formation of trapped NAPL in a water-wetting system was limited by a high water discharge rate that caused the vadose zone to be very wet before the NAPL front reached the low permeability layers. This high water saturation resulted in a continuous wetting phase providing infiltrating water with high water relative permeability during Stage I. It should be noted that for all cases, volatilization and diffusion (or evaporation) to the atmosphere were responsible for removing NAPL saturation in the upper region of the domain where residual NAPL is a high portion of total NAPL. This implies that in the field, volatile organic contaminants are likely to volatilize into the gas phase, resulting in smaller residual NAPL saturation than observed in the laboratory. These results indicate that in the field, residual NAPL saturation is strongly dependent on water saturation and NAPL volatility. In addition, due to the dependence of the form of NAPL upon soil moisture dynamics, prediction of NAPL distribution in the vadose zone requires initial water distribution and infiltration data in heterogeneous porous media. The development of a conceptual model for NAPL distribution and cost-effective remediation strategies at the Hanford CT Site has been complicated because a significant amount of CT was not accounted for by the current inventory model (DOE, 2003). Although site-specific conclusions cannot be drawn quantitatively, this study gives important insights into DNAPL distribution in a thick heterogeneous vadose zone similar to conditions at the Hanford Site. First, NAPL mass across the water table might be reduced significantly by NAPL evaporation to the atmosphere. In this study, 42.7% of the total CT mass discharged moved out of the model domain in the gas phase, with a large fraction moving to the atmosphere due to vapor diffusion (Case 4). Barometric pumping had a marginal effect (Case 5). Since the current inventory model at the Hanford CT Site includes only the barometric pumping effect for mass lost to the atmosphere, taking into account NAPL diffusion to the atmosphere and density-driven gas advection in a three-dimensional domain is essential to develop a conceptual model for the distribution and migration of NAPL, especially for NAPL spill events at the ground surface. Although the composition of NAPL mixtures was constant over time in this study, the mole fraction of CT in NAPL mixtures under high NAPL evaporation rates will decrease rapidly, causing NAPL evaporation rates to be reduced. Therefore, the boundary layer thickness, the source zone size, water and NAPL infiltration rates, and the change of NAPL composition over time must be considered to investigate NAPL fate after

H. Yoon et al. / Journal of Contaminant Hydrology 90 (2007) 159–183

181

release. Second, high water discharge rates enhanced and the presence of the low permeability PPlz layer slowed the vertical migration of NAPL. Local heterogeneities in the PPlz, in particular discontinuities of the low permeability soil class, may allow more NAPL to migrate through the PPlz layer, increasing the amount of NAPL mass transport across the water table. In strongly layered heterogeneous porous media in a deep vadose zone, the effect of local heterogeneities was enhanced by high water saturation due to water discharge events. Since there are three wastewater and CT discharge facilities where discharge rates and source zone areas are different, proper infiltration rates for different discharge facilities must be considered. Third, the formation of trapped NAPL is not likely to be a dominant process at the Hanford Site assuming a strongly waterwetting system, despite high water and NAPL discharge rates. Instead, NAPL likely exists in residual and free forms. However, under varying mixture compositions and interfacial tensions for a long term contamination, the effects of wettability change and hysteresis in capillary pressure– saturation relationships require further investigation. Acknowledgments Support for this work was provided by the U.S. Department of Energy–Environmental Science Management Program Project No. 70045. We thank Mark D. White for helping us to obtain and modify STOMP code, Mart Oostrom for the helpful remark during manuscript preparation, and Robert J. Lenhard and Mark L. Rockhold for their assistance in obtaining articles and reports. We also acknowledge the effort of the Associate Editor and Reviewers, whose comments have led to significant improvement of our manuscript. References Broholm, M.M., Christophersen, M., Maier, U., Stenby, E.H., Höhener, P., Kjeldsen, P., 2005. Compositional evolution of the emplaced fuel source in the vadose zone field experiment at airbase Værløse, Denmark. Environ. Sci. Technol. 39, 8251–8263. Brooks, R.H., Corey, A.T., 1964. Hydraulic properties of porous media. Hydrology Papers, vol. 3. Colorado State University, Fort Collins, CO. Choi, J.-W., Tillman, F.D., Smith, J.A., 2002. Relative importance of gas-phase diffusive and advective trichloroethene fluxes in the unsaturated zone under natural conditions. Environ. Sci. Technol. 36, 3157–3164. Deutsch, C.V., Journel, A.G., 1998. GSLIB Geostatistical Software Library and User's Guide. Oxford University Press, NY. DOE, 2001. A national roadmap for vadose zone science and technology. Understanding, Monitoring, and Predicting Contaminant Fate and Transport in the Unsaturated Zone. Rev. 0. DOE/ID-10871. DOE, 2003. Alternatives for Carbon Tetrachloride (CT) Source Term Location. DE-RA26-03NT41735. National Energy Technology Laboratory. Ellerd, M.G., Massmann, J.W., Schwaegler, D.P., Rohay, V.J., 1999. Enhancements for passive vapor extraction: the Hanford study. Ground Water 37 (3), 427–437. Essaid, H.I., Herkelrath, W.N., Hess, K.M., 1993. Simulation of fluid distributions observed at a crude-oil spill site incorporating hysteresis, oil entrapment, and spatial variability of hydraulic properties. Water Resour. Res. 29 (6), 1753–1770. Hofstee, C., Dane, J.H., Hill, W.E., 1997. Three-fluid retention in porous media involving water, PCE, and air. J. Contam. Hydrol. 25, 235–247. Hofstee, C., Oostrom, M., Dane, J.H., Walker, R.W., 1998. Infiltration and redistribution of perchloroethylene in partially saturated, stratified porous media. J. Contam. Hydrol. 34, 293–313. Hoitink, D.J., Ramsdell, J.V., Burk, K.W., Shaw, W.J., 2002. Hanford Site Climatological Data Summary 2001 With Historical Data. PNNL-13859. PNNL, Richland, WA. Huang, C.H., 1997. Pasquill's influence: on the evaporation from various liquids into the atmosphere. J. Appl. Meteorol. 36, 1021–1026. Jellali, S., Benremita, H., Muntzer, P., Razakarisoa, O., Schafer, G., 2003. A large-scale experiment on mass transfer of trichloroethylene from the unsaturated zone of a sandy aquifer to its interfaces. J. Contam. Hydrol. 60 (1–2), 31–53.

182

H. Yoon et al. / Journal of Contaminant Hydrology 90 (2007) 159–183

Kaluarachchi, J.J., Parker, J.C., 1992. Multiphase flow with a simplified model for oil entrapment. Transp. Porous Media 7 (1), 1–14. Keller, A.A., Chen, M.J., 2003. Effect of spreading coefficient on three-phase relative permeability of nonaqueous phase liquids. Water Resour. Res. 39 (10), 1288. doi:10.1029/2003WR002071. Khaleel, R., Freeman, E.J., 1995. Variability and scaling of hydraulic properties for 200 area soils, Hanford Site. Report WHC-EP-0883. Westinghouse Hanford Company, Richland, WA. Kincaid, C.T., Bergeron, M.P., Cole, C.R., Freshley, M.D., Hassig, N.L., Johnson, V.G., Kaplan, D.I., Serne, R.J., Streile, G.P., Strenge, D.L., Thorne, P.D., Vail, L.W., Whyatt, G.A., Wurstner, S.K., 1998. Composite Analysis for Low-Lever Waste Disposal in the 200 area Plateau of the Hanford Site. PNNL-11800. PNNL, Richland, WA. Last, G.V., Rohay, V.J., 1993. Refined Conceptual Model for the Volatile Organic Compounds — Arid Integrated Demonstration and 200 West Area Carbon Tetrachloride Expedited Response Action. PNL-8597. PNNL, Richland, WA. Lenhard, R.J., Parker, J.C., Mishra, S., 1989. On the correspondence between Brooks–Corey and Vangenuchten models. J. Irrig. Drain. Eng. ASCE 115 (4), 744–751. Lenhard, R.J., Oostrom, M., Simmons, C.S., White, M.D., 1995. Investigation of density-dependent gas advection of trichloroethylene — experiment and a model validation exercise. J. Contam. Hydrol. 19 (1), 47–67. Lenhard, R.J., Oostrom, M., Dane, J.H., 2004. A constitutive model for air–NAPL–water flow in the vadose zone accounting for residual NAPL in strongly water-wet porous media. J. Contam. Hydrol. 73 (1–4), 281–304. Mendoza, C.A., Frind, E.O., 1990a. Advective–dispersive transport of dense organic vapors in the unsaturated zone. 1. Model development. Water Resour. Res. 26 (3), 379–387. Mendoza, C.A., Frind, E.O., 1990b. Advective–dispersive transport of dense organic vapors in the unsaturated zone. 2. Sensitivity analysis. Water Resour. Res. 26 (3), 388–398. Millington, R., Quirk, J., 1961. Permeability of porous solids. Trans. Faraday Soc. 57, 1200–1207. Neeper, D.A., 2002. Investigation of the vadose zone using barometric pressure cycles. J. Contam. Hydrol. 54 (1–2), 135–162. Oostrom, M., Lenhard, R.J., 1998. Comparison of relative permeability–saturation–pressure parametric models for infiltration and redistribution of a light nonaqueous-phase liquid in a sandy porous media. Adv. Water Resour. 21, 145–157. Oostrom, M., Lenhard, R.J., 2003. Carbon tetrachloride flow behavior in unsaturated Hanford caliche material: an investigation of residual nonaqueous phase liquids. Vadose Zone J. 2, 25–33. Oostrom, M., Hofstee, C., Lenhard, R.J., Wietsma, T.W., 2003a. Flow behavior and residual saturation formation of liquid carbon tetrachloride in unsaturated heterogeneous porous media. J. Contam. Hydrol. 64, 93–112. Oostrom, M., Rockhold, M.L., Thorne, P.D., Last, G.V., Truex, M.J., 2003b. Three Dimensional Modeling of DNAPL Movement and Redistribution in the Subsurface of the 216-Z-9 Trench at the Hanford Site. PNNL, Richland, WA. Oostrom, M., Lenhard, R.J., van Geel, P., White, M.D., 2005. Comparison of models describing residual NAPL formation in the vadose zone. Vadose Zone J. 4, 163–174. Pankow, J.F., Cherry, J.A., 1996. Dense Chlorinated Solvents and Other DNAPLs in Groundwater. Waterloo Press, Waterloo, Ontario, Canada. Piepho, M.G., 1996. Numerical Analysis of Carbon Tetrachloride Movement in the Saturated and Unsaturated Zones in the 200 West area, Hanford Site. BHI-00459 Rev. 0. Bechtel Hanford Inc., Richland, WA. Poston, T.M., Hanf, R.W., Dirkes, R.L., Morasch, L.F., 2004. Hanford Site Environmental Report for Calendar Year 2003. PNNL-14687. PNNL, Richland, WA. Poulsen, M.M., Kueper, B.H., 1992. A field experiment to study the behavior of tetrachloroethylene in unsaturated porous media. Environ. Sci. Technol. 26, 889–895. Rohay, V.J., 1996. Field Tests of Passive Soil Vapor Extraction Systems at the Hanford Site, Washington. BHI-00766 Rev. 0. Bechtel Hanford Inc., Richland, WA. Rohay, V.J., 2000. Performance evaluation report for soil vapor extraction operations at the carbon tetrachloride site, February 1992–September 1999. BHI-00720 Rev. 4. Bechtel Hanford Inc., Richland, WA. Rohay, V.J., McMahon, W.J., 1996. Airflow Modeling Report for Vapor Extraction Operations at the 200-ZP-2 Operable Unit (Carbon Tetrachloride Expedited Response Action). BHI-00882. Rev. 0. Bechtel Hanford Inc., Richland, WA. Serne, R.J., Mitroshkov, A.V., Serne, J.N., Bjornstad, B.N., LeGore, V.L., Last, G.V., Schaef, H.T., O'Hara, M.J., Smith, S.C., Williams, B.A., Brown, C.W., Lindenmeier, C.W., Lanigan, D.C., Parker, K.E., Zachara, J.M., Horton, D.G., Kutnyakov, I.V., Burke, D.B., Clayton, R.E., 2002. Characterization of Vadose Zone Sediment: Uncontaminated RCRA Borehole Core Samples and Composite Samples. PNNL-13757. PNNL, Richland, WA. Sleep, B.E., Sykes, J.F., 1989. Modeling the transport of volatile organics in variably saturated media. Water Resour. Res. 25, 81–92. Sleep, B.E., Sykes, J.F., 1993. Compositional simulation of groundwater contamination by organic compounds. 2. Model applications. Water Resour. Res. 29, 1709–1718.

H. Yoon et al. / Journal of Contaminant Hydrology 90 (2007) 159–183

183

Swanson, L.C., Rohay, V.J., Faurote, J.M., 1999. Hydrogeologic Conceptual Model for the Carbon Tetrachloride and Uranium/Technetium Plumes in the 200 West Area: 1994 through 1999 Update. BHI-01311 Rev. 0. Bechtel Hanford Inc., Richland, WA. Thomson, N.R., Sykes, J.F., Van Vliet, D., 1997. A numerical investigation into factors affecting gas and aqueous phase plumes in the subsurface. J. Contam. Hydrol. 28 (1–2), 39–70. Tillman, F.D., Choi, J.-W., Smith, J.A., 2003. A comparison of two methods for estimating the total flux of volatile organic compounds from the subsurface under natural conditions. Water Resour. Res. 39 (10), 1284–1294. Truex, M.J., Murray, C.J., Cole, C.R., Cameron, R.J., Johnson, M.D., Skeen, R.S., Johnson, C.D., 2001. Assessment of Carbon Tetrachloride Groundwater Transport in Support of the Hanford Carbon Tetrachloride Innovative Technology Demonstration Program. PNNL-13560. PNNL, Richland, WA. Van Geel, P.J., Roy, S.D., 2002. A proposed model to include a residual NAPL saturation in a hysteretic capillary pressure–saturation relationship. J. Contam. Hydrol. 58, 79–110. Van Geel, P.J., Sykes, J.F., 1997. The importance of fluid entrapment, saturation hysteresis and residual saturations on the distribution of a lighter-than-water non-aqueous phase liquid in a variably saturated sand medium. J. Contam. Hydrol. 25, 249–270. White, M.D., Oostrom, M., 2000. STOMP, Subsurface transport over multiple phases. Version 2.0. Theory Guide. PNNL12030. Pacific Northwest National Lab., Richland, WA. White, M.D., Oostrom, M., 2004. STOMP, Subsurface transport over multiple phases. Version 3.1. User's Guide. PNNL14478. Pacific Northwest National Lab., Richland, WA. White, M.D., Oostrom, M., Lenhard, R.J., 1995. Modeling fluid-flow and transport in variably saturated porous-media with the stomp simulator 1. Nonvolatile 3-phase model description. Adv. Water Resour. 18 (6), 353–364. White, M.D., Oostrom, M., Lenhard, R.J., 2004. A practical model for nonvolatile NAPL-aqueous flow in variably saturated porous media, distinguishing mobile, residual, and entrapped NAPL. Ground Water 42 (5), 734–746. Wilkins, M.D., Abriola, L.M., Pennell, K.D., 1995. An experimental investigation of rate-limited nonaqueous phase liquid volatilization in unsaturated porous media: Steady state mass transfer. Water Resour. Res. 31 (9), 2159–2172. William, B.A., Bjornstad, B.N., Schalla, R., Webber, W.D., 2002. Revised Hydrogeology for the Suprabasalt Aquifer System, 200-West Area and Vicinity, Hanford Site, Washington. PNNL-13858. PNNL, Richland, WA. Wipfler, E.L., van der Zee, S.E.A.T.M., 2001. A set of constitutive relationships accounting for residual NAPL in the unsaturated zone. J. Contam. Hydrol. 50, 53–77. Yeh, T.-C., Ye, M., Khaleel, R., 2005. Estimation of effective unsaturated hydraulic conductivity tensor using spatial moments of observed moisture plume. Water Resour. Res. 41, W03014. doi:10.1029/2004WR003736. Yonge, D., Hossain, A., Cameron, R., Ford, H., Storey, C., 1996. Hanford Soil Partitioning and Vapor Extraction Study. BHI-00861 Rev. 0. Bechtel Hanford Inc., Richland, WA. Yoon, H., Valocchi, A.J., Werth, C.J., 2003. Modeling the influence of water content on soil vapor extraction. Vadose Zone J. 2, 368–381. Zalidis, G.C., Wallace, R.B., Voice, T.C., 1998. Influence of initial water saturation on the residual saturation of an organic liquid in the vadose zone. Water Resour. Manag. 12 (2), 82–93. Zhu, J., Sykes, J.F., 2000. Stochastic simulations of NAPL mass transport in variably saturated heterogeneous porous media. Transp. Porous Media 39 (3), 289–314.