Parameters that control the cleanup of fractured permeable aquifers

Available online at www.sciencedirect.com

Journal of Contaminant Hydrology 96 (2008) 128 – 149 www.elsevier.com/locate/jconhyd

Parameters that control the cleanup of fractured permeable aquifers Hillel Rubin a,⁎, Sharon Yaniv a , Martin Spiller b , Jürgen Köngeter b b

a Faculty of Civil and Environmental Engineering, Technion, Israel Institute of Technology, Haifa 32000, Israel Institute of Hydraulic Engineering and Water Resources Management, Aachen University (RWTH Aachen), D-52056 Aachen, Germany

Received 14 September 2004; received in revised form 15 October 2007; accepted 19 October 2007 Available online 13 November 2007

Abstract This study develops a modeling approach for simulating and evaluating entrapped light nonaqueous-phase liquid (light NAPL– LNAPL) dissolution and transport of the solute in a fractured permeable aquifer (FPA). The term FPA refers to an aquifer made of porous blocks of high permeability that embed fractures. The fracture network is part of the domain characterized by high permeability and negligible storage. Previous studies show that sandstone aquifers often represent FPAs. The basic model developed in this study is a two-dimensional (2-D) model of permeable blocks that embed oblique equidistant fractures with constant aperture and orientation. According to this model, two major parameters govern NAPL dissolution and transport of the solute. These parameters are: 1) the dimensionless interphase mass transfer coefficient, Kf0, and 2) the mobility number, NM0. These parameters represent measures of heterogeneity affecting flow, NAPL dissolution, and transport of the solute in the domain. The parameter Kf0 refers to the rate at which organic mass is transferred from the NAPL into the water phase. The parameter NM0 represents the ratio of flow through the porous blocks to flow through the fracture network in regions free of entrapped NAPL. It also provides a measure of groundwater flow bypassing regions contaminated by entrapped NAPL. In regions contaminated by entrapped NAPL our simulations have often indicated very low permeability of the porous blocks, enabling a significant increase of the fracture flow at the expense of the permeable block flow. Two types of constitutive relationships also affect the rate of FPA cleanup: 1) the relationship between the saturation of the entrapped NAPL and the permeability of the porous blocks, and 2) the relationships representing effects of the entrapped NAPL saturation and the permeable block flow velocity on rates of interphase mass transfer. This study provides basic tools for evaluating the characteristics of pump-and-treat cleanup of FPAs by referring to sets of parameters and constitutive relationships typical of FPAs. The numerical simulations carried out in this study show that at high initial saturation of the entrapped NAPL, during initial stages of the FPA cleanup the contaminant concentration increases, but later it decreases. This phenomenon originates from significant groundwater bypassing the NAPL entrapped in the permeable blocks via the fracture network. © 2007 Elsevier B.V. All rights reserved. Keywords: Fractured porous formation; Fractured porous media; Fractured permeable formation (FPF); Fractured permeable aquifer (FPA); Aquifer cleanup; NAPL contamination; Pump-and-treat

1. Introduction

⁎ Corresponding author. Tel.: +972 4 829 2306; fax: +972 4 822 898. E-mail address: [email protected] (H. Rubin). 0169-7722/$ - see front matter © 2007 Elsevier B.V. All rights reserved. doi:10.1016/j.jconhyd.2007.10.007

This study originates from field observations and laboratory measurements (Rubin and Braester, 2000) performed in a part of the Coastal Plain Aquifer (CPA) in Israel. Entrapped kerosene, which is a light nonaqueous-

H. Rubin et al. / Journal of Contaminant Hydrology 96 (2008) 128–149

Nomenclature B C Cav Cb Cnv C⁎ ⁎ Cav Cb⁎ Cs⁎ d D i imax ibl ifr J J0 j ja jc Kb KB Kb0 Kf Kf0 Kt Kt0 kf krn krt krw krwB kα m Mf NM NM0 n nc np Qb Qb0 Qf Qf0 QR Qt qb qB

distance between adjacent fracture intersections (L) normalized solute concentration in the fracture flow normalized flux average solute concentration in the cross-section normalized solute concentration in the permeable block flow equilibrium volumetric concentration of solute solute concentration in the fracture flow (ML− 3) flux average solute concentration in the cross-section (ML− 3) solute concentration in the permeable block flow (ML− 3) equilibrium concentration of solute (ML− 3) grain diameter (L) free liquid diffusivity (L2/T) number of the vertical grid point total number of vertical grid points number of the permeable block number of the fracture segment hydraulic gradient hydraulic gradient in regions free of entrapped NAPL number of the longitudinal grid point number of the grid point for the calculation of Cav number of the downstream-end longitudinal grid point hydraulic conductivity of permeable blocks (LT− 1) value of Kb at the fracture segment (LT− 1) value of Kb in regions free of entrapped NAPL (LT− 1) dimensionless interphase mass transfer coefficient initial value of Kf average cross-sectional hydraulic conductivity (LT− 1) value of Kt in regions free of entrapped NAPL (LT− 1) lumped mass transfer coefficient (T− 1) relative NAPL permeability cross-sectional relative water permeability relative water permeability value of krw at the fracture segment number of the fracture segment nodal point for calculating Cav number of the time step mobility of the fracture segment (L2T− 1) mobility number value of NM in regions free of entrapped NAPL power coefficient used for permeability calculations number of contaminated sections coefficient for calculating permeability rate of permeable block flow (L2T− 1) value of Qb in regions free of entrapped NAPL (L2T− 1) rate of fracture segment flow (L2T− 1) value of Qf in regions free of entrapped NAPL (L2T− 1) ratio of Qb to Qf in a vertical cross-section total flow rate flowing through the subdomain cross-section (L2T− 1) specific discharge of the permeable block flow (LT− 1) value of qb entering the fracture segment (LT− 1)

129

130

H. Rubin et al. / Journal of Contaminant Hydrology 96 (2008) 128–149

qb0 value of qb in regions free of entrapped NAPL (LT− 1) qbx specific discharge in the longitudinal direction (LT− 1) qby specific discharge in the vertical direction (LT− 1) qben specific discharge entering the fracture elementary volume (LT− 1) qbex specific discharge leaving the fracture elementary volume (LT− 1) qr ratio of qb to qb0 qt average cross-sectional specific discharge (LT− 1) Re Reynolds number Sm irreducible water saturation Sn saturation of entrapped NAPL Sn0 initial value of Sn Snav average value of Sn in the subdomain cross-section Sw water saturation Sw0 initial water saturation Sw⁎ effective water saturation Sh Sherwood number T dimensionless time needed to carry out the FPA cleanup t dimensionless time t⁎ time (T) Vb interstitial flow velocity in the permeable block (LT− 1) Vb0 value of Vb in regions free of entrapped NAPL (LT− 1) vr ratio of Vb to Vb0 vr0 initial value of vr W dimensionless volume of water that should be treated x dimensionless longitudinal coordinate x⁎ longitudinal coordinate (L) y dimensionless vertical coordinate y⁎ vertical coordinate (L) α integer number βi (i = 0–2) coefficients used for calculating Sh and Kf Δt dimensionless time step Δx dimensionless longitudinal interval Δy dimensionless vertical interval θ orientation angle of fracture segments λ pore size distribution index ρn density of NAPL (ML− 3) ϕb porosity of the porous block

phase liquid (light NAPL–LNAPL), contaminates this aquifer. In 1983, people in the area complained that their tap water tasted of kerosene. In 1986, field studies (Kanfi, 1986) showed the existence of a kerosene layer up to 70 cm thick floating on top of the water table. At that time the water table was at an elevation of approximately 2 m below sea level. As a coastal aquifer, it was subject to risks of saltwater intrusion. In winter 1991/2 intensive natural recharge raised the water table to an elevation of 4 m above the sea level. The raised water table led to an almost complete disappearance of the formerly floating kerosene

layer, and to the entrapment of the kerosene in the top layers of the aquifer. The lithology shows that the water table is 40 m below the soil surface, and that the aquifer is a sandstone formation with minor quantities of sandy silt and sandy clay lenses. Pumping tests (Averbach, 1988) indicated that the aquifer displays high transmissivity (about 4000 m2/d) and hydraulic conductivity (about 100 m/d). However, the permeability of core samples was 1 to 2 orders of magnitude smaller than the magnitude suggested by the pumping tests. This information suggests

H. Rubin et al. / Journal of Contaminant Hydrology 96 (2008) 128–149

131

Fig. 1. Conceptualization of NAPL dissolution in a fractured permeable aquifer.

that the formation is a fractured porous medium, even though the pumping test analysis did not require reference to dual-porosity (Barenblatt et al., 1960). To describe such a formation, Birkhölzer et al. (1993a,b) suggested the term “Fractured Permeable Formation” (FPF), which refers to a fractured formation with significant block flow, and provided quantitative meaning to this term by applying the definition of the mobility number, as shown later in this paper. Here we adopt the term fractured permeable aquifers (FPAs) to refer to aquifers that comprise of FPFs. Tests of formation cores (Rubin and Braester, 2000) indicated the presence of entrapped kerosene in the permeable blocks of the top layers of the CPA's contaminated area. As a result, the water authorities are searching for various ways to contain or remediate the contaminated aquifer. Contemporary computer codes, as described in the literature (Miller et al., 1998) and in the webpages of commercial firms, can simulate virtually all flow and transport possibilities typical of dual-porosity formations. However, the present study, which concentrates on preferential flow and its changes during aquifer cleanup, may provide important information for optimizing aquifer cleanup operations. FPAs are aquifers that display significant transmittance through both a primary porosity of the porous blocks and a secondary porosity of the fracture network. In such aquifers the secondary porosity, which typically amounts to no more than a few percent of the total porosity, can substantially contribute to the effective conductivity (Pettijohn et al., 1987) of the aquifer. In FPAs free of entrapped NAPL, the mobility number represents the ratio of flow through the permeable block to the flow through

the fracture network. In FPAs reviewed by Rubin et al. (2001), the mobility number of FPAs typically is associated with values in the range of 0.1–25. Such aquifers are composed of fractured sandstone, like those found in Germany (Andres and Georgotas, 1978) and along the Coastal Plain of Israel (Kanfi, 1986, Gvirtzman et al., 1988). With respect to the CPA in particular, Gvirtzman et al. (1988) documented evidence of slow flow through the porous matrix and fast flow through the fractures using environmental tritium as a natural tracer. They found that approximately 25% of the total rain-flow infiltration occurred through the fissures, and 75% occurred through the permeable blocks. According to the U.S. National Research Council (NCR) (1994), NAPL contamination in groundwater comprises a major threat to satisfactory supply of drinking water. Therefore, measures of cleanup and remediation should be quantified to allow for water use according to engineering and economic feasibility. Within the context of this paper, quantitative measures are needed for the effective cleanup of FPAs. The studies of Schwille (1984, 1988) provide basic information concerning NAPL migration in porous and fractured media. Various field studies (e.g., Kanfi, 1986, Rubin and Braester, 2000) provide basic information about NAPL flow and entrapment in the part of the CPA that is a FPA. Studies of NAPL entrapment in homogeneous and heterogeneous porous media (e.g., Powers et al., 1991, 1992, 1994; Chatiz et al., 1983; Mackay et al., 1985; Mackay and Cherry, 1989; Mercer and Cohen, 1990; Unger et al., 1998; Frind et al., 1999; Zhang and Brusseau, 1999; Brusseau et al., 2000a,b, Nambi and Powers, 2000, 2003) are also helpful in evaluating NAPL contamination and FPA cleanup.

132

H. Rubin et al. / Journal of Contaminant Hydrology 96 (2008) 128–149

Fig. 2. The basic conceptual model: (a) macroscale; (b) subdomain; (c) elementary fracture volume.

Theoretical and practical case studies have shown that traditional pump-and-treat methods are often inefficient for removing residual NAPL (Mackay and Cherry, 1989, Powers et al., 1991, U.S. National Research Council, 1994). Nevertheless, water authorities continue using pump-and-treat methods for aquifer restoration purposes

(U.S. National Research Council, 1994, U.S. Environmental Protection Agency, 1996, 2001). Causes contributing to the inefficiencies of pump-and-treat cleanup may include nonequilibrium interphase partitioning at high water velocities or with decreasing interfacial contact area between the entrapped NAPL and the

H. Rubin et al. / Journal of Contaminant Hydrology 96 (2008) 128–149

flowing groundwater (Miller et al., 1990, Powers et al., 1992, 1994, Imhoff et al., 1994), and dilution and bypassing effects resulting from heterogeneous permeability and NAPL distribution (Anderson et al., 1992a, Brusseau 1992). Various versions of pump-and-treat involve the use of hydraulic barriers in an effort to avoid the contaminants being spread into regions of fresh groundwater. Various studies have developed simplified models for assessing the efficiency of pump-and-treat methods (Anderson et al., 1992b, Hatfield and Stauffer, 1993, Zaidel and Russo, 1993, Rubin et al., 1997, 2001, 2004). This paper extends the approach of Rubin et al. (1997, 2001, 2004), but considers how the reduced block permeability and porosity resulting from residual entrapped NAPL may affect the cleanup. Because of the NAPL dissolution, variations of the permeable block conductivity change the preferential flow, and in particular, the division of the aquifer discharge between the blocks and the fracture network. This study seeks to identify the specific characteristics of water flow, NAPL dissolution, and solute transport that arise in pump-and-treat remediation of FPAs, and in the long-term dissolution of NAPLs in such aquifers. We intend to develop a simplified method for quantifying the effects of aqueous flow bypassing the entrapped NAPL ganglia and the rate of interphase mass transfer. The effects of reduced rates of interphase mass transfer, as well as those due to permeability and porosity changes associated with the dissolution of the entrapped NAPL during FPA cleanup are emphasized. 2. The conceptual model Fig. 1 provides a schematic of a domain that comprises part of a FPA contaminated by entrapped NAPL. This domain represents a heterogeneous porous medium with particular features. The fracture network is part of the domain characterized by high permeability and negligible storage. Because of the negligible storage of the fracture network, quantities of NAPL entrapped within the fractures are very small. Groundwater flowing through the fracture network therefore flushes out the small NAPL quantities entrapped in the fractures quickly. It should be noted that in many fractured rock systems, flow is dominated by fracture flow and thus these sites are very difficult to remediate. In the case of FPAs, however, significant quantities of NAPL are trapped in the porous block formation and the NAPL trapped in the fractures can be neglected. This study applies the idealized model of Rubin et al. (1997, 2001, 2004) and extends its uses. Fig. 2 describes

133

Table 1 Assumptions applied in this study Assumption Type of assumption number 1 2 3 4 5 6 7 8 9

Simulated domain is two-dimensional (2-D). The fracture network incorporates two sets of equidistant fractures with constant aperture and orientation. Storage of the fracture network is negligible. Only porous blocks incorporate entrapped NAPL. Effects of solute diffusion are negligible. Solute dispersion is attributed only to mixing between the permeable block and fracture flows. Blocks are made of a homogeneous permeable matrix. Flow in the permeable blocks is mainly in the longitudinal x⁎ direction. Total discharge flowing through the matrix and the fractures is kept constant through the simulations; but the division of the discharge between the blocks and the fractures is subject to changes due to the dissolution of the entrapped NAPL.

this model, which is geometrically similar to models applied in previous studies referring to solute transport in fractured porous media (e.g., Berkowitz et al., 1988, Birkhölzer et al., 1993a,b, Rubin and Buddemeier, 1996). The model represents a two-dimensional FPA, incorporating homogeneous and isotropic permeable blocks. These blocks embed two sets of parallel, equidistant fractures with constant aperture and orientation. The distance between two adjacent fracture intersections measured along the fractures is B, and the angle of fracture orientation, that is, the angle between the fractures and the general flow direction, is θ. Part of the formation contains an immobile residual entrapped NAPL. As implied in the discussion above, initially we ignore the NAPL entrapped within the fractures and consider the NAPL entrapped only in the permeable blocks. The study concerns the pump-and-treat cleanup procedure whereby a uniform horizontal flow enters the contaminated part of the domain. Two strips of injection and pumping wells can yield such conditions, provided all wells pump at equal rates (Powers et al., 1991). Following the short time period of aquifer flushing, the domain arrives at “apparent steady-state conditions” (Rubin et al., 1997). Simulations of Rubin et al. (2001) show that for purposes of aquifer cleanup calculations, it is possible to neglect the aquifer flushing effects on the aquifer cleanup because of the short duration of the flushing time period. However, this study does not require that we neglect the effects of the flushing time period. During the apparent steady-state, the flowing groundwater gradually dissolves the entrapped NAPL,

134

H. Rubin et al. / Journal of Contaminant Hydrology 96 (2008) 128–149

decreasing its saturation over an extended period. In a homogeneous porous medium, the entrapped NAPL's saturation may be nonuniform (e.g., Geller and Hunt, 1993). In a heterogeneous porous medium, the heterogeneity leads to a nonuniform distribution of the entrapped NAPL, which may significantly affect the nonuniform NAPL dissolution rates. In such cases parts of the water flow bypass the entrapped NAPL ganglia. This study considers a simplified domain with heterogeneity represented by fractures of negligible storage that allow water flow bypassing the NAPL ganglia entrapped within the homogeneous permeable blocks. Here the nonuniformity of NAPL dissolution originates from the fracture network, which may transfer water with low contaminant concentration into regions of high NAPL saturation. This effect leads to nonuniform variations in the interphase mass transfer and in the increased concentration of solutes along the groundwater streamlines with the reduced entrapped NAPL saturation. Decreasing NAPL saturation affects the permeability and porosity of the porous block, and the rate of interphase mass transfer from the entrapped NAPL to the mobile water phase. The study evaluates the parameters controlling the FPA pump-and-treat cleanup. In regions contaminated by entrapped NAPL, our simulations concern whether the reduced permeability of the porous blocks enables a significant increase of the fracture flow at the expense of the permeable block flow. In the real world, different types of heterogeneous porous media are possible. Therefore, we expect that there may be cases in which quantifying NAPL entrapment and dissolution is more complicated than the case analyzed by this study. 3. Basic formulas Table 1 presents a list of the different assumptions used in developing the simplified conceptual model of this study. The following paragraphs discuss these assumptions. A 2-D model is not well suited for representing scenarios typical of point source contamination, such as the bypassing of a NAPL-contaminated zone (e.g., Frind et al., 1999); these scenarios require 3-D modeling. However, in cases of large-scale contamination, as is common at airports, airbases, and transportation terminals and is observed in the CPA area, a 2-D model may describe phenomena taking place in a substantial part of the domain. The assumption of a fracture network with equidistant fractures of constant aperture simplifies the model, but ignores some phenomena. In the real world, NAPL may be entrapped in the tighter parts of a fracture network with variable

aperture, leading to flow channeled within the fracture network that may bypass the entrapped NAPL. However, when one considers a large-scale NAPL entrapment such as the case of the CPA, the assumptions represented in Table 1 do not cause major errors in characterizing FPA cleanup. Because of the symmetry of the 2-D domain shown in Fig. 2(a), simulations and analysis of this study refer to the subdomain shown in Fig. 2(b). This subdomain incorporates fracture segments, permeable blocks, and matrix sections. The set of successive fracture segments composes the “fracture network” of the subdomain. The term “permeable block” refers to the permeable medium occupying the space between adjacent fracture segments. The term “matrix section” refers to the rectangular space incorporating two halves of adjacent blocks and a single embedded fracture segment. In each permeable block the cross-section at its centerline represents the interface between two adjacent matrix sections. It is possible to assume two different cleanup procedures. The first consists of injecting and pumping water under a constant head difference between the entrance and the exit cross-sections of the contaminated domain. This choice leads to increasing aquifer flow during the cleanup, because NAPL dissolution increases overall aquifer permeability. The second procedure consists of constant discharges of the injection and pumping wells, that is, a total discharge with constant flow through the contaminated aquifer. From a practical viewpoint the difference between these two procedural choices is minor. Hence, calculations targeting the efficiency of the cleanup procedure should be consistent with either one of these choices. In the remainder of this study, only the second choice is considered, namely, cleanup subject to constant aquifer discharge. Through each cross-section of the 2-D subdomain shown in Fig. 2(b), which is perpendicular to the flow direction, a two-dimensional volumetric discharge Qt flows in the longitudinal direction of the x⁎-coordinate. This discharge incorporates the fracture segment flow rate, Qf, and the single permeable block flow rate, Qb. Note that Qf represents flow in the direction of the fracture segment, whereas Qb represents flow mainly in the direction of the hydraulic gradient. Transverse flow of the blocks in the y⁎-direction is very small because of the horizontal top and bottom boundary conditions of the domain shown in Fig. 2(a) and the hydraulic gradient that is in the x⁎-direction. Some of our numerical simulations allowing longitudinal and transverse flows have justified the preceding phrasing. In parts of the domain that are free of entrapped NAPL, Qf and Qb have constant values, Qf0 and Qb0, respectively. In the

H. Rubin et al. / Journal of Contaminant Hydrology 96 (2008) 128–149

contaminated part of the domain, the relation among the three types of flow rates is: Z B sin h qb dy⁎ ð1Þ Qt ¼ Qf þ Qb ; Qb ¼ 0

where y⁎ is the vertical coordinate and qb is the specific discharge of the permeable block flow, which is mainly in the x⁎ direction, that is, qb ≈ qbx. Fig. 2(c) shows an elementary fracture volume applied to describe solute transport through the fracture network. In this figure, C⁎ is the solute concentration in the fracture flow, Cb⁎ is the solute concentration of the permeable block flow entering the fracture segment, x⁎ and y⁎ are the horizontal and vertical coordinates, respectively, and x′ is a local coordinate extended along the fracture centerline. It is assumed that solute transport through the fracture occurs solely by advection. The schematics of fluxes shown in Fig. 2(c) neglects effects of diffusion, longitudinal dispersion, and mass storage changes within the fracture. Previous studies, reviewed by Birkhölzer et al. (1993a), have shown that these assumptions are reasonable for transport in FPAs, where the ratio of flow through the permeable block to flow through the fractures is NM0 = 0.1–25. This ratio refers specifically to sandstone aquifers. In a particular part of the Coastal Plain of Israel, Gvirtzman et al. (1988) found NM0 = 3. This study originates from evidence that the aquifer is contaminated by kerosene. Kerosene is a mixture of many organic compounds, most of them olefins. Fig. 2(c) refers to a single-species NAPL, but the basic approach of this study may consider multiple-species systems. According to Fig. 2(c), mixing between the permeable block flow and the fracture flow affects solute transport in the fracture network. The following paragraphs represent the development of the mathematical model, which is based on mass balances in the fracture and permeable block flows. By considering steady-state flow, employing the foregoing assumptions and neglecting all second-order terms, the mass balance obtained on the water phase and the solute in the elementary fracture volume are, respectively: ∂Qf dxVþ ðqbex  qben Þdy⁎ ¼ 0 ∂xV

ð2Þ

  ∂ðQf C ⁎ Þ dxVþ qbex C ⁎  qben Cb⁎ dy⁎ ¼ 0; ∂xV

ð3Þ

where qben and qbex are the block-specific discharges entering and leaving the elementary fracture volume, respectively. Because of the negligible storage of the

135

fracture elementary volume, there are no time-dependent terms in Eqs. (2) and (3). As shown in Appendix A, Eqs. (2) and (3) yield a single differential Eq. (A3) that concerns the solute in the fracture network. The following dimensionless coordinates, time, and concentrations are applied to simplify Eq. (A3): x¼

x⁎ y⁎ t ⁎ Vb0 C⁎ C⁎ ; y¼ ; t¼ ; C ¼ ⁎ ; Cb ¼ b⁎ ; B cos h B sin h B cos h Cs Cs

ð4Þ where Cs⁎ is the equilibrium concentration of solute in the water phase, Vb0 is the interstitial permeable block flow velocity in parts of the domain that are free of entrapped NAPL, and t⁎ is time. Later in this paper, the study considers time variations of quantities while referring to the permeable block flow. In Eq. (4), variables with an asterisk superscript represent the relevant dimensional quantity. By introducing the variables of Eq. (4) into Eq. (A3) we get: ∂C þ NM C ¼ NM Cb ; ∂x

ð5Þ

where NM is the local mobility number, defined as NM ¼

qB B sin h ; Qf

ð6Þ

Qf is the local fracture flow rate, and qB (which is a particular value of qb) is the block-specific discharge entering the fracture segment. In Eq. (6), both Qf and qB are values of the fracture flow rate and the block-specific discharge at the fracture segment point at which the value of NM is calculated. In parts of the domain contaminated by entrapped NAPL, Qf and qB may be subject to changes due to NAPL dissolution. Therefore, in such parts of the domain the mobility number may vary along the fracture network. In parts of the domain that are free of entrapped NAPL the mobility number is: NM0 ¼

qb0 B sin h : Qf 0

ð7Þ

Eq. (7) indicates that in NAPL-free parts of the domain, the mobility number, NM0, represents the ratio of the permeable block flow rate, Qb0, to the fracture flow rate, Qf0. However, according to Eq. (5), in contaminated parts of the domain the mobility number, NM, does not represent this ratio of discharges. As shown in Appendix B, the parameter QR, which is equal to the product of NM0 and the average relative block permeability, represents the ratio between these flow rates in all parts of the domain. We may consider NM0 as

136

H. Rubin et al. / Journal of Contaminant Hydrology 96 (2008) 128–149

a basic measure of the preferential flow in the domain. To provide sufficient perspective on the effect of the mobility number on FPA cleanup, that is, the effect of groundwater flow bypassing the fracture network, we carry out simulations at NM0 = 0.5, 2, 5. We find that at NM0-values higher than 5, our model suffers from numerical dispersion, resulting from the reduced hydrodynamic dispersion due to reduced mixing between the fracture and permeable block flows. It should be noted that our model refers to ranges of mobility numbers, in which mixing between the fracture and permeable block flows is the only source of solute dispersion. Further, simulations at NM0-values lower than 0.5 with domain parts subject to high saturation of NAPL entrapment the ratio of Qb to Qf becomes very small. Then diffusion of organic solute from the porous blocks into the fracture flow should not be neglected. In other words, the entrapped NAPL reduces the permeable block flow rate, and thereby it reduces the effective mobility number. Therefore, simulations with our model for NM0 ≈ 0.1 should be avoided. Appendix B incorporates details concerning the effect of block permeability changes on the distribution of the fracture flow rate and the specific block discharge in the domain. Various expressions developed in Appendix B are helpful for quantifying the effects of NAPL dissolution and solute transport in the domain. Various studies provide approximate expressions for the effect of the entrapped NAPL on the relative water permeability, krw, of the porous medium. In Appendix C, we consider how such expressions can be useful for this study. As two fluid phases are present in the domain, the following association obtains: Sw þ Sn ¼ 1;

ð8Þ

where Sn is the NAPL saturation. The study of Birkhölzer et al. (1993a,b) has shown that mixing between the fracture flow and the permeable block flow generally causes overall dispersion typical of FPAs. According to the simulations of Birkhölzer et al. (1993b), by neglecting dispersion in the block, the flow in FPAs approaches the state of plug flow only at mobility numbers larger than 25. If the mobility number is smaller than 25 then the value of dispersivity caused by mixing between the fracture and the block flows is significantly larger than values of dispersivity in the block flow itself. Simulations of this study suggest mobility values up to 5. Therefore, it is possible to calculate NAPL dissolution and solute transport (conservation of solute in the water phase) in the block flow with negligible contaminant dispersion in the blocks.

Under such conditions, material balance for an elementary volume of the permeable block yields ∂ð/b Sw Þ ∂qbx ∂qby þ ⁎ þ ⁎ ¼0 ∂t ⁎ ∂x ∂y

ð9Þ

for the water phase. Material balance for the organic mass dissolved in the water phase yields:         ∂ /b Sw Cb⁎ ∂ qbx Cb⁎ ∂ qby Cb⁎ þ þ ¼ kf Cs⁎  Cb⁎ : ∂t ⁎ ∂x⁎ ∂y⁎ ð10Þ In Eqs. (9) and (10), qbx and qby are components of the specific discharge in the x⁎ and y ⁎ directions, respectively; ϕb is the porosity of the permeable blocks; and kf is the lumped coefficient of the mass transfer from the NAPL phase to the water phase. It is common to apply the Sherwood number to study the effect of the domain and the entrapped NAPL properties on the value of kf. The definition of the Sherwood number is: Sh ¼

kf d 2 ; D

ð11Þ

where Sh is the Sherwood number, d is the characteristic grain diameter, and D is the free liquid diffusivity. Various experimental studies (e.g., Miller et al., 1990, Imhoff et al., 1994, Powers et al., 1992, 1994, Saba and Illangasekare, 2000, Nambi and Powers, 2003) have shown that the Sherwood number is dependent on the entrapped NAPL saturation and the Reynolds number. The Sherwood number also depends on some other quantities of the domain, which are not affected by the entrapped NAPL dissolution. According to those studies (see e.g., Mayer and Miller, 1996), it is possible to represent the value of Sh by: Sh ¼ b0 Reb1 Snb2 :

ð12Þ

Here, βi (i = 0–2) are coefficients, and Re is the Reynolds number based on the interstitial block flow velocity. The coefficient β0 is a dimensionless quantity that depends on the porous medium grain size and size of the contaminated region. Unger et al. (1998) noticed that at Re = 0, Eq. (12) yields Sh= 0, and hence they added a coefficient to avoid such an error in their simulations. However, we may assume that Eq. (12) is applicable to simulating aquifer cleanup by pump-and-treat, which requires significant flow velocities. Various studies (e.g., Guiguer and Frind, 1994, Unger et al., 1998) carried out simulations with the kf-value affected only by the entrapped NAPL saturation. After introducing their correction to Eq. (12) that avoids the value of Sh vanishing at Re = 0, Unger et al. (1998)

H. Rubin et al. / Journal of Contaminant Hydrology 96 (2008) 128–149

found that the mean dissolution time was particularly sensitive to the value of β2, and the uncertainty of the simulations was sensitive to the Reynolds number. Substituting Eq. (9) into Eq. (10), noting that qby is very small (namely, qb ≈ qbx), and applying the dimensionless quantities of Eq. (4), the following equation of solute transport in the blocks obtains: ∂Cb ∂Cb þ vr ¼ Kf ð1  Cb Þ; ð13Þ ∂t ∂x where vr is the relative interstitial block flow velocity, that is, the ratio of the local interstitial block flow velocity to the velocity in parts of the domain that are free of entrapped NAPL, and Kf is the dimensionless interphase mass transfer coefficient. The expressions representing these parameters are, respectively: vr ¼

Vb kf B cos h ; Kf ¼ : qb0 Sw Vb0

ð14Þ

Eq. (13) shows that the block flow collects solutes from the entrapped NAPL. When this flow enters the fracture segment, it mixes with the fracture flow, as shown in Fig. 2(c). According to Eq. (14), Kf is a measure of the rate of interphase mass transfer, which is proportional to the lumped mass transfer coefficient, and is affected by the specific discharge of the block flow, the water saturation of the block, and the length scale of the fractures, namely, the quantity Bcosθ. According to Eq. (12), changes in Sn, Sw, and vr during the NAPL dissolution lead to changes in the interphase mass transfer coefficient, as:  Kf ¼ Kf 0

vr vr0

b1 

Sn Sn0

 b2 

 Sw0 ; Sw

ð15Þ

where the subscript 0 refers to initial values of the relevant quantity in parts of the domain contaminated by entrapped NAPL. Fig. C2 in Appendix C implies that during the NAPL dissolution, the effect of changing vr / vr0 on the value of Kf is substantially smaller than the effect of the decreasing value of Sn. Therefore, we ignore the effect of Re on the value of Kf. Further, we assume β2 = 2/3. This value of β2 fits the assumption of spherical ganglia of the entrapped NAPL, whose number is kept constant during the NAPL dissolution. Based on all data in the above-cited references about measured values of Sherwood numbers, we consider the relevant range of initial values for the interphase mass transfer coefficient to be 0.1 ≤ Kf0 ≤ 10. Simulations of Rubin et al. (2001) showed that for Kf0 values above 10 the system approach equilibrium partitioning. For Kf0 values smaller than 0.1, and a domain of limited length, the aquifer cleanup starts

137

at a very small flux average outflow contaminant concentration, and it extends for a very long time. Similar to the derivation of Eq. (10), by referring to an elementary volume of the porous block, the conservation of the total amount of organic material (NAPL and solute) in the pores of the block yields: ∂Sn ∂Cb⁎ þ q ¼ 0; ð16Þ b ∂t ⁎ ∂x⁎ where ρn is the NAPL density. Introducing the dimensionless quantities of Eq. (4) into Eq. (16) yields:

/b qn

∂Sn ∂Cb þ qr Cnv ¼ 0: ∂t ∂x

ð17Þ

Here qr is the ratio of qb to qb0, and Cnv is the volumetric equilibrium concentration of the organic compound, which is given by: Cnv ¼

Cs⁎ : qn

ð18Þ

This dimensionless quantity is needed to connect the decreasing entrapped NAPL saturation with the increasing organic solute concentration in the permeable block flow. Due to the symmetry of the domain shown in Fig. 2(a), simulations of the aquifer cleanup pertain to the subdomain of Fig. 2(b). Therefore, such simulations incorporate calculations of average cross-sectional NAPL saturation in the permeable blocks of the subdomain of Fig. 1(b), and average cross-sectional flux concentrations of the solute in that subdomain. The average crosssectional NAPL saturation, Snav, and flux average concentration, C⁎av, of the solute (in a cross-section perpendicular to the flow direction) are, respectively: Snav ¼

1 B sin h

Z

⁎ ⁎ ¼ Qf C þ Cav

¼

B sin h

Sn dy⁎

ð19Þ

0

R B sin h 0

Qt

qb Cb⁎ dy⁎

Qf CCs⁎ þ Cs⁎ B sin h Qt

R1 0

qb Cb dy

:

ð20Þ

By substituting Eq. (1) into Eq. (20) and performing some simple mathematical operations, we get:  R 1 h C þ 0 qb BQsin Cb dy f Cav ¼ ; ð21Þ 1 þ QR where Cav is the dimensionless flux average concentration of the dissolved contaminant in the flowing water.

138

H. Rubin et al. / Journal of Contaminant Hydrology 96 (2008) 128–149

Introducing Eqs. (B8) and (B11) into Eq. (21) yields:  Z 1 1 Cav ¼ Cb krw dy : C þ NM0 R1 1 þ NM0 0 krw dy 0

Cb = 0 at t = 0. Under such conditions, integrating Eq. (24) between starting and completing the FPA cleanup yields:

ð22Þ

ð25Þ

This expression and the expression for the average NAPL saturation of the block cross-section are employed to analyze the characteristics of the FPA cleanup with quantities typical of a continuum. According to Eq. (19) the average NAPL saturation of the block is: Z 1 Snav ¼ Sn dy: ð23Þ

   Z T 1 Sn0 nc ¼ 1 þ Cav dt ; Cnv NM0 0 x¼nc

where T is the dimensionless time required to complete the aquifer cleanup. Eq. (25) represents the identity between the total amount of initially entrapped NAPL and the total amount of solute passing through the downstream crosssection of the subdomain between starting and completing the cleanup.

0

4. The numerical model It is important to check the contaminant's overall mass balance during the simulation of the pump-andtreat cleanup. Such a balance implies that the rate of removing the organic material (NAPL and dissolved organic compound) from the contaminated subdomain is identical with the solute flux at the downstream end of the subdomain. By applying dimensionless quantities, the following expression is obtained for this identity: ∂ ∂t

Z 0

nc

Z

1

ðSn þ Sw Cnv Cb Þdydx

0

  1 þ 1þ Cnv ðCav Þx¼nc ¼ 0: NM0

ð24Þ

As the initial quantity of organic compounds dissolved in the water phase is very small, we may let

A finite difference grid of small squares is applied to solve the basic equations with dimensionless quantities that were developed in the preceding sections. The model covers physical domains with no limit over the fracture orientation. However, as shown by Eq. (4), normalizing the x⁎ and y⁎ coordinates is achieved by applying the length scales Bcosθ and Bsinθ, respectively. This leads to a fracture orientation of 45° only in the dimensionless domain. Therefore, the numerical calculation with the applied grid easily incorporates the finite difference intervals of the permeable blocks and those of the fracture network. The grid intervals shown in Fig. 3 are Δx = Δy. Our preliminary numerical experiments have led to using Δx = Δy = 0.1 and Δt = 0.05. Referring to the symbols given in this figure, the simulated domain

Fig. 3. Schematic description of the numerical grid with enumeration of grid points, permeable blocks, fracture segments, and sections.

H. Rubin et al. / Journal of Contaminant Hydrology 96 (2008) 128–149

incorporates nc sections contaminated by entrapped NAPL. The region of interest is that of the subdomain of the nc contaminated sections in which NAPL is subject to dissolution and the solute is subject to advection and dispersion. Fig. 3 shows that values of i (i = 1,…, imax) represent locations of grid nodal points in the vertical direction. Values of j represent locations of grid nodal points in the longitudinal direction. The part of the domain contaminated by entrapped NAPL incorporates blocks and fractures up to the longitudinal nodal points given by jc = [nc(imax − 1) + 1]. The symbol ifr is applied to enumerate sections and fracture segments. The symbol ibl is applied to enumerate blocks. Values of k (k = 1,…, imax) represent numbers of nodal points of each fracture segment. The fracture network is divided into consecutive segments and the finite difference approximation is applied to Eq. (5) to calculate the dimensionless fracture flow concentration, C, at the k + 1 nodal point of the ifr fracture segment at the m time step:   Dx ðmÞ Dx ðmÞ ðmÞ ðmÞ Ckþ1 1 þ NMkþ1 ¼ Ck 1  NMk 2 2 h i Dx ðmÞ ðmÞ ðmÞ ðmÞ NMkþ1 Cbkþ1 þ NMk Cbk ; þ ð26Þ 2 where the superscript m represents the time step number, Cbk and Cbk + 1 are values of the dimensionless block flow concentration, Cb entering the nodal points k and k + 1 of the ifr fracture segment, respectively, and NMk and NMk + 1 are values of the mobility number at nodal points k and k + 1 of the ifr fracture segment, respectively. Variations of vr and Kf in a single time step are assumed to be very minor. Thus, the finite difference scheme, both forward in time and backward in space, can be applied to get the following implicit–explicit finite difference approximation of Eq. (13):   1 ðmÞ 1 ðmÞ ðmþ1Þ ðmÞ ðmÞ Dt þ K Dt ¼ Cbi;j 1  Kf i;j Dt Cbi;j 1 þ vri;j Dx 2 f i;j 2  ðmþ1Þ ðmÞ Dt ðmÞ þCbi;j1 vri;j þ Kf i;j Dt: ð27Þ Dx The finite difference approximation of Eq. (17) is applied to calculate the entrapped NAPL saturation at the midpoint of the longitudinal interval as: h i Þ ðmÞ ðmÞ ðmÞ ðmÞ ðmÞ Dt : ¼ S  C q C  q C Snðmþ1 nv ri;jþ1 bi; jþ1 ni; jþ1=2 ri; j bi; j i; jþ1=2 Dx ð28Þ

139

The flux average cross-sectional concentration (at the cross-section perpendicular to the flow direction and represented by the longitudinal nodal points ja) is calculated by applying the finite difference approximation of Eq. (22), trapezoidal numerical integration, and proper reference to the crossed fracture segment. Consider an integer number α satisfying the following expression for the longitudinal nodal point ja, where Cav is calculated: aðimax  1Þ þ 1bja bða þ 1Þðimax  1Þ þ 1:

ð29Þ

Under the condition of Eq. (29), the cross-section chosen for calculating Cav crosses the fracture segment at the segment nodal point kα, which is ka ¼ ja  aðimax  1Þ  1;

ð30Þ

and the finite difference approximation of Eq. (22) yields Cav ¼ ½1 þ DxNM0 INT ðkrw Þ1 ½Ckr þ DxNM0 INT ðkrw Cb Þ;

ð31Þ where INT ðKrw Þ ¼

 1 krw1 þ krwkb þ krwka þ krwimax 2 kX imax 1 a 1 X þ krwi þ krwi i¼ka þ1

i¼2

INT ðkwr Cb Þ ¼

1h

ð32Þ

ðkrw Cb Þ1 þðkrw Cb Þka þðkrw Cb Þka þðkrw Cb Þimax 2 Ximax 1 Xka 1 ðkrw Cb Þi i¼k þ1 ðkrw Cb Þi : þ i¼2

i

a

ð33Þ The subscripts in Eqs. (31)–(33) refer to vertical positions of the fracture and block nodal points whose common longitudinal position is ja, and values of all quantities are calculated at the same time step. The subscripts kα− and kα+ of Cb refer to values upstream and downstream of the fracture segment, respectively. Eventually, values of Ckα and Cbkα+ are identical. If the integer α satisfies the expression aðimax  1Þ þ 1 ¼ ja ;

ð34Þ

then the finite difference approximation of Eq. (31) collapses to: (

Cav ¼

"

#)1 imax 1 X 1 ðkrw1 þ krwimax Þ þ krwi 2 i¼2 ( " #) 1 h i imax X 1 Ckb þ DxNM0 ðkrw Cb Þ1 þðkrw Cb Þimax þ ðkrw Cb Þi : 2 i¼2 1 þ DxNM0

ð35Þ

140

H. Rubin et al. / Journal of Contaminant Hydrology 96 (2008) 128–149

Again, subscripts refer to vertical positions of the nodal points whose common longitudinal position is ja, and all quantities' values are calculated at the same time step. By using cell NAPL saturations, Eq. (23) leads to the following finite difference expression for Snav: Snav ¼ Dy

imax 1 X

Sniþ1=2 :

ð36Þ

i¼1

The finite difference approximation of Eq. (B14) yields the following expression for the average crosssectional relative permeability of the FPA subdomain: " # imax X 1 krt ¼ þ Dy 0:5ðkrw1 þ krwimax Þ þ krwi : ð37Þ NM0 i¼2 Expressions for Cav, Snav, and krt may represent the contaminated domain by quantities borrowed from the

continuum approach. Variations of these quantities during the NAPL dissolution may provide information concerning the feasibility of the FPA cleanup. According to the information reported by Rubin et al. (2001), the present study may refer to sandstone aquifer with properties similar to those of the CPA part contaminated by entrapped NAPL. In this FPA the block hydraulic conductivity is about Kb0 = 10− 5 m/s and the porosity is about ϕb = 0.2. The assumed longitudinal projection of the fracture spacing is about Bcosθ = 0.5 m. In parts of the domain that are free of entrapped NAPL, the assumed average hydraulic gradient imposed by the pump-and-treat procedure is about 0.01. Then the permeable block-specific discharge is approximately qb0 ≅ 10− 7 m/s and the flow velocity is Vb0 ≅ 5 × 10− 7 m/s. Under such conditions, one dimensionless time unit of the variable t represents 106 s, or about 11.6 days. It should be noted that the present study

Fig. 4. Examples representing temporal changes of exit Cav-values: (a) NM0 = 2, Sn0 = 0.3; (b) NM0 = 0.5, Sn0 = 0.3.

H. Rubin et al. / Journal of Contaminant Hydrology 96 (2008) 128–149

only incorporates analysis and evaluations based on carrying out numerical calculations. Such results require validation and verification based on experimental and possibly field results that will follow this study. 5. Numerical simulation results Simulations were carried out with subdomains like that shown by Fig. 3 with up to 12 consecutive sections. In our study, we have followed changes in time and along the fractures of Qf-values and solute concentration, C. In the permeable blocks, we have followed changes in distributions of entrapped NAPL saturation, contaminant concentration, Cb, hydraulic conductivity, specific discharge and interstitial flow velocity. Our simulations have provided information about NAPL disappearance from the subdomain. It is possible to apply averaged bulk quantities with terminology of the continuum approach to represent the domain properties, for instance, bulk average NAPL saturation, bulk average value of Kf, etc. However, the continuum approach usually cannot be applied to characterize NAPL dissolution in the subdomain of Fig. 3, which is subject to high initial entrapped NAPL concentration. Fig. 4 shows the temporal changes of the exit Cav-value at the downstream end of the cross-section of the subdomain, where x = nc. The figure describes the effect of the reduced block permeability on the increasing exit value of Cav during the initial stage of the aquifer cleanup. In this paragraph we explain the origin of this phenomenon. During the initial stage of the remediation, a very significant part of the groundwater discharge flows through the fractures, which are free of entrapped NAPL, namely at the domain exit crosssection Qf is large. Later, due to NAPL dissolution the value of Qb increases, and the value of Qf decreases (as our simulations refer to constant value of the total groundwater discharge). However, it should be noted that values of Qf and Qb are subject to changes in time and also along the domain. Also the specific discharge of the permeable block flow qb is subject to changes in time and place. The increasing value of Qb leads to increasing solute mass discharge at the exit crosssection of the domain. This process is ended with a maximum value of the flux average concentration at the exit cross-section of the domain, because at a certain time the decreased entrapped NAPL saturation should lead to reducing flux average solute concentration at the exit cross-section of the domain. The phenomenon of increasing and then decreasing exit flux average concentration is expected in a short domain, which

141

incorporates a small number of permeable blocks. In long domains, which include a large number of permeable blocks, the mixing between Qf and Qb at each fracture segment may lead to complete loading of the fracture flow with solute concentration of equilibrium. However, as discussed in the next paragraph some other parameters affect the phenomenon of increasing and then decreasing exit flux average solute concentration. The phenomenon of increasing and later decreasing Cav-values occurs only when the permeability of the block is reduced by the entrapped NAPL, namely in cases of using “Corey” and “Parker” models (see Appendix C). Further, this phenomenon is most significant at average values of Kf0, namely, at Kf0 ≈ 1. At high Kf0-values of Kf0 ≈ 10, the high rate of interphase mass transfer overcomes the effect of reduced block permeability, and the flowing groundwater is loaded with solutes to its full capacity (equilibrium). At very low Kf0-values of Kf0 ≈ 0.1, the exit value of Cav is low, and the reduced block permeability has a minor effect on this value. However, the effect of the reduced permeability on temporal changes of Cav-values at the exit cross-section is more significant at low values of the mobility number, namely in Fig. 4b than in Fig. 4a. The total area underneath each curve of Fig. 4a, or b, is represented by the LHS of the expression: Z T  Sn0 nc  Cav dt ¼ : ð38Þ 0 x¼nc 1 þ N1M0 Cnv Eq. (38) represents the identity of the total solute outflow with the amount of material initially entrapped as NAPL within the porous blocks. The simulations of this study preserved the identity of Eq. (38) with an error of up to 8%. Further, at each time step the accuracy of the numerical simulations was controlled by checking the vanishing value of the LHS of this equation: d dt

Z

nc 0

Z

1

ðSn þ Sw Cnv Cb Þdydx

0

  1 þ 1þ Cnv ðCav Þx¼nc ¼ 0: NM0

ð39Þ

Eq. (39) represents the overall organic material mass conservation in the subdomain subject to the numerical simulation. Fig. 4 provides some information about the complete cleanup time. However, it should be noted that due to the scale of depicting Fig. 4 some cleanup curves show

142

H. Rubin et al. / Journal of Contaminant Hydrology 96 (2008) 128–149

an abrupt end of the cleanup. However, always the cleanup completion takes more time than may visually be implied by Fig. 4. Two major quantities provide information relevant to the cleanup efficiency: • The dimensionless treatment time for complete cleanup, T, and • The dimensionless total volume of water used to dissolve and flush the originally entrapped NAPL (the treated water volume), W. The latter quantity, W, represents the ratio of the quantities of water used during the entire aquifer cleanup to the water incorporated within a single permeable block free of entrapped NAPL. The relation between the dimensionless remediation time and volume is:   1 W ¼ 1þ T: NM0

ð40Þ

According to Eq. (40), for identical values of the dimensionless cleanup time, T, the dimensionless water volume, W, is twice as large for NM0 = 0.5 than for NM0 = 2. The numerical results represented by Fig. 4 show that often the value of T is smaller for a case of NM0 = 0.5 than for a case with NM0 = 2. However, the water volume needed for complete remediation is always larger for the case of NM0 = 0.5 while all other quantities are identical in both cases. The decreased block permeability due to the entrapped NAPL, increases the time and volume of water needed to achieve complete cleanup. Therefore, in cases of using the “No K model”, values of T and W increase almost linearly with increases in the initial entrapped NAPL saturation, whereas the “Parker model” observes a nonlinear increase in these quantities with an increased Sn0-value. This effect originates from the significantly nonlinear dependence of permeability on NAPL saturation as implied by the “Parker model”. Note that in Fig. 4, we consider fixed values of Qb0. Under such conditions, a high mobility number suggests the total flow rate of the aquifer is small and vice versa. Fig. 4(b), which depicts the results corresponding to a low mobility number, shows shorter cleanup times than those of Fig. 4(a). However, the decrease in cleanup time is not inversely proportional to the increased aquifer flow rate. Therefore, the total volume of water needed to complete the treatment increases with a decrease in the mobility number. Quantitative presentation of this effect is obtained by substituting the results shown by Fig. 4 into Eq. (40). Simulations involving

Fig. 5. Comparing the numerical simulation results of this study according to the “Corey model” with NM = 0.5 with the experimental results of Geller and Hunt (1993) and Nambi and Powers (2000).

small Sn0-values, namely, Sn0 ≤ 0.1, showed that differences in the time and the water volume needed to achieve complete cleanup predicted by the “No K model” were quite similar to those predicted by the “Parker model”. However, in cases of high Sn0-values, namely, Sn0 ≈ 0.3, the differences between the time and the water volume needed to complete cleanup predicted by the “No K model” and the “Parker model” were significant. Again, this originates from the nonlinear relationship between permeability and NAPL saturation in the “Parker model”. We have found the total cleanup time and the water volume consumed during the FPA cleanup were almost inversely proportional to the logarithm of Kf0. Also, in cases of low mobility numbers the increasing value of Kf0 was more effective than in cases of high mobility numbers. In practice, the pump-and-treat procedure should probably apply a compromise between two opposing needs: 1) the minimum time period necessary to clean up the FPA, and 2) the minimum amount of treated water needed to achieve cleanup. The initially increasing exit value of Cav according to the “Parker model” implies that the rate of FPA cleanup is affected by low block permeability. Under such conditions low flushing water discharge can be considered to decrease the amount of water bypassing the blocks via the fracture network, implying a lower value of qb and in turn a higher value of Kf. Thus, in the final stages of the FPA cleanup, reduced water discharges increase the value of Kf. According to the present study, in FPAs the initially increasing and then decreasing solute concentration at

H. Rubin et al. / Journal of Contaminant Hydrology 96 (2008) 128–149

the exit cross-section is predicted. This phenomenon originates from the reduced permeability caused by the entrapped NAPL and the fractures that allow groundwater flow to bypass the entrapped NAPL ganglia. Some experimental studies have shown the phenomenon of initially increasing and then decreasing outflow contaminant concentrations. Although the experimental setups used in these studies are different from the model of FPA of this study, the phenomenon observed there also originates from water flow bypassing the domain part contaminated by entrapped NAPL. Geller and Hunt (1993) carried out experiments with columns of a uniform porous medium composed of glass beads, in which entrapped NAPL was initially entrapped in a small part of the domain. Due to the decreased permeability originating from the entrapped NAPL, initially there was very small water discharge flowing through the domain part contaminated by entrapped NAPL. At that time, dissolution of the entrapped NAPL increased the hydraulic permeability and thereby increased the water flow through the contaminated part, leading to increasing flux average solute concentration at the exit cross-section of the experimental setup. Later, the reduced NAPL saturation led to decreasing exit flux average solute concentration. Nambi and Powers (2000, 2003) carried out experiments with fine sand columns, in which lenses of coarse sand contaminated with entrapped NAPL were inserted. Their experiments started with water flowing mainly through the fine sand, which was free of entrapped NAPL, and finally the water discharge was flowing mainly through the coarse sand lens. Therefore, they also observed the phenomenon of increasing and later decreasing flux average contaminant concentration at the exit cross-section of the experimental setup. Fig. 5 presents the experimental results of Geller and Hunt (1993) and Nambi and Powers (2000, 2003). To facilitate comparison, we convert our numerical results with the “Corey model”, NM0 = 0.5, and Sn0 = 0.3 into dimensionless outflow contaminant concentration versus pore volumes that fit the size of the experimental setups of Geller and Hunt (1993) and Nambi and Powers (2000, 2003) and present these transformed results in Fig. 5. The figure shows some qualitative similarity between the numerical results of this study and the experimental results. However, there are significant differences between the porous media used in the experimental setups of the two sources of data, and none of them is quantitatively relevant to this study. Therefore, complete validation of the modeling approach of this study requires additional sets of experimental results.

143

6. Summary and conclusions This study develops a simplified model of entrapped NAPL dissolution and solute transport in FPAs. A FPA represents a particular type of heterogeneous aquifer in which the fracture network is characterized by high permeability and negligible storage. FPAs are fractured aquifers defined by the mobility number range of 0.1 ≤ NM0 ≤ 25. By applying dimensionless coordinates and variables, this study identifies two dimensionless parameters, Kf0 and NM0, which measure the effect of aquifer heterogeneity on flow, NAPL dissolution, and solute transport in the domain. Besides these parameters, two types of constitutive relationships affect aquifer cleanup: 1) the relationship between NAPL saturation and permeability of the blocks, and 2) the dependence of the lump mass transfer coefficient, kf, on NAPL saturation and the block flow velocity. Within the permeable blocks, the initial residual NAPL saturation is taken to be uniform. The basic formulation of the model considers a steady and constant total flow rate across the domain. The preceding discussion suggests there is an initial uniform division of flow over the domain between permeable blocks and the fracture network, as NAPL entrapment is initially uniform. However, it is expected that the flow through the permeable blocks will gradually and nonuniformly dissolve the entrapped NAPL. Dissolving the entrapped NAPL increases the permeability and the porosity of the porous blocks. The nonuniform decrease of the NAPL saturation changes the division of the total aquifer flow between the porous blocks and the fracture segments, that is, the characteristics of the preferential flow in the domain. The decreasing NAPL saturation decreases the local interphase mass transfer coefficient, because of changes in the interfacial contact areas and changes in the interstitial block flow velocity. The dimensionless interphase mass transfer coefficient is a measure of interphase mass transfer intensity. The mobility number is a measure of mechanisms available for water flow bypassing the entrapped NAPL ganglia. Both parameters control the preferential flow, NAPL dissolution, and solute transport in the domain. Expressions were developed to follow overall mass balance in the domain subject to pump-and-treat cleanup. These expressions characterize the water phase flow, dissolution of NAPL in the permeable blocks, transport of solutes in the permeable block, and the fracture flows, as well as mixing between the permeable block and fracture flows. This paper emphasizes using formulas with dimensionless quantities to evaluate NAPL dissolution and transport

144

H. Rubin et al. / Journal of Contaminant Hydrology 96 (2008) 128–149

phenomena in the domain. The paper emphasizes the numerical simulation results representing the effect of NAPL dissolution on changes of Cav-values at the outflow cross-section of the domain. The results show that at high initial saturations of the entrapped NAPL, the contaminant concentration increases during initial stages of the FPA cleanup, but later it decreases. This phenomenon originates from significant groundwater flow bypassing the NAPL entrapped in the permeable blocks via the fracture network. Acknowledgements The first author (H. Rubin) appreciates the support of the College of Engineering, and in particular, the School of Civil and Environmental Engineering at Cornell University, which allowed him to visit the School in the 2003–04 academic year as the Mary Sheppard B. Upson Visiting Professor. During this year he carried out parts of this research. This study has been partially supported by The Grand Water Research Institute (GWRI) and the Fund for Research Cooperation between Technion and Haifa University. It was also supported by RWTH Aachen and grant no 1573/8-1 from the German Research Foundation (DFG). The authors thank the Editor-in-Chief Emil Frind and the anonymous referees for their constructive comments and suggestions. Appendix A. Solute transport in the fractures The relations between the domain coordinate intervals and the interval dx′ of the local coordinate of the elementary fracture volume, shown in Fig. 1(c), are: dx⁎ ¼ dx Vcos h;

dy⁎ ¼ dx Vsin h:

ðA1Þ

We multiply Eq. (2) by C⁎ and introduce it into Eq. (3); then we apply Eq. (A1) to obtain:   ∂C ⁎ dxVþ qben C ⁎  Cb⁎ dx Vsin h ¼ 0: ∂xV Eq. (A2) can be rewritten as:

Qf

  ∂C ⁎ Qf ⁎ þ qben tan h C ⁎  Cb⁎ ¼ 0: ∂x

ðA2Þ

NAPL dissolution is probably not uniform in the contaminated region. Therefore, the specific discharge, qb, as well as interstitial block flow velocity, Vb, is not necessarily uniformly distributed in the contaminated blocks of the domain. Eq. (5) is obtained by introducing the dimensionless variables of Eq. (4) into Eq. (A3). Appendix B. The distribution of discharges The hydraulic gradient, J, applied in the x⁎-direction varies along the subdomain of Fig. 2(b), but we assume that the value of J does not change in the vertical direction. This assumption originates from the boundary conditions (flat top and bottom) of the domain that do not allow significant flow in the y⁎-direction. This assumption extends the approximation of Dupuit for the flow through a FPA contaminated by entrapped NAPL. Therefore, the following expression is obtained:   R B sin h Qt ¼ J Mf cos h þ 0 Kb dy⁎   R1 ¼ J Mf cos h þ B sin hKb0 0 krw dy ;

where Kb0 is the hydraulic conductivity of the blocks, which are free of entrapped NAPL, krw is the local relative permeability of the block, and Mf is the mobility of the fracture segment, that is, the ratio of fracture flow rate to the hydraulic gradient imposed along the fracture segment. As the fractures are assumed to be almost free of entrapped NAPL, the value of Mf is constant; often, the value of Mf is calculated based upon the observation that the fracture flow is similar to viscous flow between parallel plates. During NAPL dissolution, the value of Kb increases, but Qt is constant in all cross-sections. Therefore, according to Eq. (B1), the value of J decreases during aquifer cleanup. Eq. (B1) is applied to express J and J0 (the hydraulic gradient and its values in regions free of entrapped NAPL, respectively). Using these expressions yields the following expressions for Qf, qB, and qb0: Qt Mf cos h

Qf ¼

Mf cos h þ B sin hKb0

ðA3Þ

The assumed initial conditions are of uniform distribution of entrapped NAPL saturation in the NAPLcontaminated matrix domain. This saturation reduces the permeability of the permeable blocks in the contaminated part of the domain. Later, the water saturation and the block permeability increase because of the entrapped NAPL dissolution, at which point

ðB1Þ

qB ¼

Qt Kb0 krwB Mf cos h þ B sin hKb0

R1 0

R1 0

ðB2Þ

krw dy

krw dy

where

krwB ¼

KB Kb0

ðB3Þ qb0 ¼

Qt Kb0 : Mf cos h þ B sin hKb0

ðB4Þ

H. Rubin et al. / Journal of Contaminant Hydrology 96 (2008) 128–149

145

In Eq. (B3) the term Kb0 krwB represents the value of KB. This quantity is the hydraulic conductivity at the point located to the left of the vertical cross-section crossing the fracture segment; krwB is the relative hydraulic conductivity at that point. Eqs. (6) and (B2)–(B4) lead to another expression for the mobility number: NM ¼

BKB tan h BKb0 krwB tan h ¼ : Mf Mf

ðB5Þ

In parts of the domain that are free of entrapped NAPL, the mobility number is: NM0 ¼

BKb0 tan h : Mf

ðB6Þ

Eq. (B6) is applied to represent Mf. Using this expression with Eqs. (B3) and (B4) yields: ! qB 1 þ NM0 ¼ krwB ; R1 qb0 1 þ NM0 0 krw dy KB NM ¼ NM0 krwB ; krwB ¼ : Kb0

ðB7Þ

Eqs. (B2)–(B7) also imply that QR, the ratio of the permeable block flow rate to the fracture segment flow rate, is: Z QR ¼ NM0

1

ðB8Þ

krw dy: 0

In regions free of entrapped NAPL, QR = NM0. As streamlines are almost horizontal, Eq. (B7) leads to the following expression: !

qb 1 þ NM0 qr ¼ ckrw ; R1 qb0 1 þ NM0 0 krw dy

ðB9Þ

where qr is the local relative specific discharge, and krw is the local relative water permeability. Flow in the blocks is not entirely horizontal, particularly if preferential dissolution pathways develop in the permeable matrix. Therefore, Eq. (B9) incorporates the approximate identity instead of an equal sign. Eq. (B9) leads to the following expression for the local relative interstitial block flow velocity: qr 1 þ NM0 vr ¼ c R1 Sw 1 þ NM0 0 krw dy

!

 krw : Sw

ðB10Þ

Fig. C1. Curves of relative water permeability versus water saturation according to various models.

It is possible to define the average cross-sectional specific discharge of the subdomain of Fig. 1(b) as: qt ¼

Qt : B sin h

ðB11Þ

The connection between the hydraulic gradient and the average cross-sectional specific discharge yields the following expression for the average cross-sectional hydraulic conductivity of the subdomain in regions free of entrapped NAPL: Kt0 ¼

Mf þ Kb0 : B tan h

ðB12Þ

Using the same approach for regions contaminated with entrapped NAPL, the average cross-sectional hydraulic conductivity is given by: Mf 1 Kt ¼ þ B tan h B sin h

Z

B sin h

Kb dy⁎ :

ðB13Þ

0

Dividing Eqs. (B12) and (B13) by Kb0 and applying the dimensionless quantities of Eq. (4) yields: krt0 ¼

1 1 þ 1; krt ¼ þ NM0 NM0

Z

1

krw dy;

ðB14Þ

0

where krt and krt0 are the relative cross-sectional permeabilities for regions with and without entrapped NAPL, respectively. During aquifer cleanup, the following relationships hold: NM YNM0 and krt Ykrt0 ; where Sn ; Sw Y0; 1:

ðB15Þ

146

H. Rubin et al. / Journal of Contaminant Hydrology 96 (2008) 128–149

krw ¼ 1:0;

Appendix C. Calculating effects of permeability changes The constitutive relationships between block permeability and NAPL saturation should be used to apply the expressions developed in Appendix B. Various studies of the 1950's and 1960's (Irmay, 1954; Corey, 1954, 1957; Wyllie, 1962, Hausenberg and Zaslavsky, 1963; Topp and Miller, 1965) showed by theoretical and experimental means that  n krw ¼ Sw⁎ ; ðC1Þ where n is a power coefficient equal to 3 or 4, depending on the type of soil, and Sw⁎ is the effective or reduced saturation of water, which is given by Sw⁎ ¼

Sw  Sm ; 1  Sm

ðC2Þ

where Sm is the irreducible water saturation. Brooks and Corey (1964) have suggested the following expression for n: n¼

2 þ 3k ; k

at

Sw b0:8;

ðC5Þ

and (here we have introduced some minor changes, to avoid inconsistency) krw ¼

ðSw  0:16Þ2 ; 0:41

krn ¼

ð0:8  Sw Þ2 0:41

ðC6Þ at

0:16 V Sw V 0:18;

where krw and krn are relative water and NAPL permeability, respectively, and Sw is the water saturation. The value of Sm implied by Eqs. (C4) and (C5) is 0.16, which is an average value for various types of porous media. Eq. (C6) suggests that the maximum residual saturation of entrapped immobile NAPL is 0.2. However, laboratory and field measurements have shown the entrapped NAPL saturation may be significantly larger than this value. On the other hand, it is possible to apply Eq. (C6) to some soils, provided that the residual NAPL saturation is low.

ðC3Þ

where λ is the pore size distribution index. Therefore, for λ → ∞ the power coefficient n approaches the value 3, and it increases with decreasing values of λ. Whitaker (1977) argued that n should be in the range of 2 ≤ n ≤ 5. Demond and Roberts (1993) showed that n = 3 fitted data of sand contaminated by xylene. Geller and Hunt (1993) also applied n = 3. However, Frind et al. (1999) applied n = 4 for simulating dense NAPL (DNAPL) dissolution at the Canadian Forces Base Borden. Similarly, Nambi and Powers (2003) applied n = 4 for coarse sand. This study uses n = 3 for Eq. (C1) as the “Corey model” to represent a moderate effect of water saturation on block permeability. Theoretically, the value of Sm is close to 0.0 for a perfect water-wet system. However, experimental studies with different technologies and different cores (e.g., Tovar, 1997) and various types of granular porous media (e.g., Bryant and Johnson, 2003, Fleury and Deflandre, 2003, Chen, 2004) obtained values of 0.05 ≤ Sm ≤ 0.4. Low values of Sm were found for granular material with uniform grain size distribution and particle size larger than 5 mm (Chen, 2004). Wilson et al. (1990) claimed the irreducible water saturation generally lies in the range of 0.05–0.15. Faust (1985) performed simulations of a onedimensional waterflood, in which he considered: krw ¼ 0:0;

krn ¼ 0:0

krn ¼ 1:0

at

Sw b0:16

ðC4Þ

Fig. C2. Curves of relative block flow velocity in NAPL-contaminated FPAs versus water saturation: (a) according to the “Corey model", (b) according to the bParker model".

H. Rubin et al. / Journal of Contaminant Hydrology 96 (2008) 128–149

Some simulations of this study follow the approach of Faust (1985) and consider no effect of entrapped NAPL saturation on the matrix permeability and porosity. This type of simulation is called “No K model”. Parker et al. (1987) applied the analysis of VanGenuchten (1980) and Mualem (1978, 1984) to air– water systems. Arguing that the permeability of the wetting fluid phase (water) does not depend on the type of the nonwetting fluid phase (air or NAPL); they developed an expression for relative water permeability that depends on capillary pressure: n h inp o2 krw ¼ Sw⁎1=2 1  1  Sw⁎1=np ;

ðC7Þ

where np is an empirical constant whose value according to Parker et al. (1987) is 0.4565 for sand and 0.4624 for clay. However, the value of np may be very different for different soils. In Eq. (C7) we use np = 0.46 in “Parker Model” simulations concerning significant dependence of block permeability on water saturation in the range of interest, namely, 0.7 ≤ Sw ≤ 1.0. Mayer and Miller (1996) and Unger et al. (1998) applied a similar approach using the “Parker model”. Fig. C1 shows the different curves of permeability versus water saturation, which are obtained by using the models of Faust (1985), Corey (1954, 1957), and Parker et al. (1987). We have applied Eq. (B10) to depict the effect of changes in values of Sw and NM0 on the normalized block flow velocity, vr, as shown in Fig. C2. In this figure, we consider a uniform distribution of water saturation in the permeable block cross-section. The figure shows ranges of vr-values that are expected according to both the “Corey model” with n = 3 and the “Parker model”. For 0.7 ≤ Sw ≤ 1.0, the range of vr-values implied by the “Corey model” is 0.5 ≤ vr ≤ 1.0, and the range implied by the “Parker model” is 0.1 ≤ vr ≤ 1.0. References Anderson, M.R., Johnson, R.L., Pankow, J.F., 1992a. Dissolution of dense chlorinated solvents into groundwater. 1. Dissolution from a well-defined source. Ground Water 30 (2), 250–256. Anderson, M.R., Johnson, R.L., Pankow, J.F., 1992b. Dissolution of dense chlorinated solvents into groundwater. 3. Modeling contaminant plumes from fingers and pools of solvent. Environmental Science Technology 26 (5), 901–907. Andres, G., Georgotas, N., 1978. Hydrogeologische Studien zum Grundwasserhaushalt und zur Stoffbilanz im Maineinzugsgebiet. Schriftenreihe Bayr. Landesamt für Wasserwirtschaft, München, Germany. Averbach, Y., 1988. Analysis of Pumping Tests of Barka 1 Well. TahalWater Planning of Israel. Tel-Aviv, Israel.

147

Barenblatt, G.I., Zheltov, Y.P., Kochina, L.N., 1960. Basic concepts in the theory of seepage of homogeneous liquids in fissured rocks. PMM-Sov. Applied Mathematics and Mechanics 852–864. Berkowitz, B., Bear, J., Braester, C., 1988. Continuum models for contaminant transport in fractured porous formations. Water Resources Research 24 (8), 1225–1236. Birkhölzer, J., Rubin, H., Daniels, H., Rouvè, G., 1993a. Contaminant transport in fractured permeable formation. 1. Parametric evaluation and analytical solution. Journal of Hydrology 144, 1–33. Birkhölzer, J., Rubin, H., Daniels, H., Rouvè, G., 1993b. Contaminant transport in fractured permeable formation. 1. Numerical solution. Journal of Hydrology 144, 35–58. Brooks, R.H., Corey, A.T., 1964. Hydraulic properties of porous media. Hydrology Papers. Colorado State University, Fort Collins. Brusseau, M.L., 1992. Rate-limited mass transport of organic solutes in porous media that contain immobile immiscible organic liquid. Water Resources Research 28 (1), 33–45. Brusseau, M.L., Zhang, Z., Nelson, N.T., Cain, B., Tick, G.R., Oostrom, M., 2002a. Dissolution of nonuniformly distributed immiscible liquid: intermediate-scale experiments and mathematical modeling. Environmental Science and Technology 36 (5), 1033–1041. Brusseau, M.L., Nelson, N.T., Oostrom, M., Zhang, Z., Johnson, G.R., Wietsma, T.W., 2002b. Influence of heterogeneity and sampling method on aqueous concentrations associated with NAPL dissolution. Environmental Science and Technology 34 (17), 3657–3664. Bryant, S., Johnson, A., 2003. Wetting phase connectivity and irreducible saturation in simple granular media. Journal of Colloid and Interface Science 263 (2), 572–579. Chatiz, I., Morrow, N.R., Lim, H.T., 1983. Magnitude and detailed structure of residual oil saturation. Society of Petroleum Engineers Journal 311–326 April. Chen, R.G., 2004. Measurement and numerical simulation of moisture transport by capillary gravity and diffusion in porous potash beds, M.Sc. thesis, Dept. of Mechanical Engineering, Univ. of Saskatchewan, Saskatoon, SK, Canada. Corey, A.T., 1954. The interrelation between gas and oil relative permeabilities. Producer's Monthly 1 (19), 38–41. Corey, A.T., 1957. Measurement of water and air permeability in unsaturated soils. Proceedings — Soil Science Society of America 21, 7–10. Demond, A.H., Roberts, P.V., 1993. Estimation of two-phase relative permeability relationships for organic liquid contaminants. Water Resources Research 29 (4), 1081–1090. Faust, C.R., 1985. Transport of immiscible fluids within and below the unsaturated zone — a numerical model. Water Resources Research 21 (4), 587–596. Fleury, M., Deflandre, F., 2003. Quantitative evaluation of porous media wettability using NMR relaxometry. Magnetic resonance imaging 21 (3–4), 385–387. Frind, E.O., Molson, J.W., Schirmer, M., Guiguer, N., 1999. Dissolution and mass transfer of multiple organics under field conditions: the Borden emplaced source. Water Resources Research 35 (3), 683–694. Geller, J.T., Hunt, J.R., 1993. Mass transfer from nonaqueous phase organic liquids in water-saturated porous medium. Water Resources Research 29 (4), 833–846. Guiguer, N., Frind, E.O., 1994. Dissolution and mass transfer processes for residual organics in the saturated groundwater zone. In: Dracos, Th., Stauffer, F. (Eds.), Transport and Reactive Processes in Aquifers. Belkama, Rotterdam, pp. 475–480.

148

H. Rubin et al. / Journal of Contaminant Hydrology 96 (2008) 128–149

Gvirtzman, H., Magaritz, M., Kanfi, Y., Carmi, I., 1988. Matrix and fissure water movement through unsaturated calcareous sandstone. Transport In Porous Media 3, 343–356. Hatfield, K., Stauffer, T.B., 1993. Transport in porous media containing residual hydrocarbon. 1. Model. ASCE Journal of Environmental Engineering 119 (3), 540–558. Hausenberg, L., Zaslavsky, D., 1963. The effect of size of water stable aggregates on field capacity (in Hebrew). Department of Civil Engineering Technion — Israel Institute of Technology, Haifa. Publ. No. 35. Imhoff, P.T., Jaffe, P.R., Pinder, G.F., 1994. An experimental study of complete dissolution of a nonaqueous phase liquid in saturated porous media. Water Resources Research 30 (2), 307–320. Irmay, S., 1954. On the hydraulic conductivity of unsaturated soils. Transactions — American Geophysical Union 35, 463–468. Kanfi, Y., 1986. Groundwater Contamination by Oil in the Coastal Plain Aquifer of Israel. Ministry of Agriculture, Water Commission, Tel-Aviv. Mackay, D.M., Cherry, J.A., 1989. Groundwater contamination: pump-and-treat remediation. Environmental Science and Technology 23 (16), 630–636. Mackay, D.M., Roberts, P.V., Cherry, J.A., 1985. Transport of organic contaminants in groundwater. Environmental Science and Technology 19 (5), 384–392. Mayer, A.S., Miller, C.T., 1996. The influence of mass transfer characteristics and porous media heterogeneity on nonaqueous phase dissolution. Water Resources Research 32 (6), 1551–1567. Mercer, J.W., Cohen, R.M., 1990. A review of immiscible fluids in the subsurface: properties, models, characterization and remediation. Journal of Contaminant Hydrology 6, 107–163. Miller, C.T., Poirier-McNeill, M.M., Mayer, A.S., 1990. Dissolution of trapped nonaqueous phase liquids: mass transfer characteristics. Water Resources Research 26, 2783–2796. Miller, C.T., Christakos, G., Imhoff, P.T., McBridge, J.F., Pedit, J.A., 1998. Multiphase flow and transport modeling in heterogeneous porous media: challenges and approaches. Advances in Water Resources 21 (2), 77–120. Mualem, Y., 1978. Hydraulic conductivity of unsaturated porous media generalized macroscopic approach. Water resources research 14 (2), 325–334. Mualem, Y., 1984. Prediction of the soil boundary wetting curve. Soil Science 137 (6), 379–390. Nambi, I.M., Powers, S.E., 2000. NAPL dissolution in heterogeneous systems: an experimental study in simplified heterogeneous system. Journal of Contaminant Hydrology 44 (2), 161–184. Nambi, I.M., Powers, S.E., 2003. Mass transfer correlations for nonaqueous phase liquid dissolution from regions with high initial saturations. Water Resources Research 39 (2), 1030–1040. Parker, J.C., Lenhard, R.J., Kuppusamy, T., 1987. A parametric model for constitutive properties governing multiphase flow in porous media. Water Resources Research 23 (4), 618–624. Pettijohn, F.J., Potter, P.E., Siever, R., 1987. Sand and Sandstone. Springer, New York. Powers, S.E., Loureiro, C.O., Abriola, L.M., Weber Jr., W.J., 1991. Theoretical study of the significance of nonequilibrium dissolution of nonaqueous phase liquids in subsurface systems. Water Resources Research 27 (4), 463–477. Powers, S.E., Abriola, L.M., Weber Jr., W.J., 1992. An experimental investigation of nonaqueous phase liquid dissolution in saturated subsurface systems: steady state mass transfer rates. Water Resources Research 28 (10), 2691–2706.

Powers, S.E., Abriola, L.M., Weber Jr., W.J., 1994. An experimental investigation of nonaqueous phase liquid dissolution in saturated subsurface systems: transient mass transfer rates. Water Resources Research 30 (2), 321–332. Rubin, H., Braester, C., 2000. Field Measurements, Laboratory Tests, and Theoretical Analysis of Oil Contamination in the Coastal Plain Aquifer of Israel (in Hebrew). Department of Civil Engineering, Technion — Israel Institute of Technology, Haifa, Israel. Rubin, H., Buddemeier, R.W., 1996. Transverse dispersion of contaminants in fractured permeable formations. Journal of Hydrology 176, 133–151. Rubin, H., Rathfelder, K., Abriola, L.M., 1997. Modeling quasi-steady NAPL dissolution in permeable fractured media. ASCE Journal of Environmental Engineering 123 (3), 205–216. Rubin, H., Rathfelder, K., Abriola, L.M., Spiller, M., Demny, G., Köngeter, J., 2001. The effect of fractures on the reclamation of NAPL contaminated aquifers. International Austrian Technion Society Symposium, Preserving the Quality of our Water Resources, April 23–25, 2001, Vienna, Austria. Rubin, H., Rathfelder, K., Spiller, M., Köngeter, J., 2004. Continuum quantifying remediation of fractured permeable formations. ASCE Journal of Environmental Engineering 130 (11), 1345–1356. Saba, T., Illangasekare, T., 2000. Effect of groundwater flow dimensionality on mass transfer from entrapped nonaqueous phase liquid contaminants. Water Resources Research 36 (4), 971–979. Schwille, F., 1984. Migration of organic fluids immiscible with water in the unsaturated zone. In: Yaron, B., et al. (Ed.), Pollutants in Porous Media. Springer, Berlin, pp. 27–48. Schwille, F., 1988. Dense Chlorinated Solvents in Porous and Fractured Media. Lewis Publishers, Chelsea, MI. translated by J.F. Pankow. Topp, G.C., Miller, E.E., 1965. Hysteretic conductivity calculation versus measurement. 1964 Meeting of American Society of Agronomy, Agronomy Abstracts. Tovar, R.A., 1997. Measurements of relative permeability for steamwater flow in porous media, M.Sc. thesis, Department of Petroleum Engineering, Stanford University, Stanford, CA. Unger, A.J.A., Forsyth, P.A., Sudicky, E.A., 1998. Influence of alternative dissolution models and subsurface heterogeneity on DNAPL disappearance times. Journal of Contaminant Hydrology 30, 217242. U.S. Environmental Protection Agency (EPA), 1996. Pump-and-Treat Ground-Water Remediation: A Guide for Decision Makers and Practitioners. EPA 625/R-95/005, Washington, D.C. U.S. Environmental Protection Agency (EPA), 2001. A Citizen's Guide to Pump-and-Treat. EPA 542-F-01-025, Washington D.C. U.S. National Research Council (NCR), 1994. Alternatives for Ground Water Cleanup. National Academy Press, Washington D.C. Van Genuchten, M.T., 1980. A closed form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Science Society Am. Journal 44, 892–898. Whitaker, S., 1977. Simultaneous heat, mass, and momentum transfer in porous media: a theory of drying. Advances in Heat Transf 13, 119–203. Wilson, J.L., Conrad, S.H., Mason, W.R., Peplinski, W., Hagan, E., 1990. Laboratory Investigations of Residual Organic Liquids from Spills Leaks and Disposal of Hazardous Wastes in Groundwater. U.S. EPA, Washington, D.C.. EPA/600/6-90/004.

H. Rubin et al. / Journal of Contaminant Hydrology 96 (2008) 128–149 Wyllie, M.R.J., 1962. Relative permeability in oil production handbook. Reservoir Engineering, vol. II. McGraw-Hill, New York. 1962. Zaidel, J., Russo, D., 1993. Analytical models of steady state organic species transport in the vadose zone with kinetically controlled volatilization and dissolution. Water Resources Research 29 (10), 3343–3356.

149

Zhang, Z., Brusseau, M.L., 1999. Nonideal transport of reactive solutes in heterogeneous porous media simulating regional-scale behavior of a trichloroethene plume during pump-and-treat remediation. Water Resources Research 35 (10), 2921–2935.