Interphase mass transfer between fluids in subsurface formations: A review

Berken Agaoglu; Nadim K. Copty; Traugott Scheytt; Reinhard Hinkelmann

This paper presents a review of the state-of-the-art on interphase mass transfer between immiscible fluids in porous media with focus on the factors that have significant influence on this process. In total close to 300 papers were reviewed focusing to a large extent on the literature relating to NAPL contamination of the subsurface. The large body of work available on this topic was organized according to the length scale of the conducted studies, namely the pore, meso and field scales. The interrelation of interphase mass transfer at these different scales is highlighted. To gain further insight into interphase mass transfer, published studies were discussed and evaluated in terms of the governing flow configurations defined in terms of the wettability and mobility of the different phases. Such organization of the existing literature enables the identification of the interfacial domains that would have significant impact on interphase mass transfer. Available modeling approaches at the various length scales are discussed with regard to current knowledge on the physics of this process. Future research directions are also suggested.

Advances in Water Resources 79 (2015) 162–194 Contents lists available at ScienceDirect Advances in Water Resources journal homepage: www.elsevier.com/locate/advwatres Review Interphase mass transfer between fluids in subsurface formations: A review Berken Agaoglu a,c,⇑, Nadim K. Copty b, Traugott Scheytt a, Reinhard Hinkelmann c a Institute of Applied Geosciences, Berlin Institute of Technology, Sekr. BH 3-2, Ernst-Reuter-Platz 1, 10587 Berlin, Germany Institute of Environmental Sciences, Bogazici University, 34342 Bebek, Istanbul, Turkey c Water Resources Management and Modeling of Hydrosystems, Department of Civil Engineering, Berlin Institute of Technology, TIB 1-B14, Gustav-Meyer-Allee 25, 13355 Berlin, Germany b a r t i c l e i n f o Article history: Received 24 March 2014 Received in revised form 15 February 2015 Accepted 16 February 2015 Available online 23 February 2015 Keywords: Review Interphase mass transfer in porous media Multiphase flow NAPL Dissolution Mathematical models a b s t r a c t This paper presents a review of the state-of-the-art on interphase mass transfer between immiscible fluids in porous media with focus on the factors that have significant influence on this process. In total close to 300 papers were reviewed focusing to a large extent on the literature relating to NAPL contamination of the subsurface. The large body of work available on this topic was organized according to the length scale of the conducted studies, namely the pore, meso and field scales. The interrelation of interphase mass transfer at these different scales is highlighted. To gain further insight into interphase mass transfer, published studies were discussed and evaluated in terms of the governing flow configurations defined in terms of the wettability and mobility of the different phases. Such organization of the existing literature enables the identification of the interfacial domains that would have significant impact on interphase mass transfer. Available modeling approaches at the various length scales are discussed with regard to current knowledge on the physics of this process. Future research directions are also suggested. Ó 2015 Elsevier Ltd. All rights reserved. Contents 1. 2. 3. Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Fundamentals of interphase mass transfer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.1. Ionic strength and temperature . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2. Multi-component NAPLs and aging . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3. Cosolvents and surfactants . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Mechanics of interphase mass transfer. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1. Pore scale . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2. Meso scale . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3. Field scale . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4. Factors influencing interphase mass transfer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4.1. Interfacial area . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4.2. Hysteresis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4.3. Entrapment . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4.4. Wettability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4.5. Mean grain diameter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4.6. Pore size distribution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4.7. NAPL mobilization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4.8. Desorption. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4.9. Density and viscosity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ⇑ Corresponding author. E-mail address: Berken@hotmail.com (B. Agaoglu). http://dx.doi.org/10.1016/j.advwatres.2015.02.009 0309-1708/Ó 2015 Elsevier Ltd. All rights reserved. 164 165 167 167 167 168 168 170 171 172 172 174 175 176 176 177 177 177 178 B. Agaoglu et al. / Advances in Water Resources 79 (2015) 162–194 4. 5. Mathematical modeling of interphase mass transfer. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.1. Analytical models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2. Pore scale numerical models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2.1. Pore network models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2.2. Other pore scale models. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3. REV-based numerical models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3.1. Single phase models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3.2. Semi-multiphase models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3.3. Multiphase models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Summary and conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 163 178 178 180 180 182 183 183 183 184 187 188 Nomenclature Roman letters A absolute interfacial area (L2) Aia the air–water interfacial area normalized by porous medium volume (1/L) Aij cross section area of the pore throat between pore i and j (L2) ANW interfacial area between wetting and non-wetting phase per unit pore volume (1/L) Ana;i interfacial area in pore i (L2) AsN interfacial area between solid and non-wetting phase per unit pore volume (AsN ¼ 0 if Sw ¼ 1) (1/L) Aw cross sectional area in water filled pore throat (L2) a specific interfacial area (1/L) awn effective specific interfacial area (meniscus interfacial area per unit volume of porous media (1/L) C solute concentration (M/L3) Ca solute concentration at the center of each subdivision (M/L3) C a;i solute concentration in pore i (M/L3) S solute equilibrium concentration (M/L3) Ca Cb concentration in bulk phase (M/L3) C eq solute equilibrium concentration (M/L3) Ci solute concentration in pore i (M/L3) in solute concentration entering pore i (M/L3) Ci solute concentration leaving pore i (M/L3) C out i in solute concentration entering pore j (M/L3) Cj C NW constant of integration C out effluent concentration (M/L3) Cs solute equilibrium concentration (M/L3) C ðsÞ solute equilibrium concentration (M/L3) solute concentration entering pore throat (M/L3) C in t D molecular diffusion coefficient (L2/T) Dm molecular diffusion coefficient (L2/T) Dy dispersion coefficient in y direction (L2/T) ds diameter of spherical blobs (L) dm grain diameter (L) E interphase mass transfer rate divided by the interfacial area Ed energy dissipation factor (Eq. (10)) normalized average effluent conc. (Eq. (25)) Ed Ewn interfacial area production rate (1/L T) ei separation distance between nodes (L) ewn strength of change in specific interfacial area (1/L) F ðjÞ force acting on position x of component j (M L/T2) F ij mass flux between pore i and j (M/L2 T) Fðxj xi ; yj yi ; zj zi Þ erfc type transport function (open format can be found in [220]) Gij pore conductance between pores i and j (T L4/M) Gjk interaction strength parameter between component j and k (Eq. (29)) GTP Ii Kf Kl s K kf 1;2 l l lt L M M ðtÞ Mi M0 m N N Ni Ni n Pij Pcrit c Pe Q ij Qw ij Q aij q q Re Ri Rt SA Sai Si Sm Sw s Sh Sh0 ti U Ux V a;i Vi Vs Vs ganglia to pool ratio interphase mass transfer term ratio of NAPL density to the contaminant conc. in the flushing solution in stream tube i (Eq. (19)) mass transfer rate coefficient (i.e., mass transfer coefficient lumped with interfacial area) (1/T) average hydraulic conductivity of domain (L/T) mass transfer coefficient (L/T) length of hypothetical film thickness (L) length of the pore throat (Eq. (23)) (L) length of the pore throat (L) length in x direction (L) solute mass transferred between phases (M) NAPL mass at time t (M) number of interfaces connected to pore i initial NAPL mass (M) parameter determined with respect to uniformity coefficient of the medium (Eq. (7)) number of subdivisions (Eqs. (13) and (14)) number of adjacent pores to pore i (Eq. (22)) number of pores connected to pore i molar flux of component i (Eq. (30)) (mol/L2 T) parameter determined with respect to uniformity coefficient of the medium (Eq. (7)) pressure difference between pore i and j (M/L T2) critical disjoining pressure (M/L T2) Peclet number flow rate between pores i and j (L3/T) flow rate of water between pore i and j (L3/T) flow rate of phase a between pore i and j (L3/T) superficial aqueous phase velocity (L/T) average Darcy velocity (L/T) Reynolds number retardation factor of stream tube i pore throat radius (L) specific surface area of porous medium (L2/M) local saturation of phase a in pore i trajectory average NAPL saturation of stream tube i monolayer saturation wetting phase saturation normalized surface area of the porous media Sherwood number modified Sherwood number nonreactive travel time in stream tube I (T) velocity (L/T) velocity in x direction (L/T) volume of pore containing water (L3) volume of pore I (L3) volume of NAPL in shape of spherical blobs per unit volume of medium (Eq. (12)) (L3/L3) volume of site (pore i) (Eq. (23)) (L3) 164 x x; y; z B. Agaoglu et al. / Advances in Water Resources 79 (2015) 162–194 geometric parameter related to pore wall curvature that is positive for convex shape (Eq. (28)) lengths in x, y and z-directions (L) Greek letters a parameter determined with respect to uniformity coefficient of the medium (Eq. (7)) b mass depletion exponent b1 velocity dependence parameter b2 mass depletion exponent jeff mass transfer rate coefficient (1/T) j0 parameter found by fitting the analytical solution to numerical results Pmax maximum disjoining pressure (M/LT2) 1. Introduction Multiphase flow in porous media is encountered in numerous natural and industrial systems. Water seepage through the vadose zone during precipitation events is just one example of the many natural systems involving the flow of multiple phases. Applications involving multiphase flow include fuel cells, oil recovery, non-aqueous phase liquid (NAPL) contamination/remediation, CO2 injection into deep saline aquifers, and numerous industrial applications involving drying of porous media such as coatings, food, paper, textile, wood, ceramics, building materials, granular materials, electronic devices and pharmaceuticals [1–5]. A fundamental process that is common to all these multiphase systems is interphase mass transfer, which is the transfer of components across the interface separating the different phases. Given the centrality of interphase mass transfer to many applications, significant efforts across a number of disciplines have been directed to characterize this process. Earlier efforts were mostly in relation to petroleum exploration and in the field of chemical engineering, dealing mostly with engineered packed bed systems [6,7]. Starting in the late twentieth century when NAPL contamination of the subsurface was recognized as one of the most challenging environmental problem, significant research was conducted to elucidate the factors influencing interphase mass transfer such as NAPL dissolution. Much of the advances of these research efforts consisted of laboratory experiments that led to the definition of empirical expressions of interphase mass transfer (e.g., [8,9]). These expressions were mostly specific to the conditions present during the experiments and a unified interphase mass transfer theory was beyond reach. Major limitations have been the inability to define the spatial distribution of the multiphase system, including the distribution of the interfacial areas and the transient flow in the surrounding regions. In the past two decades significant developments in experimental technologies such as non-invasive imaging techniques and pore scale modeling have meant that the multiphase system within the porous medium can be more accurately characterized. In this work we present a review of the state-of-the-art on interphase mass transfer in porous media with focus on the factors that have significant influence on that process. The molecular basis of interphase mass transfer is presented first, followed by the mechanics of interphase mass transfer when multiple fluids exist in porous media. In the latter section the existing literature is grouped under different ‘‘fluid-porous medium’’ configurations. As interphase mass transfer is directly influenced by the contact area and the flow properties (velocity, path) of the solvent phase, such distinction between different two-phase flow conditions in porous media enables a comparative evaluation and further insight into this process. Table 1 presents 5 different flow configurations that can exist in two phase Interfacial tension (M/T2) interfacial tension between non-wetting and wetting phases (M/T2) si reactive travel time in stream tube i (T) / porosity of the porous medium (Eq. (10)) /sNW contact angle (°) UNW area under capillary pressure–saturation curve ½UNW ðSw ÞD;I;D area under capillary pressure–saturation curve for primary drainage, secondary and all subsequent imbibition, and main and all subsequent drainage conditions. W ratio of the effective to total specific interfacial area WðjÞ effective mass given as a function of local density (Eq. (29)) r rNW Table 1 Classification of two phase flow in porous media for the evaluation of interphase mass transfer. Acronym Dissolving phase Solvent phase Examples INtoMW Immobile, Nonwetting Immobile, Wetting Mobile, Nonwetting Mobile, Wetting Mobile, Wetting Non-wetting NAPL dissolution Air sparging IWtoMN MNtoIW MWtoIN CHEM Mobile, Nonwetting Immobile, Wetting Immobile, Nonwetting Mobility and wettability of phases is altered through the use of miscibility enhancing chemical agents CO2 injection Gas exsolution Enhanced NAPL remediation flow in porous media. These are; INtoMW: Mass transfer from the Immobile Non-wetting phase to the Mobile Wetting phase; IWtoMN: Mass transfer from the Immobile Wetting phase to the Mobile Non-wetting phase; MNtoIW: Mass transfer from the Mobile Non-wetting phase to the Immobile Wetting phase; MWtoIN: Mass transfer from the Mobile Wetting phase to the Immobile Non-wetting phase; CHEM: Mass transfer when miscibility-enhancing CHEMical agents are present. Here after these flow configurations will be referred to with their acronym defined above and in Table 1. ‘‘INtoMW’’ is encountered in applications such as NAPL dissolution in subsurface systems contaminated by non-wetting NAPLs whereas ‘‘IWtoMN’’ can be observed during remediation activities such as air sparging, NAPL dissolution of wetting NAPLs and applications involving drying of porous medium. ‘‘MNtoIW’’ occurs during in CO2 injection into saline aquifers and ‘‘MWtoIN’’ is observed in, for example, gas exsolution applications. ‘‘CHEM’’ occurs in enhanced oil recovery and enhanced NAPL remediation activities. The advantages of the categorization adopted in this review paper are further discussed in Section 3.1. It should be mentioned that both ‘‘MNtoIW’’ and ‘‘MWtoIN’’ have been investigated only recently and as a results there are only a handful of reported studies on these configurations. Therefore they are mentioned only briefly in this review. Although the research on interphase mass transfer relating to NAPL contamination of the subsurface constitutes the foundation of this review, the discussion of the physics of this process enables the transition of the findings to other applications with similar conditions. Within the environmental field, interphase mass transfer has been commonly investigated at three different scales; (i) pore scale, (ii) meso-scale, (iii) field scale [2,10,11]. Pore scale is the scale where the flow properties at discrete pore bodies/throats can be observed. The characteristic length of the meso-scale mostly B. Agaoglu et al. / Advances in Water Resources 79 (2015) 162–194 ranges from 1 to 5 cm while the field scale ranges from a decimeter to hundreds of meters where the influence of the large scale heterogeneity becomes significant. Although interphase mass transfer at these different scales is intrinsically related, the experimental and mathematical tools used for their investigation can be significantly different. Here we review the research on all three scales and attempt to shed light on the relationship between the different scales. For each of the scales, the discussion is also framed in terms of the different flow configurations defined above. The rest of this paper is organized as follows: Section 2 presents basic definitions of mass transfer and its fundamentals including its basis in molecular chemistry. Section 3 discusses the mechanics of interphase mass transfer in porous media with an emphasis on scale issues. Mathematical models developed for the simulation of interphase mass transfer in porous media are presented in Section 4. All mathematical terms appearing in these models are defined in the nomenclature. Section 5 presents a critical assessment of the reviewed works and identifies areas for future research. 2. Fundamentals of interphase mass transfer The partitioning of one type of molecule into a different phase environment is determined by the combined effect of the intermolecular forces namely: the London dispersive forces, dipole– dipole forces and dipole induced forces [12–16]. These forces are much weaker than the intramolecular forces such as covalent bonds [12]. The partitioning can be defined thermodynamically in terms of the chemical potential of the compound. If the chemical potentials of a compound in two adjacent systems (in this case; phases) are the same, this means that equilibrium is present and no net interphase mass transfer of that compound would occur [16]. A commonly adopted approach to determine the interphase mass transfer rate is to define it as the product of a coefficient of interphase mass transfer that is unique for the system conditions and the concentration difference (or chemical potential) between the interface, where it is assumed to be at solubility, and the bulk solvent phase where sufficient mixing can be assumed to occur allowing for the definition of a representative concentration over some specific volume. Since the interphase mass transfer is function of the interfacial area across which the mass transfer is occurring, the interfacial area is also incorporated in this approach resulting in [8]: @M ¼ kf 1 AðC s C Þ @t ð1Þ Various theories have been developed to relate interphase mass transfer to system parameters based on the type of the boundary conditions existing at the interface. These theories, which include the film theory, boundary layer theory, penetration theory and surface renewal theory, were mostly developed for an air–water interface which in nature often exhibits turbulence conditions especially at the air side [17–20]. As a result the influence of turbulence on diffusivity is considered to be significant in the development of these theories. The film theory, which is based on bottleneck boundary conditions (i.e., abrupt decrease of diffusivity adjacent to an interface), assumes that a stagnant thin layer develops at the interface between the two phases where molecular diffusion is the governing mass transport process. Consequently, the mathematical description of interphase mass transfer proposed by the film theory expresses the interphase mass transfer coefficient in terms of the diffusion coefficient and a hypothetical film thickness value l, @M D ¼ ðC s C Þ A@t l ð2Þ The boundary layer theory on the other hand assumes a gradual decrease in diffusivity (i.e., the sum of molecular and turbulent 165 diffusivities) with approach to the interface, while the interphase mass transfer coefficient is defined through the Schmidt number [20]. Detailed description of these theories can be found in [13,21,22]. It has been well shown in the literature that one of the most influential parameters on interphase mass transfer is the velocity. The dependence of interphase mass transfer coefficient on the velocity of the solvent phase was first investigated by Frössling [23] and Ranz and Marshall [24]. Friedlander [25] derived an analytical solution for evaporation of a single droplet and cylinder whereas Dobry and Finn [26] conducted laboratory experiments to determine the impact of velocity on mass transfer. Kusik and Happel [27] developed an analytical solution for interphase mass transfer in particle beds. Williamson et al. [28] investigated the velocity dependence of the dissolution of benzoic acid spheres in a bed of particles flushed with water, which is analogous to dissolution from non-wetting fluid trapped in a packed bed to a flowing wetting fluid. These studies have consistently shown that the interphase mass transfer coefficient increases with velocity (e.g., the Ranz–Marshall equation). However, there is no agreement on a universal relationship as it is system dependent. The influence of the system complexities arising from the porous structure of the medium, on particle to fluid interphase mass transfer was first examined by Pffefer and Happel [29] who focused on the influence of the void ratio of the porous media. Wilson and Geankoplis [30] investigated the influence of axial mixing in packed beds on interphase mass transfer. Dwidedi and Upadhyay [6] reviewed the prior data on interphase mass transfer in packed beds and showed that it is inversely proportional to the void fraction of the medium. One of the earliest investigations focusing on NAPL contamination in porous media was conducted by van der Waarden et al. [31]. The authors showed that NAPL dissolution from unsaturated porous media is a much more complex process than dissolution of solid particles such as benzoic acid spheres in porous media even though there are similarities between the two. The main reason for that difference is the complex mechanics of phase formation during multiphase flow in porous media, which are discussed in detail in Section 3 of this manuscript. The uncertainty in the distribution of flow/phases due to the tortuosity of natural porous media and the complex nature of the multiphase flow field renders the use of the physical theories (e.g., the determination of the film thickness for film theory) very difficult if not impossible. As a result the mass transfer coefficients have been usually fit to data with respect to known system properties that exist in flow in porous media (e.g., porosity, uniformity coefficient, density, average velocity, etc.). The difficulty in observation of interfacial area in porous media has also led researchers to lump the kf 1 and A into a single coefficient K l (e.g., [32]). Fig. 1 illustrates a typical NAPL contamination scenario. Fig. 1(a) schematically depicts an idealized condition where the NAPL blobs are entrapped uniformly in a uniform porous medium. The dissolved concentration with distance for two different flow velocities (low/high Peclet numbers) for this idealized case is presented in Fig. 1(b) and (c). The dominant transport process in Fig. 1(b) (low Peclet number) is diffusion which results in a concentration equal to the aqueous solubility within the NAPL zone. In Fig. 1(c) (high Peclet number) the higher velocity means that the contact time between the flowing fluid and the NAPL is insufficient to yield concentrations equal to the solubility. As a result the effluent concentration is lower than the aqueous solubility at the end of the domain. However, the higher velocity condition leads to higher concentration gradients between the interface and bulk fluid, resulting in a higher interphase mass transfer rate. Fig. 1 depicts interphase mass transfer at different scales. Pore scale investigations focus on discrete pore scale parameters that are presented idealistically in Fig. 1(a). Meso-scale investigations 166 B. Agaoglu et al. / Advances in Water Resources 79 (2015) 162–194 Fig. 1. Typical NAPL contamination (based on [33]); (a) schematic of uniformly distributed NAPL blobs (shaded area), flow is from left to right; (b) concentration as a function of distance; for low Peclet number and (c) concentration as a function of distance; for high Peclet number. focus on the average dissolution behavior of some porous medium volume such as the domain presented in Fig. 1(a). Field scale investigations incorporate large scale flow attributes such as the spatial variability of the permeability distribution in the field and the corresponding bypass conditions that may occur. Initial attempts to mathematically describe interphase mass transfer for NAPL dissolution in porous media were derived based on the steady state one-dimensional advection/dispersion equation coupled with a mass release term into the flowing water similar to Eq. (1) [8,34,35]. By further assuming that advection is the dominant solute transport process [9], the solute transport equation can be rewritten as: @C q ¼ kf 2 aðC s C Þ @x ð3Þ The above equation was first proposed to back calculate an interphase mass transfer coefficient from 1D experimental data observed at the meso-scale. Here a refers to the specific interfacial area [L2/L3], or the interfacial area per unit porous medium volume whereas A in Eq. (1) refers to the absolute interfacial area between the two phases [L2]. In both approaches C refers to bulk concentration [M/L3] in the aqueous phase which is a function of location x, and the interphase mass transfer is defined as first order rate limited via kf ’s which are system unique interphase mass transfer coefficients and have units of [L/T]. The latter equation can be solved for kf 2 analytically yielding; kf 2 ¼ C eff ln 1 La Cs q ð4Þ B. Agaoglu et al. / Advances in Water Resources 79 (2015) 162–194 Here C eff refers to the concentration at the downstream end of the NAPL zone (x = L). Although Eq. (4) has a number of simplifying assumptions, it has been used to calculate the interphase mass transfer coefficients from experimental data such as column experiments (e.g., [9,36–38]) and to develop upscaled mass transfer coefficients at the field scale [38]. On the other hand, the definition of the mass transfer given in Eq. (1), has rarely been used for the interpretation of experimental data (e.g., [32] with lumped approach) because it requires the estimates of the concentration distribution and interfacial areas which are generally not available. Eq. (1) however can be used with data obtained from pore scale numerical models (e.g., pore network models; [39]) where the explicit concentration distribution in the domain, which can be used to determine the representative concentration, is available. Eq. (4) might further suggest that kf 2 is a linear function of the velocity q; however, this is not the case because, as shown in Fig. 1(b) and (c), increases in velocity may also decrease C eff . The relationship between interphase mass transfer and velocity will be further discussed in Section 3. It is evident that interphase mass transfer is controlled by the thermodynamic driving force and the physics of advective and diffusive transport. The influence of the chemical properties of the existing phases (i.e., thermodynamic driving force) is further addressed below in this section where as the mechanics of phase formation and its influence on the physics is discussed in Section 3. 2.1. Ionic strength and temperature The mass transfer of organic compounds from the non-aqueous phase to the aqueous phase can be significantly reduced by the ionic strength of the aqueous solution. Setchenow [40] presented a simple equation to predict the solubility of organic compounds for different ionic strengths. This is particularly relevant to groundwater remediation activities since groundwater contaminated by organic compounds in the vicinity of landfills and industrial areas often has high ionic strengths. Lesage and Brown [41] showed that for certain multi-component NAPLs, the solubility of all components decreased by 2–3 folds with flushing of 1 M NaCl aqueous solution. However, it is also reported that salting-in effect may sometime occur for certain organic compounds ion combinations such as naphthols in sodium fluoride solutions [42]. The temperature of the multiphase system can influence interphase mass transfer through different pathways. Increase in temperature will lead to increase in the solubility of organic species in aqueous solutions [43,44]. The diffusion coefficient of the organic species also increases with increasing temperature [45]. Consequently, the mass transfer coefficient of organic species from the non-aqueous phase to the aqueous phase would also increase with increasing temperature. Temperature may have even further implications on mass transfer because high temperature gradients in the aqueous phase might cause convection and further mixing of the aqueous phase leading to an increase in the mass transfer rate. However, this was reported to have insignificant influence within typical aquifer conditions [46]. Sorption onto the solid surface and interfacial tension between organic and aqueous phase which has direct impact on the mobilization of organic phase, are also affected by temperature change [43,47]. 2.2. Multi-component NAPLs and aging When multi-component NAPLs are present, the dissolution of individual organic compounds into the aqueous phase interfere with each other. These interactions occur both on the aqueous side and NAPL sides of the interface. The diffusion coefficient of an organic compound dissolved in the aqueous phase would depend 167 on whether or not other compounds are also present in the aqueous solution. Pfannkuch [17] pointed out that for multi-component NAPLs, the more soluble organic compounds will impact the diffusivity of the less soluble ones and, consequently, its interphase mass transfer rate. Moreover, if the change of the composition at the interface with dissolution cannot be compensated with diffusion of the organic compounds from the bulk NAPL to the interface, the interphase mass transfer rate of each compound will change since the composition of NAPL at the interface will be different than the initial conditions [48]. Furthermore, multi-component complex NAPLs often behave as non-ideal liquids due to the nonequal attraction between different molecules and molecules of the same type. This non-ideal behavior which would lead to activity coefficients that are significantly different from unity further complicates the rate of interphase mass transfer [49,50]. For multi-component NAPLs that behave as ideal solutions, solubilities are defined in terms of Raoult’s Law which expresses the effective solubility of each component as the product of the solubility of that component and the mole faction of that component in the NAPL. Luthy et al. [51] investigated the effect of NAPL aging on interphase mass transfer and showed that a single coal tar droplet develops a wrinkled, skin-like formation after as little as 3 days of direct water contact. This formation is likely to occur in case of multi-component NAPL comprising of high molecular weight organic compounds. For crude oil it was reported that these surface deformations are due to adsorption of high molecular weight polar molecules at the crude oil–water interface [52]. The interphase mass transfer rate of naphthalene from coal tar was shown to decrease after the occurrence of the skin [51]. 2.3. Cosolvents and surfactants Various chemical agents have often been considered for the enhanced recovery of NAPLs from the subsurface. The two primary recovery mechanisms are enhanced dissolution and/or mobilization of the NAPL mass. The most known agents in the literature are cosolvents and surfactants, and, to a lesser extent cyclodextrins and humic substances [53,54]. In this work we cover only the first two agents. Although both cosolvents and surfactants can enhance interphase mass transfer, the chemical processes that lead to enhancement of dissolution are different [55]. The majority of the studies cited in this section are based on batch tests and meso-scale column experiments. Cosolvents are chemical compounds that exhibit both polar and non-polar behavior which eventually leads to increased miscibility of phases. Cosolvents can primarily partition into the aqueous phase (e.g., ethanol) or into the NAPL (e.g., methanol) and in both cases the dilute solution assumption does not hold for both of the existing phases. This partitioning influences the diffusion coefficients of each component in each phase [56] and there are few methods to predict the diffusion coefficients in ternary phase systems (e.g., [57]). The change in the diffusive flux, in turn, influences interphase mass transfer. Imhoff et al. [58] demonstrated the non-linear alteration of the mass transfer coefficients with increasing cosolvent mass fraction. The equilibrium compositions of the phases can be calculated from thermodynamic principles which are incorporated in the computer codes UNIQUAC and UNIFAC [59,60]. On the other hand, the Hand’s method, which defines the phase equilibria as straight lines on a log–log plot, has been used in several multiphase multi-component models due to its easier applicability [61–63]. The Hand’s model input includes the constants that define the ternary phase diagram (i.e., height of the binodal curve, plait point) which can be acquired via batch experiments [64]. For multi-component NAPLs, it was shown that the increase in ethanol content of the aqueous solution leads to enhanced solubility of each component but at different proportions 168 B. Agaoglu et al. / Advances in Water Resources 79 (2015) 162–194 [54]. A simple model to predict the equilibrium solubility of organic compounds in aqueous phase in the presence of cosolvents was also developed by Yalkowsky [65] whereby the enhanced solubility is proportional to the octanol–water partition coefficient of the organic compounds. In NAPL remediation the use of cosolvents also influences the interfacial tension, viscosity, density and microbial activity [66–68]. As will be discussed in Section 3, alteration of these properties may also have a significant indirect impact on the interphase mass transfer process. Surfactants are surface active materials that enhance the dissolution of organic compounds by encapsulating them in micelles which exhibit a polar outer and a non-polar inner surface [69,70]. Several mechanisms have been proposed to represent the miceller solubilization process such as: oil diffusion model; micelle diffusion–adsorption model; micelle disassociation–reformation model [71–73]. Unlike the log-normal dependency of solubility values of organic compounds to cosolvent amount as proposed by Yalkowsky [65], surfactant solubilization is linearly dependent on the surfactant amount. Moreover, the enhanced solubility does not exhibit any preferences between organic compounds from a multi-component NAPL, whereas cosolvents dissolve the compounds with higher octanol–water partition coefficient at higher rates [54,55,72]. Similar to cosolvents, the phase behavior of the multiphase system can be defined through ternary phase diagrams for surfactant solutions [74]. Surfactant micelles do not form unless a critical surfactant concentration is achieved (referred to as the critical micelle concentration, CMC). This means that the surfactant has minimal effect on the solubility below the CMC [75]. The addition of surfactant at concentrations below its CMC however yields substantial decline in interfacial tension [76]. It has also been shown that the interphase mass transfer coefficient decreases with increasing surfactant concentration above CMC [72,73,75,76]. Grimberg et al. [72] proposed that this is due to the retardation of the diffusion of solute saturated micelles from the NAPL–aqueous interface towards the bulk solution. This phenomenon was attributed to the accumulation of surfactants near the NAPL–aqueous phase interface retarding the diffusion of micelles containing organic compounds [77]. The partitioning of surfactant micelles into the NAPL has been also shown to influence the enhanced interphase mass flux [78,79]. However Lee et al. [80] proposed that the impact of partitioning into the NAPL be taken into account together with the solubilization capacity. Effect of pH and ionic strength on surfactant enhanced dissolution has also been investigated for certain surfactant, organic compound couples. For example Park [81] found that lower pH leads to higher solubilization of pentachlorophenol (PCP) due to the greater hydrophobicity of PCP, whereas increased ionic strength leads to a decrease in solubilization. investigations are applied for each computational node (i.e., grid block) (e.g., [38,84]). Representative elementary volume (REV) is the volume over which a pore scale phenomena, in this case interphase mass transfer, is considered uniform. In practice the grid discretization of REV-based numerical models forces REV to be defined at the size of grid block. Grid block size can range from a size of 1 cm3 to model a 1 L column to a size of 10,000 cm3 to model 1000 m3 field scale problem. Hence, an REV in one study may be even greater than the whole domain of another study. This is further elaborated in the following sections. 3.1. Pore scale The investigation of interphase mass transfer at the pore scale requires information on the spatial distribution of the interfacial area and relevant flow paths within individual pores. Recent experimental techniques have identified three different domains of interfacial area when two phases are present in a porous domain: (1) capillary (or meniscus) area; (2) thin film area; (3) grain surface roughness area (Fig. 2) [85–87]. The capillary area is the interface that forms at the meniscus of a pore whereas the film area is formed when the wetting phase covering the soil particle is surrounded by a non-wetting phase. A third domain, referred to as the grain surface roughness area, forms when the wetting phase saturation approaches zero such that the grain surfaces are covered with a very thin layer of wetting fluid that follows the roughness (troughs and ridges) present on surfaces. Interphase mass transfer can potentially occur through these domains in all of the flow configurations described in Table 1. Fig. 3 demonstrates the possible flow pathways that will lead to interphase mass transfer in two different cases of wettability. In the first case, where the mobile phase constitutes the wetting phase (Fig. 3(a); ‘‘INtoMW’’), both meniscus and film areas can contribute to interphase mass transfer. The impact of film interfacial area on total interphase mass transfer is proportional to the flow rate through the corner which is a function of the 4th power of the corner equivalent radius (Poisuielle’s flow). The reduction of the non-wetting blob via dissolution induces shrinking. As 3. Mechanics of interphase mass transfer In this section we review studies that investigated interphase mass transfer at the pore, meso and field scales. The discussion of each scale includes also the influence of different flow configurations, which were described in Table 1. Pore scale studies (e.g., micro-models, pore network models, [82,83]) have been used to examine the implications of key parameters at the meso-scale. This is because pore scale investigations enable the control of each parameter separately (e.g., individual pore throats, pore sizes, etc.). On the other hand interphase mass transfer studies at meso-scale have been used to evaluate the mass transfer at the field scale. The most common method for that has been to simulate field scale scenarios with use of REV-based numerical models where the interphase mass transfer formulations developed from meso-scale Fig. 2. Schematic of the meniscus and film area around grain particles. Electronmicroscopy visualization of surface roughness can be found in [85]. B. Agaoglu et al. / Advances in Water Resources 79 (2015) 162–194 169 Fig. 3. Shrinking of a non-wetting blob (a) and wetting blob (b) via dissolution. Numbers represent the progression of dissolution in time. L and w are the film length and width respectively [82,89]. dissolution continues, the length of the corner diminishes. However, the equivalent radius of the corner is not significantly influenced by this shrinking because the non-wetting phase approaches a more stable spherical shape with on-going dissolution. Once the blob takes the shape of a sphere, further dissolution leads to the reduction of the equivalent radius of the corner (position 5 in Fig. 3(a)). This means that the flow velocity will increase in the film layer only after the blob takes a spherical shape and decreases in radius. To date there is no experimental study that has presented the flow velocities occurring through film layers due to experimental difficulties at this small scale. Even though Poisuielle’s equation can provide an estimate of the flow velocity, the influence of surface roughness of grains on the flow cannot be accurately determined. Fig. 3(b) depicts the case when the mobile phase constitutes the non-wetting phase (i.e., ‘‘IWtoMN’’). Under these conditions, most of the interphase mass transfer occurs from the film area. As the mobile phase transports the dissolved solutes from the interfacial film area, this film area is depleted and the wetting phase residing in the meniscus region spreads itself to the film area. This leads eventually to the withdrawal of the meniscus (from position 3 to 1 in Fig. 3(b)) until the pressure of the mobile non-wetting phase overcomes the capillary pressure and breaks the wetting phase. At this instant the equivalent radius of the pore throat center (where the mobile phase starts to flow through) increases significantly because the wetting phase spreads itself over the grain surface to form a more thermodynamically suitable shape. The formations represented in Fig. 3 have been shown qualitatively in micro-model experiments conducted by Sahloul et al. [82]. Zhao and Ioannidis [88] used pore network modeling to demonstrate that the continuity of the NAPL film, which is a function of the wettability of the medium, is the most influential factor on the interphase mass transfer rate for NAPL wet-media. The benefit of defining the flow configurations as in Table 1 is clearly observed in the discussion presented above along with Fig. 3. For example in ‘‘IWtoMN’’ (i.e., Fig. 3(b)) the dominant interfacial domain in terms of interphase mass transfer is the film area as it is in contact with mobile water regions whereas for ‘‘INtoMW’’ (i.e., Fig. 3(a)) mass transfer is mostly through the meniscus area as the flow rate in the film region is very low and therefore the contribution of these regions to interphase mass transfer is limited although some studies suggest that the film area may also contribute to the interphase mass transfer especially at high pressure gradients [83,90]. It should be noted to date there is no quantitative data on the thickness of film layer presented in Fig. 3(a) for real porous media which is the parameter governing the flow rate within this film layer. However, X-ray microtomography images taken by Schnaar and Brusseau [37] showed that the film thickness appearing in ‘‘INtoMW’’ (i.e., Fig. 3(a)) is extremely small. Furthermore, the curvature and surface roughness of the grain particles may even further impede the flow within the film layer. In accordance with these observations, Bradford et al. [89] experimentally showed that the interphase mass transfer in NAPL-wet media is much greater than the water-wet media providing further evidence that the flow configurations depicted in Fig. 3 can be useful in identifying which components of the interfacial areas contribute the most to interphase mass transfer ‘‘MNtoIW’’ can be visualized from Fig. 3(b) if dissolution occurred from mobile phase to immobile wetting phase. ‘‘MNtoIW’’ appears in high pressure CO2 injection (i.e., supercritical phase) applications [91,92]. ‘‘MWtoIN’’ can be observed also from Fig. 3(a) if the dissolution were to occur from the mobile phase to the immobile non-wetting blob. ‘‘MWtoIN’’ was investigated within gas exsolution applications where for example the injected wetting phase (water) is super saturated with CO2, leading to the expansion of the non-wetting CO2 gas bubble (blob) with dissolution from the water phase to the gas phase [93,94]. In both of these applications mobile phase is the dissolving one. The implication of this condition is that now the velocity of the mobile phase is responsible for the magnitude of interphase mass transfer as it supplies the dissolving compounds to the interface. This is different than the first two configurations (‘‘INtoMW’’ and ‘‘IWtoMN’’) where the velocity of the solvent phase is responsible for the magnitude of interphase mass transfer. This is the reason why both ‘‘MNtoIW’’ and ‘‘MWtoIN’’ are categorized as new flow configurations even though they have the same influential interfacial domains with ‘‘INtoMW’’ and ‘‘IWtoMN’’ respectively. It must be noted that both ‘‘MNtoIW’’ and ‘‘MWtoIN’’ have not been studied in detail in terms of dissolution in porous media as these subjects are relatively new in the field. ‘‘CHEM’’ is observed in enhanced NAPL remediation/oil recovery applications and depending on the wettability of the system either Fig. 3(a) or (b) will appear. However at very high amounts of chemical agents the phases become fully miscible in both conditions. 170 B. Agaoglu et al. / Advances in Water Resources 79 (2015) 162–194 In two-phase flow both components dissolve to some extent into the counter phases. However, in most real applications one of the components is generally of significance and therefore, the dissolution of the other component is ignored if it does not have significant impact on system behavior (e.g., the resistance of dissolving solutes in the dissolving phase is negligible due to very low amount of mutual dissolution; dilute solution). For example in NAPL remediation investigations, the dissolution of water molecules into NAPL is ignored in experimental studies as well as in numerical studies (e.g., pore network models, single phase REV models, semi-multiphase REV models) as it does not influence fate of organic compounds in the system. However, especially in enhanced NAPL remediation applications, even though the fate of a certain compound is sought (i.e., fate of organic compound), interphase mass transfer of all components into both phases should be considered as the composition of both phases may significantly change the system behavior in some instances [21,60,63]. 3.2. Meso scale Most of the data and observations on interphase mass transfer in porous media are derived from the column or tank experiments. The complexity of flow in porous media prevents the use of physical models such as film layer or boundary layer theory and as a result many of the experimental and numerical investigations developed correlations between interphase mass transfer coefficients and system parameters. The information gathered at this scale has been used in numerical models defined in terms of Representative Elementary Volumes (REV). Earlier investigations were based on derivation of interphase mass transfer coefficients from column and tank experiments with use of Eq. (4) or similar approaches (e.g., [8,9,35,95]). As noted in Section 2, the unavailability of the interfacial area was compensated by lumping the area term into the mass transfer rate coefficient, Kl. Recent advances in non-invasive imaging technologies have in some instances led to relaxation of this assumption, enabling the direct use of Eq. (4) (e.g., [37]). In addition investigations based on pore network models also permit the use of Eq. (1) because these models can provide the concentration distribution within the porous media domain (e.g., [39]). After the mass transfer coefficient is determined from experimental data, a Sherwood number (Sh ¼ kf dm =D) is defined and correlated against the system properties such as porosity, velocity, and uniformity coefficient. The modified Sherwood number 0 2 (Sh ¼ K l dm =D) has been used in cases when accurate estimates of the interfacial area are not available. Below are two examples of these Sherwood formulations developed by Pfannkuch [17] (Eq. (5)) and Nambi and Powers [36] (Eq. (6)). Sh ¼ 0:55 þ 0:025Pe1:5 ð5Þ Re0:61 Sh ¼ 37:2S1:24 n ð6Þ 0 Pfannkuch developed the above Sherwood formulation (Eq. (5)) from the results of column experiments conducted by Hoffmann [96] and by using the Bowman’s equation as a base. Hoffmann’s experimental set up consists of a stagnant oil phase located in the upper half of a horizontal column filled with glass beads with water flowing horizontally along the flat contact area. Pfannkuch developed the Sherwood formulation of each experiment having differm ent velocity from Sh ¼ Ed . E refers to interphase mass transfer Cs D rate divided by the interfacial area, (which is rather clearly defined in [96]), kf is calculated, by ignoring the concentration of solute in the column along the oil–water contact area due to lack of direct measurements, with; kf ¼ @M1 @t A Cs (i.e., use of Eq. (1) assuming C = 0). Furthermore, in the original German text Hoffmann states that the solubility of domestic fuel oil in his experiments (i.e., Heiz Öl) is between 103 and 105 mg/m3 and Pfannkuch assumed this value in his calculation to be 104 mg/m3 which is another source of uncertainty to the calculated mass transfer coefficient. The main advantage of this set up was that it allowed for the estimation of the oil–water contact area, whereas in other situations the interfacial area would be a complex function of fluid and porous media properties that is difficult to measure or estimate accurately. This feature of the Pfannkuch’s Sherwood expression led to its use in several subsequent studies [33,97,98]. It is important to note that in real systems the complex NAPL distribution and complex flow fields may actually lead to different concentration profiles than those proposed for these idealized conditions. The influence of these factors on interphase mass transfer coefficient is discussed in Section 3.4. The modified Sherwood formulation of Nambi and Powers [36] (Eq. (6)) was developed from the results of small 2D experiments using Eq. (4), where the interphase mass transfer coefficient is lumped with interfacial area (i.e., mass transfer rate coefficient), and multiple regression analysis was used [9,36]. A number of studies have developed similar Sherwood formulations (e.g., [8,32,99–101]) from meso-scale experiments whereby the NAPL was mostly in the form of distributed blobs while Nambi and Powers [36] is the only study that was based on pool NAPL configuration. In recent years, several studies have also shown that the use of the modified Sherwood formulations may lead to contradicting results when used in numerical models even though they have the same input parameters [83,84,90,98,102–108]. This is mostly attributed to differences in the specific experimental conditions used in developing these models. The Sherwood formulations which consider interfacial area (e.g., [17] Eq. (5)), are also subject to this consideration. For example Pfannkuch developed one other Sherwood formulation from the results of the experiments of Zilliox [109] which had median grain size of 2000 lm, compared to 500 lm for the Hofmann data. The calculated Sherwood numbers were found to be about 10 to 100 times greater than the Sherwood numbers calculated from Hoffmann’s experiments for the same Peclet numbers. Moreover these 0 Sherwood formulations (both Sh and Sh; e.g., [8,17,105,110]) have been mostly developed from the results of column experiments. In many instances these formulations are used in numerical models with node (REV) sizes that are significantly different from the experimental scale that they have been developed which may also lead to erroneous results. This is a result of the fact that the parameters used to calculate mass transfer coefficients via Eqs. (1) and (4) change at different orders when system size changes. For example, a large grid block having a uniform NAPL blob distribution would lead to near equilibrium conditions even for high velocities whereas a small grid having the same type of NAPL distribution and flow velocity would result in rate limited behavior. If one chooses a Sherwood formulation developed from an experiment with size and conditions similar to the small grid block described above and uses it in a numerical model with larger grid sizes, this may result in under-estimation of interphase mass transfer. Furthermore, because the variability at the sub grid scale cannot be accounted for, the larger grid size would also result in a higher level of approximation. Because most of the Sherwood formulations are derived from column experiments with uniform NAPL blob distributions, accurate REV-based numerical simulations requires compatibility with this condition in the grids that contain NAPL. If this condition is not met interfacial area and flow paths will cause different dissolution results B. Agaoglu et al. / Advances in Water Resources 79 (2015) 162–194 than the predicted one. Consequently, the adoption of Sherwood formulations to real applications should be done with caution. 3.3. Field scale The salient difference observed at the field scale compared to meso-scale experiments is the complexity of the NAPL spatial distribution and the groundwater flow field. In the subsurface NAPLs can be retained as either residual, pool or both. For NAPL–water system residual saturation of NAPL refers to the case where NAPL cannot further decrease in saturation with applied hydraulic pressure (e.g., [111]) and exist in the form of single or multi-pore blobs. In pool configuration the saturation is greater than the residual yet hydraulic pressure is not enough to initiate mobilization. Such configuration is especially observed when DNAPL accumulates over impermeable media after downward leakage or when LNAPL is formed over the water table. As in other subsurface applications (e.g., groundwater extraction) the most important factor that controls interphase mass transfer at the field scale is the permeability spatial distribution at the site [84]. It has been reported that poorly accessible NAPL mass is the primary reason for the long term mass fluxes observed at field [112]. Many studies have reported that the interphase mass flux from residual NAPLs is significantly higher than pool NAPLs (e.g., [113,114]). Anderson et al. [115] showed that flow through a residual DNAPL zone leads to equilibrium concentrations which decrease with reduction in NAPL saturation. In residual NAPL regions, the interphase mass transfer rate can also decrease due to dissolution fingering. Imhoff and Miller [116] reported that the non-uniform distribution of residual NAPL leading to variation in flow velocity between pores is the primary reason for this phenomenon. On the other hand it has been also shown that such fingers are limited in length due to transverse dispersion [117–119]. For pool NAPL configurations Brusseau et al. [120] demonstrated that dilution by clean water bypassing NAPL contaminated regions is a major reason for the reduced effluent concentrations. Spatial variability in intrinsic permeability is the primary reason for the flow diversion and dilution. Variations in the relative permeability of aqueous phase in the NAPL contaminated area can also cause flow diversion and low effluent concentrations. Nambi and Powers [121] showed that for pooled NAPL in a heterogeneous media water flow lines are diverged around the surface of the contaminated region, initially inducing slower interphase mass transfer rates. At later stages when the water starts to penetrate into the contaminated area due to increase in aqueous phase relative permeability, the mass transfer rate increased substantially. Similar results were reported by Marble et al. [107] and Zhang et al. [122]. Page et al. [123] also showed that the reduction of effluent concentrations due to increased velocity is not significant compared to bypassing of NAPL. Brusseau et al. [112] reported that NAPL existing in poorly accessible porous media (i.e., lower hydraulic conductivity regions) regardless of being pool or residual leads to low effluent concentrations and long term mass fluxes. These observations were mostly obtained from tank experiments with localized NAPL pools and are discussed here because of the similarity of the initial localized NAPL configurations used in these studies with those typically observed in field applications. It is important to note that the method used for emplacing the entrapped NAPL within the porous medium is also of particular importance [9,32]. Mixing the sand with the contaminant phase which has been performed in a number of studies (such as [8,120,124]), rather than injecting the NAPL, might result in different dissolution patterns. Estimates of interphase mass transfer at the field scale have been mostly obtained from modeling studies (e.g., [38]). As also 171 observed in large scale tank experiments, these numerical models have shown that flow bypassing due to permeability and relative permeability differences is the major reason for low effluent concentrations. Using an REV-based numerical model for a domain size of 10 m 10 m 10 m Parker and Park [38] showed that the interphase mass transfer rate coefficients determined from fieldscale applications (i.e., up-scaled mass transfer rate coefficients) tend to be much lower than those derived from meso-scale experiments. The difference is attributed mostly to the higher degree of bypassing of the NAPL zone at the field scale compared to mesoscale. The authors also reported that reducing the domain size led to slight increase in the calculated interphase mass transfer coefficients. Groundwater flow velocity is another factor influencing interphase mass transfer; however, this dependency differs from the one observed at the meso-scale. Parker and Park [38] stated that the mass transfer coefficient of their field scale simulation is a function of velocity with a power of 1 whereas meso-scale studies reported a power of about 0.7 (e.g., [8,9]). These differences further show that the formulations of interphase mass transfer coefficients are dependent on the system conditions (e.g., permeability distributions, saturation, etc.) that they are developed from. The complex features of the field scale NAPL distributions are also influential on the field scale dissolution behavior. Zhu et al. [125] stated that early high concentrations are related to dissolution of residual NAPLs as they enable contact with the flow through the NAPL region. Following a steep decrease in effluent concentrations, once residual NAPLs are depleted, a period of low effluent concentrations initiates. The invalidity of interphase mass transfer expressions derived from meso-scale studies to field applications led some researchers to develop analytical interphase mass transfer formulations that are specific for field scale NAPL dissolution problems (e.g., [11,38,125]). For example Christ et al. [11] proposed that ganglia NAPL mass to pool NAPL mass ratio can be used to predict effluent concentrations in field. The authors used the results of a field scale numerical model, whose three dimensional domain was 8 m 8 m 9 m, to develop such an analytical formulation. These formulations are discussed in Section 4. Because these analytical formulations are determined from the results of numerical models, their validity is directly related to the accuracy of the numerical simulations. A particularly relevant study is that of Maji and Sudicky [84] who simulated a field scale NAPL dissolution problem (domain size of 10 m 16 m 6.7 m with a grid block size of 0.25 m 0.50 m 0.10 m) with 10 different mass transfer formulations. It was demonstrated that the time needed for full depletion of NAPL mass varied from 562 days (local equilibrium approach) to 48,085 days (formulation of [105]) between different mass transfer formulations. On the contrary Parker and Park [38] compared the depletion time of full NAPL mass for two formulations (local equilibrium and model of [8]) and stated that the difference is not significant and local equilibrium approach might be usable. In a related issue Nambi and Powers [36,121] showed that saturation above 35% would lead to local equilibrium conditions for their meso-scale NAPL dissolution experiments meaning that use of local equilibrium in a grid block of an REV-based numerical model might be a viable option in such cases. Yet this threshold saturation value that results in equilibrium conditions in a grid block is variable for different sizes of grid blocks and different types of porous media and NAPL distribution topologies. Decreasing the grid block size would decrease this threshold provided the grid block size does not decrease below possible minimum REV size value. Yet no data were reported on the dependence of this threshold saturation on grid block size for different velocities, and porous medium properties. The accuracy of numerical models with different mass transfer formulations may be checked against small 3D dissolution experiments. 172 B. Agaoglu et al. / Advances in Water Resources 79 (2015) 162–194 However it should also be noted that discrepancies between observed and simulated NAPL dissolution is likely to increase at the field compared to smaller 3D laboratory experiments due to use of the larger grid block and greater NAPL amount at the field. In field applications many contaminated sites are decades old and therefore it is likely that the residual NAPL residing in accessible media is already depleted. Therefore, assuming that NAPL containing REV’s are at pool conditions and they can be represented with local equilibrium assumption, seems to be in some instances viable. But the influence of dissolution once the saturation decreases below the threshold NAPL saturation is uncertain. As the influence of residual NAPL in accessible media can be neglected and fast dissolution can be assumed, it is also likely that in poorly accessible media the residual NAPL will lead to slow dissolution and tailing of effluent concentrations. The difference in dissolution behavior of residual NAPL in high and low hydraulic conductivity media is further discussed under Section 3.4.5. It should be emphasized again that field scale NAPL dissolution, like other real subsurface applications, is strongly controlled by permeability differences in subsurface and accurate description of field characteristics is essential for any prediction effort. Apart from the above discussion which pertains to ‘‘INtoMW’’, the use of surfactants and cosolvents or their mixtures for enhanced remediation purposes (i.e., ‘‘CHEM’’) can drastically alter the original multiphase flow conditions at field scale. Application of cosolvents/surfactants may culminate into increased miscibility whereby the dilute solution assumption is not valid anymore [63]. The density alteration of the existing phases may lead to flow stratification which can lead to bypassing and reduced overall mass transfer [126,127]. The alteration of viscosity might also induce development of fingers and influence system stability which can also lead to bypassing [128–130]. Reduction of interfacial tension may also influence the system mobility which, as reported by Agaoglu et al. [54], can significantly impact interphase mass transfer rate. Zhong et al. [131] observed the occurrence of dissolution fingering during surfactant flushing. The authors also showed that once a finger reaches the effluent end, a substantial decrease of the interphase mass transfer develops, implying significant bypassing of NAPL by surfactant solution. This can significantly influence the performance of field scale enhanced remediation activities. ‘‘IWtoMN’’ is encountered in the field as NAPL dissolution in NAPL-wet media or air sparging applications. In comparison to former, air sparging shows very unstable flow conditions due to high viscosity difference between the two phases. The injected air was shown to follow three different flow patterns, namely: bubble flow, channel flow and chamber flow. In a tank of uniformly packed glass beads Brooks et al. [132] demonstrated that channel flow occurs within glass beads smaller than or equal to 1–2 mm size whereas bigger glass beads lead to bubble flow (e.g., flow of discrete bubbles). Peterson et al. [133] pointed to chamber flow which can be identified by a significant horizontal component in case of a porous media of very fine to fine sand (d50 < 0.21 mm). Selker et al. [134] showed that a near injection region is formed after the injection of air where viscous forces dominate the buoyancy forces and lead to lateral expansion of the gas channels. Beyond this region, the injected air starts to move in a parabolic pathway upwards. The authors also observed direct vertical or near vertical movement of air bubbles in coarser sand and for typical injection pressures. The above literature focused on the flow paths due to air sparging in air–water systems with no NAPL present in the medium. Heron et al. [135] showed that air sparging is not effective for pooled DNAPLs since both DNAPL pools and the air injection point typically reside at the bottom of the aquifer and horizontal movement of air channels is limited. Berkey et al. [136] reported substantial short circuiting of air channels in fine sand (silty sand). Thomson et al. [137] and Waduga et al. [138] reported that the prediction of accurate recovery rates necessitates the accurate topology of the air channels and NAPL. The reasons for unstable movement of air are discussed in more detail in Section 3.4.3. A part from that air sparging can lead to decrease in relative permeability of the groundwater at the injection region in case of a moving contaminant plume [139]. The system under consideration in the above cited works relating to NAPL recovery was entrapped NAPL (in form of pools and entrapped blobs) in water saturated porous medium, meaning that the interphase mass transfer occurred from both the air–water interface and the air–NAPL interface. Steam injection is another groundwater remediation application that falls under ‘‘IWtoMN’’. Steam injection was first proposed by Hunt et al. [140] whereby the NAPL is vaporized during the flow of steam in the aquifer. Once the steam front reaches the NAPL zone, it instigates interphase mass transfer from the NAPL to the steam. Reverse mass transfer can occur in the form of re-condensation and accumulation of the vaporized NAPL at the steam front due to heat loss [141–143]. Several factors have been shown to impact NAPL mass transfer and recovery including, the presence of air with the steam [112], injection pressure [113] and field conditions such as soil type and well distance [144–146]. 3.4. Factors influencing interphase mass transfer The preceding section highlighted some of the main features of interphase mass transfer at the pore, meso and field scales. In this section, we review and discuss the influence of individual parameters on interphase mass transfer. The information presented here is mostly inferred from controlled pore and meso-scale experiments. 3.4.1. Interfacial area Interfacial area is one of the most influential parameter on interphase mass transfer. As depicted in Fig. 1, there are three different domains that contribute to interphase mass transfer, namely: (1) capillary (or meniscus) area, (2) thin film area, (3) grain surface roughness area. Because the grain surface roughness area occurs when saturation of wetting phase approaches zero, it is relevant to the flow configuration ‘‘IWtoMN’’. The meniscus and thin film area dominate the flow configuration ‘‘INtoMW’’. A substantial body of literature investigated the development of interfacial areas in two phase flow in porous media. The mostly used techniques to determine interfacial domains have been interfacial tracers and non-invasive imaging technologies [85–87]. Because of the resolution limits of the latter which is on the order micrometers, the third domain (grain surface roughness), which is on the order of nanoscales, has only been studied with interfacial tracers [147,148]. Non-invasive technologies such as X-ray microtomography have been shown to be capable of distinguishing between the meniscus and film areas where as interfacial tracers provide a measure of the overall interfacial area. It has been also demonstrated that the gaseous phase interfacial tracers have more accessibility to the interfaces that occur in porous media compared to aqueous tracers [86,149]. The interfacial area between water and air has been reported to increase with decreasing saturation of the wetting phase (aqueous phase) [150,151] and decreasing mean grain diameter [147,152]. The change of interfacial area was also shown to be greater for poorly sorted media [148]. It has been found that the data inferred from synchrotron X-ray micro-tomography exhibits a negative linear relationship between water saturation and interfacial area whereas interfacial tracer tests indicate a negative exponential relationship especially when the water saturation approaches zero. This was attributed to the inability of micro-tomography to measure the very thin layer of water forming on the surface roughness of grain particles [153,154]. Distinction between meniscus and thin 173 B. Agaoglu et al. / Advances in Water Resources 79 (2015) 162–194 Fig. 4. Variation of different domains of the interfacial area with wetting phase saturation. film area made by non-invasive methods demonstrated that meniscus area diminishes to zero when wetting saturation decreases to a certain saturation depending on the type of porous medium [111,153]. Generic relationships between the interfacial area of different domains as a function of wetting phase saturation (during the first drainage) are presented in Fig. 4. The interfacial area between air and water was also found to be influenced by the existence of surfactant molecules when present in the aqueous phase. It has been reported that increase in surfactants content alter the interfacial tension between air and water which eventually affects the saturation profile of the system in a non-linear manner [154]. Interfacial area was also investigated in fractured media; however, it was reported that in fractured media there is no correlation between interfacial area, saturation and aperture/aperture ratio [155]. The evolution of interfacial area during the dissolution of residual NAPL was first demonstrated via MRI for a column filled with glass beads. It was shown that interfacial area decreases linearly with saturation which decreases with dissolution [156,157]. This was also shown to be the case for real porous media [37]. It is also argued that the meniscus interfacial area is the principal domain for interphase mass transfer. However the ratio of meniscus interfacial area to film area (domain 1 to domain 2) increased as dissolution proceeded, reaching 50% once 95% of the NAPL is dissolved. The authors attributed this increase to conversion of film area to meniscus area as the dissolution shrinks the NAPL blobs. This is also depicted in Fig. 3(a) which shows that the ratio of the meniscus area increasing from position 1 to position 5. Several studies attempted to develop functional relationships for the prediction of interfacial area for an air–water system. Peng and Brusseau [148] used a Van Genuchten type equation to predict the total interfacial area (domain 1, 2 and 3) normalized by the total porous medium volume: m Aia ¼ s 1 þ ðaðSw Sm ÞÞn ð7Þ Aia ðcm1 Þ ¼ SA½ð0:9112ÞSw þ 0:9031 ð8Þ The authors found good agreement between the developed relationship and the experimental data of a particular study. Costanza-Robinson et al. [87] developed another functional relationship which focused on the prediction of interfacial area of domains 1 and 2 excluding domain 3. The authors reported that when the air–water interfacial area in a particular porous media is normalized by the geometric surface area of that medium, the relation between this normalized interfacial area and water saturation is independent of porous medium type: Consequently, the developed relationship can be used to determine the air–water interfacial area if the water saturation of a sample and the mean grain diameter of the medium (to be used in the determination of the geometric surface area, SA) are known. It was reported that this predictive method resulted in good agreement with experiments conducted with glass bead media whereas for natural media some inconsistencies emerged. Thermodynamics has been also used to investigate the interfacial area in porous media [158]. Morrow [159] presented a method for the estimation of interfacial areas through thermodynamic-based approach which can also be found in Bradford and Leij [160]. This method relates the external work applied on a closed system of porous medium comprising of two phases (which constitutes the area under Pc–Sw) to the existing interfacial area through a constant of integration (CNW) which refers to the influence of existing NAPL blobs. CNW is zero for the first drainage/imbibition and nonzero after the formation of first residual saturation: þ C NW ¼ rNW cosð/sNW ÞAsN SNW þ rNW ANW SNW UNW SNW W W W ð9Þ Dalla et al. [161] adopted the pore scale model of Hilpert and Miller [162], where the domain comprised of spherical particles, to investigate the development of interfacial area and to compare the results against the thermodynamics-based predicted interfacial area developed by Bradford and Leij [160] and the interfacial area measured by Kim et al. [150] using aqueous tracer tests. The authors showed that the thermodynamics approach of Bradford and Leij [160] predicted larger interfacial area compared to the area attained through their pore network model and aqueous tracer tests conducted by Kim et al. [150]. The dissipation of the mechanical work applied to the system (i.e., conversion of only a portion of the mechanical energy into surface energy) which is not accounted for in Bradford’s method [160], is stated to be the main reason for the higher interfacial area calculated by this method. Another reason for the low interfacial areas determined by the aqueous interfacial tracer test compared to results of Bradford and Leij [160] is reported to be the partial measuring of the interfacial areas by aqueous tracers because of the noncontinuity of the aqueous phase especially at reduced water saturations. Dobson et al. [163] showed that the Bradford and Leij [160] method culminates into the prediction of lower interfacial areas compared to experimental ones measured by the aqueous tracer test on a NAPL–water system whereby the NAPL was at residual saturation. It must be noted that in Dobson et al. [163] the aqueous phase was continuous, meaning that all interfacial areas are measured with aqueous interfacial tracers and therefore the method of Bradford and Leij [160] is expected to better fit the measured data. However, the method of Bradford, along with a similar method developed by Oostrom et al. [164] produced smaller values compared to the measured ones which lead the authors to develop a modified version of this method that considers the area under the drainage curve minus the area under the imbibition curve. The authors found good agreement between the predicted and the measured interfacial areas for different type of porous medium having different median grain size. Another predictive method derived from thermodynamics principles was developed by Grant and Gerhard [165] whereby the focus was on the meniscus interfacial area considering energy dissipation. Similar to the modified version of the Bradford method which is presented in Dobson et al. [163], this method also considers the areas under both the drainage and imbibition curves. However, this method enables the prediction of the interfacial area along the saturation curve contrary to the Dobson et al. [163] who investigated only the residual saturation. Moreover, Grant’s method [165] accounts for the effect of hysteresis in the Pc–S 174 B. Agaoglu et al. / Advances in Water Resources 79 (2015) 162–194 relationship and therefore enables the prediction of interfacial area at subsequent drainage and imbibition cycles: awn SiW ¼ W Siw Ed / ½UNW ðSw ÞD;I;D rNW ð10Þ The authors also made the assumption that capillary (meniscus) interfacial area decreases linearly with reduction of non-wetting phase saturation induced by dissolution. This was tested in a semimultiphase flow model in a companion paper [33]. It is interesting to note that the authors made the correct assumption independently of the study of Schnaar and Brusseau [37] who had presented the linear relationship between capillary (meniscus) interfacial area and saturation induced by dissolution via micro-tomography. Since thermodynamics approaches give information about the total interfacial area (i.e., domain 1 and 2) the authors used an equation, which is extracted from the results of the pore scale model of Dalla et al. [161], that gives the ratio of capillary (meniscus) interfacial area (referred to in this paper as specific interfacial area) to total interfacial area at different water saturations. However, the influence of hysteresis on this ratio is not considered in this equation. The authors set the energy dissipation parameter (Ed) to 0.21 after testing it against the interfacial area measured by Brusseau et al. [166]. It is noteworthy that the single equation representing the ratio of capillary interfacial area to total interfacial area and single valued energy dissipation chosen by the authors (i.e., Ed = 0.21) may be system dependent; further experimentation would be required to establish its general validity. The assumption that the dissolution occurs only through the capillary interfacial area is also subject to discussion which will be considered further in this review under the pore network modeling section. Occurrence of interfaces between three fluid phases in porous media is an even more complex process involving spreading coefficients [167]. Consequently, there are only a handful of studies which attempted to investigate this problem, mostly through micro-models (e.g., [167,168]). Schaefer et al. [151] measured the interfacial area between water and oil using interfacial tracers for a porous medium containing the three phases water, air and oil. The authors reported that the interfacial area between water and oil increased linearly with decreasing water saturation yet the measured values were half of the values that were measured for the two phase experiment (water–oil) at the same conditions. This result appears to contradict findings of earlier micro-model work [167,168]. 3.4.2. Hysteresis NAPL spatial distribution in the subsurface and, consequently the resulting interfacial area and flow paths have major implications on the interphase mass transfer as discussed above. The distribution of NAPL mass is typically represented in terms of constitutive relationships which express the capillary pressure and relative permeability of the phases in terms of fluid saturations. A complication that appears in porous media relating to these relationships is hysteresis. Pc–S relationship of the initial drainage was shown to be different than the following imbibition, as also subsequent drainage and imbibition cycles were shown to be distinct. This hysteric behavior is argued to be a result of lumping of too many system parameters [169]. Using micro tomography Cullingan et al. [111] also showed that the saturation – interfacial area (S–Ai) relationship exhibits hysteric behavior. As accurate interphase mass transfer predictions require accurate determination of NAPL spatial distribution and interfacial area, the hysteresis should also be adequately accounted for. Starting from thermodynamics principles and taking into account the forces acting on fluid interfaces, Hassanizadeh and Gray [169] proposed that the incorporation of interfacial areas in Pc–S relationship would likely reduce hysteresis. Initial investigations of this hypothesis were based on pore network models due to the difficulty arising in experimental observation of the interfacial areas. Reeves and Celia [170] demonstrated that the hysteric behavior can be discarded when the Pc–S relation is extended to include the Ai existing between two phases and that the Pc–S–Ai surface is unique for a particular fluid/fluid/medium system. It should be noted that the pore network model did not include snap off mechanisms. Held and Celia [171] expanded this model to embody several snap off mechanisms providing further data supporting the hypothesis of Hassanizadeh and Gray [169]. Another pore network model which confers further evidence in support of this hypothesis was developed by Joekar-Niasar et al. [172]. The authors showed that the trapping assumptions incorporated in the pore network models significantly influence the Pc–S–Ai relationship. Moreover, it was reported that the structure of the pore network influences this relationship such that one of the pore network models used in that study (namely, the tube model) could not reproduce hysteresis in the Pc–S curve during imbibition. The first experimental confirmation of dependence of interfacial area on the Pc–S curve with use of micro model was presented by Cheng et al. [173] who showed that the Pc–S–Ai surface is unique. Chen et al. [174] presented further evidence that the Pc–S–Ai relationship is a unique surface and independent of whether it is acquired through imbibition or drainage scanning curves. In these pore network studies and micro model experiments the determined interfacial area comprises only of the capillary (meniscus) interfacial area (domain 1) excluding the film area (domain 2). The inclusion of film area leads to linear dependence of interfacial area on saturation which prevents the generation of a unique Pc–S–Ai surface as was the case observed in Chen et al. [175] who determined the interfacial area with interfacial tracer tests which give information about the two interfacial domains (1 and 2) jointly. Joekar-Niasar et al. [176] developed an unstructured pore network model that resembles the topology of the micromodel used by Cheng et al. [173] to simulate their experiments. In addition to good agreement between the Pc–S–Ai surface attained from pore network model and from experiments, the authors also showed that the obtained surfaces were independent of whether they are obtained from drainage or imbibition curves as has been demonstrated earlier by Chen et al. [174]. Porter et al. [177] confirmed the hypothesis put forward by Hassanizadeh and Gray [169] using the Lattice Boltzmann (LB) method for multiphase flow (Shan and Chen type) within a glass bead domain that simulates the experiments conducted by Culigan et al. [111] (Fig. 5). The relatively minor discrepancy Fig. 5. Pc–S–Ai relationship obtained from LB model; from [177]. B. Agaoglu et al. / Advances in Water Resources 79 (2015) 162–194 between the interfacial areas determined via the LB method and Culligan et al.’s experiment was attributed to the unphysical spreading of the wetting phase during imbibition leading to increased film flow within the LB method. Porter et al. [178] demonstrated the unique surface of Pc–S–Ai with micro-tomography and also showed that this surface can be estimated from the thermodynamics method based on the area under the Pc–S curve developed by Grant and Gerhard [165]. 3.4.3. Entrapment Entrapment strongly influences the spatial distribution of phases in the porous media and the interfacial area which in turn have significant impact on interphase mass transfer. Entrapment of a fluid phase in porous media is directly related to pore scale movement of the intruding phase in the media which is directly related to the wettability of the media. In drainage the intrusion of non-wetting fluid is determined by an invasion process and piston-type motion prevails [179]. In imbibition the wetting phase moves also through the crevices, and snap-off mechanisms occur together with piston-type motion [180]. In this case the flow is less predictable as both mechanisms also rely on the surrounding topology of the phases. At the meso-scale the Pc–S curve in a porous domain provides information on the entrapped portion of the receding phase against the applied pressure on the advancing phase (i.e., the injection pressure of the advancing phase) while the pressure of receding phase is kept constant. The residual saturation represents the lowest saturation of the receding phase that can be attained. For NAPLs in water-wet media Anderson et al. [115] reported that the residual saturation will increase as the median grain size decreases. Using MRI Zhang et al. [181] showed that finer media will lead to greater residual saturation given the same wettability and grain size variance. A physical explanation for the occurrence of higher residual saturation is the increase in capillary force with decrease in grain size, and subsequently, pore throat, leading to entrapment of more NAPL blobs. Schnaar and Brusseau [37] showed that as the uniformity index (Ui = d60/d10) increases, the 175 residual saturation of NAPL (TCE) decreases for the same median grain size and wettability. Different correlations between residual saturation and porous medium properties (i.e., median grain size, wettability, pore size distribution) have been proposed by different authors. For example, the residual saturations measured by Powers et al. [9] do not have the same relationship with medium properties as Bradford et al. [89] or Schnaar and Brusseau [37] reported. One reason for these variations is that the meso-scale experiments present the impact of system parameters on entrapment jointly. Micro model experiments and pore network models provide a better picture of this process. Imbibition starts with piston type movement but in cases where the topology does not enable piston type movement snapoff would occur due to decrease in capillary pressure [180]. The two movement mechanisms have been shown to lead to different types of displacement patterns (i.e., flat-frontal; self-affine growth; bond percolation; compact cluster growth, etc.) [182]. Capillary number and contact angle have been demonstrated to have influence on these flow regimes [183]. It was also demonstrated that decreasing the aspect ratio suppresses snap-off hence, residual saturation and very low aspect ratio with high porosity may even lead to zero snap-off and residual saturation [176,182,183]. In drainage the only flow mechanism is piston type motion and therefore, pore size variance is responsible for entrapment of receding phase [179,184]. Slight occurrence of snap-off at the front of non-wetting phase is reported, yet it does not have influence as it gets connected to non-wetting phase with ongoing invasion at the front [4]. In drainage the viscosity ratio of the fluids has significant impact on entrapment of phases at different pressure levels. Lenordmand et al. [179] demonstrated that higher viscosity of the invading fluid relative to the receding fluid induces less entrapment via pore network modeling. Fig. 6 shows the different types of entrapment with respect to the capillary number and the viscosity ratio proposed by the pore network model of [179]. It was demonstrated that for unfavorable viscosity ratios (i.e., lower viscosity of the invading fluid) pore size variance has pronounced effect on the flow regime [184]. One implication of this Fig. 6. Emergence of fingering as a function capillary number and viscosity ratio based on Lenordmand et al. [179]; from [185]. 176 B. Agaoglu et al. / Advances in Water Resources 79 (2015) 162–194 phenomenon is that pore scale capillary forces dominate the development of the fingers in air sparging (Fig. 6). DNAPL entrapment in water-wet media resembles residual non-wetting phase occurrence during imbibition as water imbibes into a DNAPL region from above during downward leakage of DNAPL. On the other hand DNAPL entrapment in NAPL-wet media resembles residual wetting phase occurrence during drainage as water drains into a DNAPL region from above during downward leakage of DNAPL. For both of the cases the intrusion velocity of water into a DNAPL region impacts amount of entrapped NAPL and several authors proposed several relationships departing from the results of the pore network simulations [179,182,183]. In addition to velocity, the other parameters discussed in this section (such as contact angle, viscosities, pore size distributions, aspect ratio) determine the interfacial area at the pore scale which influence the interphase mass transfer at the meso-scale. This area is a subject of further research. 3.4.4. Wettability Wettability describes the tendency of a fluid to spread over a solid surface in the presence of another fluid. The configuration of fluids on the solid surface is affected by aqueous chemistry, solid phase mineralogy, organic matter distribution, surface roughness and charge, and contaminant aging [119]. Generally aquifer solids are accepted as water wet however non-water wetting or partially NAPL-wetting media can also be found [186]. Moreover, the wettability of a porous media can be altered through interaction with molecules of one of the phases. For example, in petroleum engineering it was reported that the adsorption of high molecular weight polar compounds existing in organic phase may lead to alteration of wettability [187]. Powers et al. [186] demonstrated that neat NAPLs such as toluene, TCE and isooctane do not alter the wettability of a quartz surface whereas complex petroleum and coal derived NAPLs can. The authors reported that creosote and coal tar makes a distinct alteration of the quartz surface causing it to become fully oil wet whereas crude oil and its derivations (e.g., gasoline, fuel oil, etc.) lead to a partially oil wet surface. The alteration of surface also has impact on the Pc–S curve whereby flushing of the NAPL contaminated area by water does not culminate into imbibition but drainage. However, during a DNAPL spill event the aquifer is already water saturated and the direct contact of the wettability-altering fluid with grain surface is limited [188]. The relative permeability of the water has also been shown to increase when the soil becomes more NAPL wetting [189,190]. Partially NAPL-wet media, which is likely to occur in cases where NAPLs have been in contact with the media longer than months [188], further complicates the formation of phases. Bradford et al. [89] confirmed this phenomenon on fractional wettable media which is a mixture of water-wet and NAPL-wet media. The authors showed that in water-wet media residual NAPL exists as single or multi-pore blobs. However, in fractional wettable media a NAPL film exist around the NAPL wet soil particles in addition to blobs. It was also reported that in fine porous media, which contains a small amount of NAPL-wet media (10%), NAPL residual saturation is lower compared to fully water wet media. This is attributed to the formation of NAPL films around NAPL wet grains which enable the movement of some of the NAPL globules through those NAPL wet pores. It was also found that further increase in the fraction of NAPL-wet media (>50%) leads to increase of residual saturation which was attributed to formation of NAPL film on many soil particles. A similar trend exists in coarse porous media which contains a small amount of NAPL-wet media (10%) such that NAPL residual saturation is lower compared to fully water-wet media. However when the fraction of the NAPL wet media increases (>50%) in coarse medium, the residual saturation does not increase further as it does in finer media. The authors also conducted dissolution experiments which demonstrated that in all types of porous media considered, increase in the NAPL wet fraction of the media leads to higher interphase mass transfer rate which was attributed to NAPL films around soil particles having a much larger interfacial area compared to non-wetting blobs retained in the pores. The authors also showed that within high fraction of NAPL-wet media (>75%), grain size is not influential on the interphase mass transfer rate, whereas for low fraction NAPL-wet media (10%), decreasing the median grain size leads to lower mass transfer rates. In addition Seyedabbasi et al. [119] showed that the dissolution fingering reduces when the medium becomes more NAPL-wet. 3.4.5. Mean grain diameter For ‘‘INtoMW’’, Mayer and Miller [191] showed that the median grain size and median NAPL blob size at residual NAPL saturation are positively correlated. Schnaar and Brusseau [192] reported that the ratio of the two values is between 1 and 2 for different types of porous media. Zhang et al. [181] also showed that the number of blobs in finer media is greater than the coarser one. Intuitively, greater mass transfer coefficients may be expected to occur in finer media (smaller median grain size) due to the existence of greater residual saturation and larger number of blobs, leading to more contact area between the flowing aqueous phase and entrapped NAPL. However, the mass transfer coefficients were found to be proportional to the median grain size in different studies such as [9,32,110,181]. One explanation of this apparent phenomenon is presented in Zhang et al. [181] which relates to the geometry of the entrapped NAPL. It was reported that the blobs occurring in coarser media are mostly spherical singlet blobs whereas in finer media they are distributed among singlet, doublet, triplet and multipore NAPL blobs. These multipore NAPL blobs induce stagnant water formations adjacent to NAPL–water interface, and as a result the flowing water tends to follow preferential flow paths without contacting the interface. In coarser media, on the other hand, the occurrence of singlets does not create stagnant water formations and the aqueous phase flows in a more distributed manner leading to more contact with the NAPL. The aforementioned flow field condition was also demonstrated using pore scale modeling (LatticeBoltzmann) by Knutson et al. [193] and in etched network experiments by Chomsurin and Werth [194]. The implication of this phenomenon on field scale is the retention of residual NAPL formations in less accessible media leading to further tailing of effluent concentrations. In air sparging the mean grain size influences the flow paths of the mobile air phase in water saturated medium. The radius of influence of the air sparging cone was shown to be influenced by the mean grain size. Reddy and Adams [195] showed that the radius of influence is greater for finer sand. On the other hand Peterson et al. [196] showed that sand with Average Modal Grain diameter (AMG, defined as the average grain size of sediment retained in the ASTM sieve with the largest volume yield) smaller than 1.3 mm leads to discrete air channels whereas AMG greater than 1.84 mm leads to a pervasive occurrence of air channels distributed as a symmetrical cone around the injection point. Between the two values a mixture of both flow types occur. Braida and Ong [197] and Rogers and Ong [198] showed that in sandy material the recovery rate decreases with decreasing mean grain diameter for the case where recovery is through volatilization of the solutes through the water–air interface (i.e., not directly from NAPL–air interface). The decrease is attributed to increased distance between the air–water interface and the NAPL–water interface owing to the greater tortuosity of the finer sand. However, their experimental set up consisted of a tank with custom made wide air channels, meaning that the results may not directly apply to air sparging. B. Agaoglu et al. / Advances in Water Resources 79 (2015) 162–194 3.4.6. Pore size distribution For ‘‘INtoMW’’, Mayer and Miller [191] showed that the pore size distribution of a porous media influences the formation of residual NAPL blobs whose size distribution can be defined in terms of a Van-Genuchten type equation. The authors did not quantify the relation between blobs and pore size distributions because all the tested porous media had overlapping pore size distributions. However, it was shown that the variance of the blob size distribution is greater than the variance of the pore size distribution. Schnaar and Brusseau [37,192] investigated this relation for residual NAPL (TCE) blobs using synchrotron X-ray microtomography. The authors showed that when the variance of the pore size distribution (defined through a uniformity index) increases, the variance of the NAPL blob size distribution increases. The number of blobs created in the non-uniform media is also found to be higher than the uniform one. One possible explanation for that can be found in Brooks et al. [132] who stated that when the grain size distribution becomes wider, the larger pores created by big grain particles are filled with small grain particles leading to the formation of smaller pore bodies. Under such conditions, the pore space will be governed by the small grain particles as observed in the X-ray images presented by [192]. Schnaar and Brusseau [37] also showed that the number of blobs with relatively higher volumes is greater in uniform media compared to non-uniform media. As discussed under the section ‘‘median grain size’’, smaller pore bodies might lead to higher number of NAPL blobs and multi-pore blobs which eventually induce NAPL bypassing and stagnant zones [181,193]. Consequently, high uniformity coefficient media may lead to smaller pore bodies and, therefore, a decrease in the mass transfer coefficients. In accordance with this argument Chomsurin and Werth [194] demonstrated on an etched network that less uniform media leads to a smaller mass transfer coefficient at the same Peclet number. Russo et al. [199] also showed that non-uniform media (hayhook soil Ui = 16, d50 = 0.26 mm) leads to 10 times smaller mass transfer coefficient than uniform media (Accusand 45/50, Ui = 1, d50 = 0.35). However in both studies the more non-uniform media is also the one with smaller median grain size which might be one other possible reason for the difference in mass transfer of these two media. Schnaar and Brusseau [37] on the other hand reported greater mass transfer coefficients for non-uniform sand in comparison to uniform sand having the same median grain size. However, greater flushing velocity was used in their dissolution experiment with the more non-uniform sand which may have also contributed to the greater mass transfer coefficient with the more non-uniform porous media. The impact of non-uniformness of the soil on the interphase mass transfer was also investigated by Powers et al. [9] who conducted column experiments using 6 different sand types whose uniformity index varied between 1 and 3. From their dissolution experiments the authors developed a Sherwood number formulation function of the uniformity index to a power of 0.41 implying that the more non-uniform media would lead to higher mass transfer coefficients. Even though the dissolution results comply with findings of Schnaar and Brusseau [37] it must be noted that the mean grain diameter is not the same among the sands used in Powers et al. [9] which makes a direct comparison not possible. In general it is reasonable to presume that non-uniform media will have smaller pore bodies due to filling of big pores with smaller grain particles in comparison with uniform media. However this might not always be the case. The development of a correlation between grain size distribution (i.e., uniformity coefficient) and pore size distribution for natural porous media would contribute substantially to the prediction of phase occurrences. Pore size distribution is also known to increase fingering when the viscosity ratio of the advancing phase to the receding phase is 177 very low [184]. Hence, the flow paths that occur during air drainage into water saturated media (e.g., during air sparging) are significantly influenced by the pore size distribution. In case of water imbibition into a NAPL pool in NAPL-wet media, pore size distribution has minor influence due to much smaller viscosity ratio between the NAPL and water. 3.4.7. NAPL mobilization NAPL mobilization in the subsurface is expected to occur when the exerted hydraulic pressure and gravitational force exceed the capillary pressure at the widest pore throat [200]. Observation of the main drainage and imbibition capillary pressure–saturation curve of a two-phase system indicates that mobilization stops when the residual saturation is reached. Mobilization is at maximum for a saturation value between the residual saturations of the two phases depending on the multiphase system [180]. DNAPL mobilization is mostly encountered either during the initial leakage period when gravitational forces induce fast downward movement or at a slower rate after the leakage and pooling over an inclined deep impermeable barrier. LNAPL mobilization on the other hand occurs mostly due to inclined water table level, exerted hydraulic gradient or changes in the water table level [200–202]. In an MRI study Johns and Gladden [157] showed that the velocity of the water during the leakage of the NAPL not only influences the entrapped NAPL amount, but also the average ganglion volume which decreases with higher water velocities. This was attributed to the higher pressured exerted on larger size ganglions due to their higher cross sectional area orthogonal to the flow direction. It must be noted that the diameters of the chosen glass beads ranged between 1 and 5 mm whereas the aqueous superficial flow velocity was very high compared to typical groundwater velocities. The influence of the velocity on average ganglion volume would likely be less significant at typical groundwater velocities. Recently Schnaar and Brusseau [37] reported that a formerly trapped NAPL blob can be mobilized due to induced blob diameter reduction induced by dissolution. The mutual mobilization of the two phases is likely to increase the mass transfer coefficient due to mechanical mixing occurring during flow of the two phases in porous media; however, no experimental observation of the alteration of the mass transfer coefficient has been reported in the literature. Mobilization becomes significant especially when remedial agents such as cosolvents and surfactants (i.e., ‘‘CHEM’’) are used. Both types of chemical agents can significantly decrease the interfacial tension between the NAPL and the aqueous phase which may lead to mobilization if the capillary force is overcome at the pore throat. Aydin et al. [47] demonstrated that mobilization is enhanced with increasing temperature due to reduced interfacial tension whereas Agaoglu et al. [54] showed that the use of intermediate levels of ethanol content in the flushing solution may lead to the development of preferential flow paths and the decrease of interphase mass transfer. 3.4.8. Desorption Zhu and Sykes [125] suggested that the dissolution of entrapped NAPL blobs (i.e., residual) can be defined through three stages. The first stage is defined by equilibrium or near equilibrium concentrations. The second stage is characterized by a rapid drop in effluent concentrations by three orders of magnitude. The third stage appears as tailing of the effluent concentrations which is associated with about 0.5% of the initial NAPL mass remaining in the aquifer. Imhoff et al. [203] stated that this last stage tailing may occur due to desorption or dissolution of isolated NAPL blobs. Gradual desorption of organic contaminants can significantly increase the time of remediation [204]. In particular, the slow diffusion of desorbed contaminants from stagnant regions of porous media (i.e., low permeability zones) and other mechanisms of 178 B. Agaoglu et al. / Advances in Water Resources 79 (2015) 162–194 kinetically slow desorption can cause long-term slow rates of solute flux which indicates that distinguishing between desorption and interphase mass transfer may be difficult in some cases. The processes and the modeling associated with desorption is, however beyond the scope of this work, and we refer the interested reader to [205–207]. 3.4.9. Density and viscosity Interphase mass transfer is influenced by the interfacial area which is also controlled by the stability of the interface during mobilization. Stability has been shown to be a function of density and viscosity differences [208]. Density difference in aquifer fluids leads to stratification of the multiphase system. Mackay et al. [209] reported that stratification forms when the density difference exceeds 1%. The impact of density is easily observed in various applications such as the use of remedial agents for groundwater remediation, air sparging activities and CO2 injection applications. The density of the remedial solution (i.e., ‘‘CHEM’’) strongly influences the formation of over- or under-riding of the flushing solution, and, consequently, the contact of fluids which has a crucial impact on the interphase mass transfer. Jawitz et al. [126] showed that the cosolvent solution which is comprised of 20% ethanol + 80% water, moves upwards during flushing due to density contrast against pure water (i.e., density difference of 2.7%) whereas the subsequent water flushing of the cosolvent solution leads to under-riding of water and trapping of cosolvent solution in capillary fringe. Furthermore, when the cosolvent solution comes in contact with NAPL, it drives the NAPL either upwards or downwards depending on its density further advancing stratification which decreases the field scale interphase mass transfer. In case of remediation of LNAPLs with densities higher than ethanol (e.g., toluene, benzene) the stratification and its negative effect on recovery can be overcome with the use of water ethanol mixtures which would decrease the density contrast that leads to stratification [210]. Grubb and Sitar [211,212] showed the occurrence of stratification when pure ethanol was flushed through DNAPL source zone. Another influence of density on interphase mass transfer is observed in high pressure CO2 injection into brine aquifer. The dissolution of the CO2 into brine aquifer at the CO2–brine interface at great amounts (due to its very high pressure) leads to density driven convection of CO2 enriched brine and collapsing of the boundary layer and, consequently, enhancing the interphase mass transfer [213,214]. 4. Mathematical modeling of interphase mass transfer A large number of mathematical models have been developed in the past several decades that focus on interphase mass transfer in porous media, particularly in relation to the dissolution of NAPLs into the surrounding flow. These models have generally included other fate and transport processes, such as advection and diffusion, because of the dependence of interphase mass transfer on these processes. The models have mostly been used as predictive tools for the estimation of contaminant concentrations in the vicinity of NAPL zones as well as in inverse mode to quantify the impact of different operational parameters and to back calculate interphase mass transfer parameters based on observed experimental data. In reviewing this body of work we have broadly classified the developed models as analytical or numerical. Analytical solutions are generally easier to apply and are in many instances used for screening purposes since they do not require detailed information about the subsurface conditions; however, their drawback is the oversimplification of the subsurface multiphase system. The scale that these models are developed for ranges from field scale applications with NAPL zones that can be on the order of 10 s of meters to meso-scale experiments with NAPL zones on the order of a few cm. Numerical models on the other hand generally provide more flexibility by relaxing some of the limiting assumptions incorporated in analytic models, yet they require more information about subsurface properties and the NAPL distribution and substantially larger effort to simulate the multiphase system. Numerical models that incorporate interphase mass transfer can be broadly divided into pore scale models and REV-based models. Pore scale models which focus on the interaction of interphase mass transfer and other multiphase flow processes at the pore scale mostly based on pore network models, are used to quantify interphase mass transfer at the meso-scale. Two other types of pore-scale models namely, Lattice Boltzmann and volume of fluid method have also been proposed for modeling interphase mass transfer at pore scale. REV-based models represent the interphase mass transfer at the continuum scale, and hence can be used to simulate interphase mass transfer at field scales. Three different types of REV-based models can be identified as: single phase models, multiphase models and an intermediate category, referred to in this review as semi-multiphase models which simulate the flow of a single phase (e.g., the aqueous phase) that is influenced by the existence of the other phase. 4.1. Analytical models The first models that focused on interphase mass transfer in porous media were developed to simulate NAPL dissolution (e.g., ‘‘INtoMW’’). These models consisted of analytical solutions of the partial differential equation (PDE) describing solute transport in porous media that is the advection–dispersion equation. The connection of the solute transport equation to interphase mass transfer is accomplished through a prescribed concentration equal to the solubility at the NAPL–water interface [34]. Hunt et al. [200] was one of the first papers that specifically considered dissolution from different NAPL configurations entrapped in porous media. For pool NAPLs oriented parallel to the flow direction, the steady-state concentration as a function of vertical distance from the pool is: 0 1 y C B C ðL; yÞ ¼ C ðsÞ erfc@ 1=2 A Dy L 2 Ux ð11Þ For a region of uniformly distributed spherical NAPL blobs, Hunt et al. [200] derived the following PDE to describe the transient state NAPL dissolution: @C ðx; t Þ kf ðds Þ 6V s ðx tÞ ¼ ½C s C ðx t Þ @x ds ðx t Þ U ð12Þ The above equation expresses the dissolved contaminant concentration gradient as a function of interfacial area, the difference between the concentration at solubility and initial concentration, and the mass transfer coefficient which is expressed as a function of the NAPL blob diameter. The authors reduced the above PDE to an ordinary differential equation (ODE) by solving it for a particular point in time knowing the volume of NAPL (Vs; assumed to be in shape of spherical blobs) and diameter of spherical blobs (ds). The above equations which are applicable at the meso-scale as they require continuous NAPL formation were also applied by the authors to a larger scale with NAPL length of 10 m. At such lengths the field scale permeability distribution and the non-continuous NAPL distribution, which are not accounted for in the above formulations, have significant impact on the interphase mass transfer. B. Agaoglu et al. / Advances in Water Resources 79 (2015) 162–194 Hunt et al. [200] also derived expressions for the transient mass transfer coefficient based on blob mass balance consideration. Both of the equations are based on the idealized assumption that NAPL is continuous in the space either as a pool or spherical blobs along the flow direction. While the above equations can be used to predict dissolved hydrocarbon concentrations near NAPL zones, a number of studies have used similar equations in an inverse mode to estimate the mass transfer coefficient of conducted meso-scale experiments [8,9,37]. These studies have also attempted to relate the estimated mass transfer coefficient to the experimental operational parameters such as flow velocity, saturation, and viscosity (Eqs. (1)–(4)). Similar approaches had been developed by Williamson et al. [28] and Wilson and Geankoplis [30] to determine the mass transfer coefficients of benzoic acid sphere dissolution in packed beds. Anderson et al. [115] used a Laplace approach to solve the governing PDE for both residual and pool NAPL dissolution at the field scale. This approach ignores the impact of field scale permeability (e.g., flow bypassing) on interphase mass transfer. Assuming the aqueous phase flowing through residual NAPL region is at solubility, the authors showed that the orthogonal position of a vertical NAPL finger with respect to groundwater flow leads to the occurrence of high concentrations at the down gradient of the finger due to advection of dissolved contaminants. On the other hand, a NAPL pool oriented parallel to groundwater flow induces lower concentrations because advection occurs parallel to the pool and dispersion is the only mechanism for the orthogonal spreading of the dissolved contaminant. Experimental observations that support these findings were discussed above in Section 2 of this review. In some instances water might penetrate the pool NAPL depending on the saturation as shown by Nambi and Powers [121]. Hatfield and Stauffer [215] developed a one dimensional analytical model for residual NAPL dissolution at the meso-scale that incorporates the mobile and immobile regions of water (due to clogging of a pore by NAPL) into the transport equation which also includes the sorption of solutes on solid surfaces. The derived partial differential equations for mobile and immobile regions were solved using a Laplace transform approach. Soerens et al. [216] also divided the porous media into two domains at the meso-scale namely: a NAPL region which represents flow paths that are in contact with NAPL and a clean region which represents flow paths that bypass the NAPL. Unlike the Hatfield model, Soerens et al. [216] assumed flow through the NAPL region. The total effluent concentration was defined as the sum of the concentration of the NAPL region multiplied by the fraction of flow through that region and the concentration of the clean region multiplied by the fraction of flow through that region. The effluent concentrations in both regions were determined through Eq. (4) whereby the interfacial area is lumped into mass transfer rate coefficient. The authors demonstrated the dependence of effluent concentrations on interphase mass transfer coefficients and other operational parameters (e.g., fraction of flow) by investigating three different scenarios through which (i) the clean region concentration was set to 0 when interphase mass transfer in the NAPL region is assumed to be rate limited, (ii) the clean region is set to rate limited interphase mass transfer when the NAPL region is assumed to be at equilibrium conditions or (iii) both regions are set as rate limited interphase mass transfer with NAPL region having a mass transfer coefficient 100 times greater than that of the clean region. Chrysikopoulos et al. [217,218] computed the dissolution of a DNAPL pool residing beneath the aquifer in 2-dimensional and 3dimensional spaces, respectively. Holman and Javandel [219] solved for the dissolution of an LNAPL pool residing on top of water table. The solutions presented in these studies were derived by solving the governing differential equation using a Laplace transformation approach. An analytical model referred to as multiple 179 analytical source superposition technique (MASST) which enables the solution for multiple NAPL sources was presented by Sale and McWhorter [220]. The model is based on the analytical solution derived by Hunt [221] for the solute transport equation incorporating one dimensional advection, three-dimensional dispersion and a source term which refers to instantaneous occurrence of solute due to interphase mass transfer. The source term was defined through a first order mass transfer term, M ¼ kðC s C a Þ. With this term the solution considers the influence of existing concentration (Ca) on source term. The authors expanded the methodology to multiple NAPL source zones. For that, a domain (here named subzone) is defined as a summation of subdivisions and the concentration at the centroid of each subdivision is found with respect to other subdivisions and their respective existing concentrations which are unknown. Writing this summation equation for each subdivision leads to N number of equations with N number of unknowns (Ca, the concentrations at the center of each subdivision) with the mass transfer coefficient of each subdivision defined by the user (Eq. (13)) [220]. Since this methodology does not incorporate the impact of field scale heterogeneity and the corresponding issues on interphase mass transfer, it is more suitable for meso-scale investigations. The authors applied it on both meso and field scale domain. N X C a xj ; yj ; zj ¼ K l ðC s C a ðxi ; yi ; zi ÞÞ F xj xi ; yj yi ; zj zi i¼1 ð13Þ where; i; j ¼ 1; 2; 3; . . . ; N ð14Þ This procedure enables MAAST to define the initial and boundary conditions in a more detailed manner including discontinuous NAPL formations a condition often encountered in the field. . Parker and Park [38] developed an analytical solution to estimate the effluent concentrations at the field scale (with a domain size of 10 m) based on the ODE defining the concentration change with respect to location (Eq. (4)). The authors assumed that the exponential solution of the governing ODE can be approximated by a linear function when the mass transfer coefficient is very small (k 1) which is often the case in the field: C out jeff L q C eq for jeff L q 1 ð15Þ ð16Þ The authors coupled this equation with the change of NAPL mass in the domain with time by correlating the mass transfer b2 M coefficient to change of mass with a power beta i:e:; MðtÞ0 : b1 M ðtÞ b2 C out j0 L q ¼ s K q C eq M0 ð17Þ The power beta, b2 ; is an empirical parameter (also named ‘‘mass depletion exponent’’) which is reported to be a function of field conditions and was obtained through fitting the analytical solution to measured data. Zhu and Sykes [125] presented a NAPL dissolution model with the effluent concentration is either a linear (L-model) or non-linear (N1-model; N2-model) function of the change of NAPL mass in the domain. The constant parameters are chosen by fitting the analytical solution to experimental results. Similar approaches were presented by Falta et al. [222]. Christ et al. [11] investigated the relation between the mass depletion exponent and aquifer 180 B. Agaoglu et al. / Advances in Water Resources 79 (2015) 162–194 contamination conditions in the field. The authors proposed that ganglia to pool ratio (GTP, the saturation of NAPL ganglia formations to saturation of pool NAPL formations in the domain) which has been investigated by Lemke et al. [223] is a better predictor of NAPL mass dissolution characteristics. Therefore, the authors suggested that the ‘‘mass depletion exponent’’ be estimated from the GTP: b ¼ 1:5 GTP0:26 ð18Þ Christ et al. [224] improved the estimation of the mass depletion exponent by considering whether the mass flux is from a NAPL pool or a ganglia region. The authors note that when residual NAPL exists together with pool NAPL, the time evolution of effluent concentrations shows two distinct trends: an initial high concentrations stage attributed to residual NAPL dissolution followed by a stage of lower concentrations resulting from pool NAPL dissolution. The authors stated that for pool formations the mass depletion exponent can be assumed to be 0.5 until the mass is almost depleted in this region whereas in regions with residual NAPL the exponent is given by the formulation proposed by Christ et al. [11] (Eq. (18)). In addition to the GTP, the flux contributions of each region, which are difficult to estimate in the field, are also input parameters in this approach. The ratio of the solute mass flux reduction to NAPL mass reduction has also been considered at the field scale as a possible parameter for estimating the duration of remediation activities [225]. Brusseau et al. [112] showed that pool NAPL formations lead to approximately 1 to 1 flux reduction against mass reduction which can be represented by the GTP approach [11] fairly well whereas uniformly distributed residual NAPL leads to the occurrence of high effluent concentrations until 80% of NAPL mass is depleted. Hence, the use of the GTP approach does not cover this behavior. Similar results were also reported by Difilippo et al. [226]. Both of the above studies also used a simple power function model to represent the mass flux reduction against NAPL mass reduction and indicated that the multi-step behavior which may appear in certain cases cannot be simulated with this approach. A recent decision tool for NAPL remediation, REMCHLOR [227], also incorporates this simple function. Stream tube analysis, which is based on a stochastic advective transport approach was also used to investigate the relationship between mass flux reduction and mass reduction in the field [228,229]. Jawitz et al. [228] showed that the increase in variance of the reactive travel time, s distribution between streamlines, which is a function of the hydraulic conductivity variance coupled with the variance of spatial distribution of NAPL leads to faster decrease of the effluent concentration: si ¼ ti Ri ¼ ti þ ti K f Si ð19Þ The required parameters of this method such as moments of reactive and non-reactive travel times can be deduced from tracer tests. It should be noted that the stream tube analysis is based on the conditions where the NAPL saturation is around or less than residual saturation such that the decrease in NAPL saturation (via dissolution) does not lead to the occurrence of new stream lines. Zhang et al. [108] used different analytical approaches including MASST and the stream tube approach to model four different 3D flow cell experiments each with different NAPL saturation and hydraulic conductivity fields. The authors showed that none of the models were able to capture all the change in effluent concentration with time. Chen et al. [230] evaluated the feasibility of the equilibrium stream tube approach for cosolvent/surfactant flushing applications (CHEM) over a 2D flow chamber where the average NAPL saturations are low (0.007–0.039). The occurrence of mobilization would prevent the use of stream tube analysis due to the change of NAPL configuration that was defined apriori through reactive tracer tests, however the authors observed no mobilization during flushing. Wood et al. [231] did not consider this phenomenon for predictions of cosolvent flushing with cosolvent content as high as 72% which likely initiates mobilization and changes the NAPL configuration defined earlier for the use of stream tube model. Braida and Ong [232] developed a simple model to predict vapor phase concentrations in air sparging applications based on the air–water interfacial area determined from the number and characteristics of air channels, the diffusive mass transfer zone (MTZ) in the aqueous phase around the gas channels and a userdefined rate limited volatilization parameter. The basis of the MTZ concept had been presented earlier by Chao et al. [233]. 4.2. Pore scale numerical models 4.2.1. Pore network models Flow in porous media can be represented as flow in a network of pore throats and pore bodies which is the founding aspect of pore network models [179]. In pore network models the flow is defined through Poisuielle’s law due to the laminar structure of flow in porous media if the domain chosen in the pore network model has a cylindrical cross section [179]. In case of rectangular shapes adapted versions of Poisuielle’s flow can be used by determining an equivalent pore radius [83,234]. Pore network models can be divided into two categories (i) quasi static models and (ii) dynamic models [235]. Dynamic models can simulate the transient behavior of the flow when the viscous forces are comparable to capillary forces [4,236]: N Vi i @Sai X þ Q aij ¼ 0 @t j¼1 ð20Þ However, the pore network models used to investigate interphase mass transfer processes are quasi static models which simulate the equilibrium states of drainage and imbibition processes when capillary forces dominate the flow which is the case in typical oil/water flow (e.g., for capillary number in the order of 107) [237]: Q ij ¼ Gij DPij ð21Þ X Q ij ¼ 0 j ¼ 1; 2; . . . N ð22Þ j Eq. (21) relates the flow through to the resistance and the pressure gradient while Eq. (22) is the conservation of mass equation. Another common feature of the pore network models investigating interphase mass transfer is that they assume the dissolving phase to be immobile, trapped in the pore throats/bodies. Therefore, displacement processes are not covered. Reviews of pore network models are presented by [4,235,238]. One of the first efforts to use pore network models for simulating NAPL dissolution was presented by Jia et al. [234]. The exact locations of NAPL blobs were obtained from visual observation of a micro model and incorporated into the pore network model. The time evolution of the NAPL constituents in the aqueous phase at the nodes which do not neighbor any NAPL–water interface was determined using the following equation with the flow velocity determined using Poisuielle’s law. Vs DC i X DC ij X Gij DPij ¼ þ DAij C ij Dt l l j j ð23Þ For nodes neighboring the NAPL–water interface it was assumed that the interface in the pore throat resembles either a flat or a cavity like shape. The equation describing the time 181 B. Agaoglu et al. / Advances in Water Resources 79 (2015) 162–194 Fig. 7. Schematic of NAPL formation in a chamber in the pore network model of Dillard et al. [90]. evolution of NAPL constituents in these neighboring nodes contains an extra first order mass transfer term with a mass transfer coefficient that is determined from the assumed interface shape (i.e., flat or cavity). The authors demonstrated that when the mass transfer process in the nodes neighboring the NAPL–water interface is well described, the results of the pore network model match quite well with experimental results. Dillard and Blunt [83] used a pore network model to simulate the NAPL dissolution experiments conducted by Powers et al. [9]. The authors used the invasion percolation method, developed by Wilkinson and Willemsen [239], to generate the fluid distributions. The pore size distribution of the network model was created by conditioning the simulated Pc–S curve to the experimental data of Powers et al. [9] by adjusting the size of pore throats. In the pore network model solute transport occurs through advection and diffusion with possible corner flow of the aqueous phase when the NAPL is trapped in a pore body or pore throat: F ij ¼ in Qw ij C t þ Dm Aw C out C in i j lt ! ð24Þ The interphase mass transfer was assumed to occur only during the corner flow of the aqueous phase (see Fig. 7) with a first order term (in this work named as Ed) is determined from the analytical approach presented by Lake [7] which is based on the solution of the advective–diffusive transport equation for flat fluid–fluid surface conditions. For a NAPL filled chamber i, solute concentration leaving the chamber through corner flow is given with respect to mass transfer coefficient: in C out ¼ C in i i þ Ed C s C i ð25Þ Unlike the work of Jia et al. [234], Dillard and Blunt [83] excluded NAPL dissolution into the aqueous phase through diffusion from the NAPL–water interface which resides at the pore body or throats facing NAPL-free pore body/throat. The authors demonstrated that the dependence of effluent concentrations to the Peclet number agrees well with experimental data for a constant NAPL volume assumption. In another pore network model the authors generated the time evolution data of effluent concentrations as the NAPL volume decreased with dissolution. However, slight discrepancy was observed between simulated and experimental data which was attributed to the complex nature of interphase mass transfer. A salient property of the difference between simulated and experimental effluent concentrations is the tailing which occurred in simulation but was not observed in the experimental data. This may be attributed to the fact that the size reduction of NAPL blobs due to dissolution instigated the mobilization of these blobs in the column experiment whereas in pore network model the NAPL blobs were set to be trapped permanently during dissolution. The column experiments conducted by Powers et al. [9] were also simulated by Zhou et al. [100] using a physics-based model which represents the porous media as bundle of parallel pores. The number of pores in a pore bundle was not found to be important as the geometry itself enables a realistic representation of the corner diffusion, the pore diffusion (i.e., the dissolution from the NAPL interface orthogonal to flow direction) and mixing of stream lines. The model was suitable to simulate column experiments which represent one-dimensional flow with relatively small cross-sectional area (e.g., Powers et al. [9] with a column diameter of 5.5 cm) because at each incremental distance from the inlet the flow can be assumed to be mixed as can be represented with the connection of pore bundles. This assumption however may not be valid for larger scale flows when flow diversion or low transverse dispersivity conditions are present, and might limit the use of this approach at larger scales. The authors also modeled interphase mass transfer in NAPL wet media to compare the results with the column experiments conducted by Parker et al. [240]. In this case the NAPL was assumed to be residing adjacent to the pore walls enabling aqueous flow through the center of the pores. Corner diffusion was assumed to be the only pattern for interphase mass transfer (however in this case from corner to center). The fraction of corner flow, b which is necessary for this model, is difficult to estimate in practice since it is a function of many properties which may preclude the use of this model in a predictive way. By matching model results to the experimental results, the authors were able to develop a relation between mass transfer coefficient and Peclet number. Dillard et al. [90] used the same pore network model presented by [83] to investigate the mass transfer rate coefficient, K l of an REV-based numerical model with respect to NAPL saturation and Peclet number. The authors also simulated DNAPL dissolution using an REV-based numerical model that incorporated the mass transfer formulation developed from their pore network model. For comparison the authors used the same REV-based numerical model with local equilibrium approach and with the mass transfer formulation developed by Powers et al. [9]. It was reported that the results were similar for all three simulations which was attributed to low velocities used in the simulations (about 1 m/d). The pore scale interphase mass transfer coefficient developed from the analytical solution presented in Lake [7], which was employed in the pore network models of [83,90,100], was also investigated by Sahloul et al. [82]. The authors solved the transport equation numerically in a corner with adequate boundary conditions to determine the relationship between mass transfer coefficient and flow conditions in that corner and showed that the analytical solution presented in Lake [7] underestimates the interphase mass transfer. It was reported that the discrepancy was caused by the chosen boundary conditions and the non-uniform velocity distribution assumed to exist in their corner model. Held and Celia [39] used the pore network model of Reeves and Celia [170] to generate NAPL distributions in a network whose pore structure was defined based on the properties of the porous medium used in the column experiments of Imhoff et al. [32] with incorporation of the reported pore size distributions. The generated fluid distributions were used to create the aqueous flow field in the network which was then used to simulate transport of dissolved constituents. No diffusive transport of solutes was incorporated in the model; only the diffusive flux appearing in the interphase mass transfer term was accounted for. Advective transport was modeled by tracking the dissolved mass explicitly in time (i.e., no matrix is solved due to discretization): V a;i X DC a;i X Q ij C a;j ¼ J j Ana;j þ t j2N j2M i i ð26Þ 182 B. Agaoglu et al. / Advances in Water Resources 79 (2015) 162–194 where; X i J i Ana;i ¼ X i Ana;i Dm S C a C a;i li ð27Þ Full depletion of NAPL in a pore body due to dissolution was followed by generation of the new aqueous velocity field. The authors excluded corner diffusion and the only pattern of interphase mass transfer was defined through mass flux from the NAPL into adjacent NAPL-free pore body via diffusion along the stagnant aqueous phase residing in the pore throat between. The authors reported excellent agreement with results reported by Inhoff et al. [32]. Yiotis et al. [2] focused on the drying of liquid phase in a porous domain which is closed from all sides except one where a fracture allowed gas phase to be purged through. The authors excluded the film formation in their pore network model and the mass flux from liquid to the gas phase at a pore containing a gas–liquid interface was determined with the stagnant gas phase assumption in that pore (i.e., a function of the molecular diffusion of liquid molecules in the gas phase and stagnant boundary layer length). In case of a liquid filled pore, that is adjacent to a fracture which has flowing gas phase, the mass flux was expressed in terms of the velocity of the gas phase in the fracture. The authors pointed to the influence of trapped wetting phase islands on interphase mass transfer. However, the exclusion of film formation, which is shown by Sahloul et al. [82] to be a dominant process in case of non-wetting phase intrusion, may have influenced the results significantly. Zhao and Ioannidis [88] investigated the same phenomena for a NAPLwet domain with the wetting NAPL phase trapped in the middle of the domain with the non-wetting water phase flowing through the surrounding pores. The authors represented interphase mass transfer from the NAPL filled pore/throat to the water filled pore/ throat with stagnant layer assumption and from the NAPL film residing at the corner to the water residing at the center of the pore/throat with a first order mass transfer coefficient which was determined with the formulation developed by Sahloul et al. [82]. The authors also incorporated the NAPL dissolution behavior in the NAPL wet domain that was investigated by [82] via micro model experiments. The NAPL filled chambers feeding the NAPL films eventually ceased as these films dissolve into the aqueous phase. The disconnection of films and the chambers was represented by assigning a critical value of capillary pressure at which NAPL connection is broken. This critical value was defined in the model as part of the input. The authors showed that the disconnection of NAPL filled chambers from anterior NAPL films lead to a sharp decrease in interphase mass transfer rates. This was attributed to the vanishing of NAPL film interfacial areas once the disconnected NAPL film depletes. Zhao and Ioannidis [241] expanded this model by determining the critical value of capillary pressure in terms of the pore radius: Pcrit ¼ Pmax xr=Rt c ð28Þ Zhao and Ioannidis [242] investigated the corner interphase mass transfer coefficient in a 3D corner geometry whereby the fluid–fluid interface (i.e., one side of the corner) is assumed to be at solubility and compared the results with the results presented by [82,83,100]. The authors considered both slip and no-slip conditions and found that representing the corner diffusion (Ed) via the 2D analytical solution of Lake [7] leads to erroneous predictions because in reality the corner has a 3D geometry which results in distinct velocity distributions that are not accounted in the 2D approach. The relation between Ed and dimensionless residence time of wetting phase in a 3D corner developed by Zhao and Ioannidis [242] was further used in a pore network model developed by Zhao and Ioannidis [93] (see Fig. 8) to investigate the CO2-gas exsolution processes (i.e., supersaturated CO2 concentration in aqueous phase dissolving into gaseous phase of CO2). In this model the commencement of gas bubbles was assigned to arbitrarily chosen sites. The authors demonstrated that the mass transfer characteristic is significantly influenced by the initiation of gas bubbles. 4.2.2. Other pore scale models Pore scale multiphase multi-component flow can also be simulated by a Shan-Chen type Lattice Boltzmann (SC-LB) model [243]. In this approach the forces acting on a node due to interaction with other nodes of other type components is added to the LB equation. This methodology accounts for the diffusion in the LB equation. The interaction force is determined as a function of the user-defined interaction strength parameter, G, and the separation distance between the corresponding nodes, e. The force acting at position x of component j due to interaction between component j and k (Gjk) for a multi-component system is given by: F ðjÞ ðx; t Þ ¼ WðjÞ ðx; tÞ k¼C X X Gjk ðjei jÞWðkÞ ðx þ ei Þei k¼1 Fig. 8. Schematic of pore network model with water flowing around a bubble. ð29Þ i Small values of G induce considerable amounts of diffusion which also means high miscibility whereas larger values lead to occurrence of purer phases. An alternative approach can be found in Swift et al. [244]. Most of the SC-LB implementations in porous media presented in this work have been conducted to investigate fluid–fluid interfaces or saturation–capillary pressure relations (e.g., [177,245,246]). Meakin and Tartakovsky [247] reported that the poor numerical stability of LB models in case of large density and viscosity contrasts prevents their use in geologically important systems (e.g., water–air flow in vadose zone) whereas incorporation of moderate density or viscosity ratios for the sake of stability leads to unrealistically high diffusion rates. Very recently new formulations of interaction force have also been presented that are suitable for use of higher density ratios [248,249]. It is important to note that LB models have very high computational costs in comparison to pore network models. However, the structure of the model enables easy parallelization [246,250]. B. Agaoglu et al. / Advances in Water Resources 79 (2015) 162–194 The Volume of Fluid method has also been used to simulate multiphase flow at pore scale. Raeini et al. [251,252] further developed the openFOAM two phase flow simulator to account for the sub-pore scale forces emerging in flow through porous medium. However solute transport and interphase mass transfer have not been yet considered in this approach. 4.3. REV-based numerical models 4.3.1. Single phase models Single phase numerical models are based on the solution of the solute transport equation in the aqueous phase while NAPL zone is assumed to be immobile. The effect of saturation on the relative permeability of the aqueous phase is also ignored. In single phase numerical models an interphase mass transfer term (e.g., first order mass transfer) is introduced into the transport equation for NAPL containing nodes. Powers et al. [35] simulated a typical NAPL (TCE) dissolution scenario (i.e., ‘‘INtoMW’’) through which the interphase mass transfer is defined as M ¼ kðC s CÞ. Depletion of the NAPL mass is accounted for by deducting the total mass of the dissolved contaminants from the remaining NAPL mass at each grid block and at each time step. Interphase mass transfer coefficients were obtained from different formulations developed by Pfannkuch [17], Wilson and Geankoplis [30], Williamson et al. [28] and Friedlander [25] that relate the mass transfer coefficient to system parameters such as flow rate, viscosity, density, porosity and saturation. The authors showed that predictions of the various formulations can vary significantly especially at high flow rates and low NAPL saturations. Borden [49] used a similar approach for multi-component NAPL (gasoline) dissolution using both equilibrium and a rate limited mass transfer approach with a mass transfer coefficient determined by matching experimental effluent concentrations. Powers et al. [95,110] developed ‘‘theta’’ and ‘‘sphere’’ models for the prediction of interphase mass transfer coefficients which were subsequently incorporated into single phase dissolution model. The sphere model incorporates the interfacial area between NAPL and aqueous phase while maintaining the Sherwood number as a function of Reynolds number only; that is, it ignores the impact of saturation, uniformity coefficient and other parameters. In theta model, the Sherwood number is a function of all parameters with the interfacial area lumped in it. The authors showed that the sphere model can simulate the effluent concentrations accurately provided the calibration parameter, F is known apriori. Baldwin and Gladden [102] solved the transport equation in the aqueous phase with NAPL saturation data that were determined using MRI data from column experiments. The authors also considered interphase mass transfer with a linear mass transfer model, a shrinking-core model and a pore-diffusion model for comparison. Significant deviations from the measured saturations were predicted by all three models. Johns and Gladden [156] improved the linear dissolution model of Baldwin and Gladden [102] by correlating the interfacial area to saturation. The NAPL saturation data used in the model were determined experimentally using MRI. The mass transfer coefficients were subsequently determined by fitting the numerical results to the experimental data. Kim and Chrysikopoulos [253] and Chrysikopoulos and Kim [101] simulated solute transport due to NAPL dissolution assuming that solute concentrations in NAPL-containing nodes of the domain were equal to the solubility value. The contributions of the sub REV parameters on the interphase mass transfer flux in the NAPL containing nodes such as flow field, saturations, interfacial area were not accounted for. Bradford et al. [104] modified the framework developed by Powers et al. [35] by including sorption. The authors 183 showed that the mass transfer formulations developed by [32,43,110] cannot simulate interphase mass transfer for different wettable media. Consequently, the authors developed a revised formulation for interphase mass transfer coefficient that considers the wettability of the media using data obtained from column experiments by Bradford et al. [89]. The developed model does not consider interfacial area separately but lumps it into the interphase mass transfer rate coefficient. Abriola et al. [75], Mayer et al. [77] and Schaerlaekens et al. [254] simulated the dissolution of residual NAPL using a surfactant solution (i.e., ‘‘CHEM’’) by incorporating the increased solubility value obtained from batch experiments in the mass transfer term. Formulations for the mass transfer coefficients were determined by fitting the numerical results to experimental data. Ji and Brusseau [255] presented a similar approach for enhanced dissolution by different chemical agents. 4.3.2. Semi-multiphase models Semi multiphase models consider the flow of only one phase that is influenced by the saturation of the other phase. For example, in case of ‘‘INtoMW’’ saturation of immobile phase has impact on the flow of the mobile aqueous phase by changing its relative permeability. Although the solute transport equation applied to the aqueous phase transport equation is mathematically similar to the single phase models, the incorporation of relative permeability in the model can fundamentally alter the flow field by causing the mobile aqueous phase bypass the NAPL which in turn can significantly impact transport simulations. Geller and Hunt [256] used this methodology to simulate multi-component NAPL dissolution. The authors used the Wyllie function to estimate the relative permeability of the aqueous phase based on the saturation and updated the interfacial area at each time step with respect to the depleted NAPL mass. The authors excluded dispersion because it is already accounted for in the mass transfer coefficient which is adopted from Hunt et al. [200]. This however also means that no dispersion is accounted at nodes containing no NAPL. Imhoff et al. [116,117] used a similar approach but included the dispersion term and used the mass transfer formulations of Powers et al. [110], Imhoff et al. [32] and Miller et al. [8] for comparison. The authors reported that the two formulations [32,110] predicted similar results but deviated from the experimentally measured results at later dissolution times. Powers et al. [257] used the flow and transport simulators, Visual MODFLOW and MT3D (Waterloo Hydrogeologic versions) in a sequential procedure to approximate NAPL dissolution. The relative permeability was defined through the Corey function [258] and used in MODFLOW to determine the velocity field which was subsequently used as an input for MT3D to simulate the solute transport. The change in saturation due to interphase mass transfer was represented as a decrease in saturation at each time step. The comparison of simulation and experimental results for high NAPL saturations demonstrated that the model is capable of capturing the main trend of the experimental data, although significant difference between the two was observed. This was attributed to the use of the local equilibrium assumption in the transport model. Frind et al. [124] followed a similar approach whereby the groundwater flow was simulated numerically to determine the hydraulic head distribution which was then used in the transport model to simulate contaminant transport. Unlike Powers et al. [257] the authors incorporated rate limited mass transfer in their model with the mass transfer coefficient determined with a new formulation developed by the authors. Saba and Illangaskare [105] also used MODFLOW and MT3D to simulate NAPL dissolution. The hydraulic conductivity of the medium in MODFLOW was defined equal to the relative permeability (estimated by Wyllie 184 B. Agaoglu et al. / Advances in Water Resources 79 (2015) 162–194 function) times the intrinsic permeability. The authors conducted a tank experiment and developed a new mass transfer formulation by regression analysis performed between effluent concentrations attained from experiment and the MT3D model. For verification the authors performed another tank experiment that had slightly more complex contamination area and simulated this experiment with MT3D that is revised to include the mass transfer formulations of Powers et al. [110] and Imhoff et al. [32] and the one developed from the tank experiment. The authors reported that the latter formulation predicted the observed contaminant plume better. The authors reported that the difference between the new mass transfer formulation and the previous ones is because the new formulation is developed from 2D dissolution experiment whereas other formulations [32,110] are developed from 1D column dissolution experiments. The significance of the dimensionality of the experiments was also evaluated by Brusseau et al. [259] who used a mass transfer formulation developed from 1D experiments to simulate 2D dissolution experiments. Contrary to Saba and Illangaskare [105], Brusseau et al. [259] stated that the simulations gave good agreement with experimental results even though the mass transfer coefficient is derived from 1D experiments. Imhoff et al. [118] used the numerical algorithm presented by Miller et al. [260] to investigate the dissolution fingering in two dimensional flow where the relative permeability was determined from the Wyllie expression whereas the rate limited mass transfer was determined from the formulation developed by Powers et al. [110]. The authors showed that increasing the variance of hydraulic conductivity leads to greater fingering, whereas increasing the transverse dispersivity lead to the suppression of fingering. Parker and Park [38] used a percolation model to generate a DNAPL saturation profile as the DNAPL leaks downward. The percolation model was followed by the MODFLOW and MT3D model to simulate the interphase mass transfer. The authors stated that modeling NAPL dissolution with this approach enables the use of very fine mesh (one million nodes for a 10 m3 domain) and time steps. Marble et al. [107] used the model presented by Brusseau et al. [259] with a new mass transfer formulation which is independent from the experimental conditions that the mass transfer coefficient is derived from. Seyedabbasi et al. [119] improved the model of [260] to account for the influence of wettability of the medium (i.e., different flow configurations) by incorporating: i) the Modified Mualem’s approach [190] which relates the relative permeability to the wettability of the media and ii) the interphase mass transfer formulation developed by Bradford et al. [104]. Mason and Kueper [261] simulated surfactant enhanced dissolution of pooled DNAPL (‘‘CHEM’’) using a revised semi multiphase approach whereby both saturation and interfacial tension control the relative permeability of the flushing solution. The DNAPL saturation distribution along the pool is calculated separately based on the capillary pressure at each location along the pool. The authors used two different approaches for the determination of interphase mass transfer. In the first approach the parameters that are required for the interphase mass transfer coefficients were determined by fitting the numerical results to experimental data. In the second approach the solubility value (Cs) appearing in the linear mass transfer term is defined as a function of time. The authors simulated the column experiments that they had conducted and reported that the second approach resulted in better predictions. Rathfelder et al. [262] modeled surfactant enhanced NAPL dissolution experiments using a modified version of the MISER model [263] whereby the use of surfactants do not reduce the interfacial tension to such level that mobilization is initiated. Apart from the rate limited mass transfer formulation which is obtained from Taylor et al. [264], the authors also considered the influence of density and viscosity alteration due to the use of surfactants on contaminant transport. The authors reported good agreement with the experimental results. 4.3.3. Multiphase models Prior to their use for environmental applications, multiphase models were mostly developed within the field of petroleum engineering. Here we specifically focus on how interphase mass transfer processes are implemented into multiphase models. For brevity we exclude the equations describing multiphase flow which have been presented in comprehensive reviews by Abriola [1] and Miller et al. [10]. The first multiphase model developed for environmental applications is the compositional model of Abriola and Pinder [265] which assumed equilibrium partitioning of the constituents between phases. Apart from Abriola and Pinder, early multiphase modeling efforts were based on immiscible flow of liquids whereby each fluid is assumed to be comprised of only one constituent (e.g., aqueous phase comprising only water) without interphase mass transfer [266–268]. Sleep and Sykes [269] presented a compositional model, solved with both the fully implicit and IMPES schemes, through which interphase mass transfer was determined according to local equilibrium assumption. White et al. [270] developed another local equilibrium compositional model STOMP using the multiphase formulations of the Richards equations [271]. Unger et al. [103] developed a three phase compositional model which allows for the interphase mass transfer to be defined as rate limited or at equilibrium. Adapted versions of mass transfer formulations developed by Mayer and Miller [99] and Guiger and Frind [272] were incorporated in the model. Application of these mass transfer formulations to 20 different NAPL contamination scenarios, each with a different hydraulic conductivity field but having the same variance, showed that the uncertainty of the time required for the dissolution of all NAPL mass is greatest for the formulation of [99] whereas the equilibrium partitioning assumption induced the lowest uncertainty. The authors also conducted a set of simulations where they scaled the relative permeability–saturation curve of each node with respect to the permeability of that node. It was demonstrated that the uncertainty of the time required for the dissolution of all the NAPL mass decreased with this scaling procedure which also induced more lateral spreading of the NAPL. No scaling was performed on the capillary pressure– saturation curve. It was also reported that decreasing the variance of the permeability field leads to a reduction of the uncertainty of the time required for the complete dissolution of the NAPL mass. Dekker and Abriola [273] modified an IMPES type model (VALOR) by scaling the capillary pressure–saturation curve and residual saturation of each node with the permeability of that node to investigate the impact of these parameters on DNAPL spill geometry. Residual saturations were determined as a log-linear function of permeability whereas the Pc–S curve was scaled by changing the a parameter of the Van-Genuchten formulation with respect to permeability. The authors demonstrated that increasing the vertical correlation length of the permeability increased the penetration depth and decreased the spread and maximum NAPL saturations. The horizontal correlation lengths on the other hand did not have a significant impact. Hinkelmann et al. [274,275] further developed the MUFTE-UG model [276] to include two phase (gas–liquid) three component (methane–water–air) conditions. This was achieved without an additional mass transfer term in the mass balance but by coupling the mole fraction of methane in the water phase to its gaseous mole fraction using Henry’s law which is essentially an equilibrium condition [277]. Falta [278] modified the T2VOC code [279] with a dual-domain approach. The author defined two distinct zones whereby the first zone contains the NAPL and the aqueous phase which is assumed to be at equilibrium. The other zone contains only the aqueous B. Agaoglu et al. / Advances in Water Resources 79 (2015) 162–194 phase and acquires the dissolved solute concentration from the surface of the first zone. The mass flux from this surface is determined with the integral of the analytical solution presented by Hunt et al. [200] which encompasses transverse dispersivity and the length of the NAPL containing zone. The NAPL pool is assumed to be immobile. The main advantage of this modeling approach is the possibility of using coarser grid cells and the reduction of computation times. Grant and Gerhard [33] used the DNAPL3D code (which solves the immiscible multiphase flow equation for a water–NAPL system) and MT3D sequentially to simulate NAPL dissolution experiments by first simulating DNAPL migration over a time step followed by contaminant transport over the same time step. The authors used three different approaches to simulate interphase mass transfer: (i) a local equilibrium approach, (ii) a rate-limited mass transfer approach that explicitly accounted for the interfacial area, and (iii) a rate-limited mass transfer approach without interfacial area (e.g., a lumped mass transfer approach). In the case of approach (ii), the authors determined the interfacial area using the thermodynamics-based approach developed by Grant and Gerhard [165] which considers the area under the Pc–S curve. The mass transfer coefficient was determined from the Peclet number as proposed by Pfannkuch [17]. In the case of lumped rate-limited mass transfer, the mass transfer rate coefficient was determined using the formulation developed by Saba and Illangaskare [105]. The authors reported that the rate-limited mass transfer approach that explicitly accounted for the interfacial area produced the best agreement with the experimental data. It must be noted that this kind of modeling of interphase mass transfer necessitates the choice of an interphase mass transfer coefficient that is independent of the interfacial area. The authors have chosen the formulation of Pfannkuch [17] which is actually adopted from the work of Hoffmann [96]. In this work the dependence of interphase mass transfer coefficient to velocity is only determined for a specific NAPL–water interfacial area meaning that the relationship between interphase mass transfer coefficient and Peclet number is only valid for the particular interfacial area used in the experiment of [96]. Furthermore, in approach (ii) the mass transfer rate always increases with interfacial area. However, this assumption is not always true as the mass transfer depends on the geometry of the entrapped NAPL and the flow field around it. Hence, greater interfacial area in some instances might lead to smaller mass transfer (see section mean grain diameter). Moreover, this method accounts for the meniscus interfacial area only, ignoring the contribution of film interfacial area to interphase mass transfer. Especially in ‘‘IWtoMN’’ interfacial area has been shown to be the dominant contributor of interphase mass transfer (see [82]). Kokkinaki et al. [97] followed a similar methodology and incorporated the interfacial area that is calculated by Grant’s method [165] into the model developed by Sleep and Sykes [269] to model interphase mass transfer that explicitly accounts for interfacial area. The authors also compared different interphase mass transfer formulations to 3D NAPL dissolution experiment conducted by Zhang et al. [108]. Kokinnaki et al. [98] developed simplified thermodynamic model and power law model for non-hysteric case by comparing the Grant’s method of interfacial area calculation to saturation parameter that appears as a surrogate for interfacial area in simple Sherwood formulation. Basu et al. [280] used the UTCHEM code to simulate DNAPL leakage and to conduct tracer tests to be used in their stream tube model. The DNAPL was flushed with water to compare the interphase mass transfer pattern with the stream tube model predictions. UTCHEM enables the use of local equilibrium or a rate limited mass transfer approach with the mass transfer coefficient determined from the formulation developed by Miller et al. [8]. The authors stated that even at high velocities, the use of local 185 equilibrium assumption and the rate limited approach both yielded the same average effluent concentrations. This was attributed to the fact that NAPL architecture had a dominant influence on the non-equilibrium conditions. Basu et al. [281] employed the T2VOC model to investigate how the first and second moments of the flux distribution on a control plane down gradient of the source zone changes with NAPL mass reduction. The mean and standard deviation of the flux distribution was found to decrease with source mass depletion by dissolution. Maji and Sudicky [84] used the Compflow model to simulate DNAPL dissolution using the hydrogeological data set of a real site in Tübingen/Germany whose parameters such as permeability, effective porosity, grain size, mineralogy and sorption coefficients were determined with high resolution. The authors generated 20 realizations of the hydraulic conductivity field using a transition probability/Markov chain (TP/MC) method that uses the Tubingen data for conditioning. 10 different mass transfer formulations were incorporated in the model to determine the interphase mass transfer coefficient (including the formulations of Miller et al. [8], Powers et al. [110], Imhoff et al. [32], Saba and Illangasekare [105], Geller and Hunt [256], Nambi and Powers [121] and 4 other Sn-based models with differing parameters). The authors reported that the model of Saba and Illangasekare [105] and Nambi and Powers [121] deviated from others in terms of time required for the dissolution of all the NAPL mass whereas the other formulations culminated in dissolution times similar to the local equilibrium approach. The authors also showed that even after significant NAPL mass reduction (99.99%), the concentrations at a plane located at 10 m down-gradient of the source can exceed regulatory targets (i.e., maximum concentration level; MCL). It was also demonstrated that the mass flux/mass reduction relationships emanating from the different realizations may vary vastly which also depends on the mass transfer formulation used. O’Carroll and Sleep [282] simulated hot water flushing of NAPL with the compositional model developed by Sleep and Sykes [269] which assumed local equilibrium conditions. It was demonstrated that increasing the temperature has significant impact on NAPL mass recovery due to alteration of viscosity However, uncontrolled NAPL downward migration might pose a greater risk at elevated temperatures. Reitsma and Kueper [63,283] developed the first fully implicit two-phase three-component (cosolvent–water–NAPL) compositional model that simulates cosolvent flushing of PCE (CHEM). The model assumes the fluids to be non-ideal whereby diffusion of one type of molecule is calculated with Maxwell– Stefan Diffusivity as a function of the composition of the multicomponent fluid. The detailed derivation of the mass flux concept and its relation to the diffusion coefficient is presented in Taylor and Krishna [56]. The phase behavior was determined by Hand’s method, which states that equilibrium phase concentrations are straight lines on a log–log plot, because of its practicality. This approach had been also integrated into UTCHEM for surfactants [62]. Interphase mass transfer flux was defined through the twofilm model and the diffusion coefficient in the boundary layer was assumed to be a constant. The interphase mass transfer rate was set equal to the molar flux times the specific interfacial area: Ii ¼ N i a ð30Þ The interfacial area was determined with respect to the saturation of that specific node with use of an equation that ignores the effect of hysteresis, developed from a pore network model [284]. The alteration of interfacial tension with respect to cosolvent content was determined based on the expression proposed by Li and Fu [285] which is also used to scale the relative permeability. The authors fit their simulation to experimental data to determine the film thickness appearing in film model. As mentioned in Section 2, the film model with a hypothetical film thickness does 186 B. Agaoglu et al. / Advances in Water Resources 79 (2015) 162–194 not apply directly to interphase mass transfer in porous media and, therefore, the determined film thickness is only valid for that particular experiment. Roeder and Falta [286] modified the UTCHEM code for cosolvent-enhanced NAPL recovery and demonstrated that fine meshed heterogeneous 3D domains can simulate unstable flow conditions which may occur during cosolvent flushing of DNAPLs. Fluid viscosities were determined with method of Grunberg whereas interfacial tension alteration as a function of cosolvent content was determined with the Li and Fu method [285]. Liang and Falta [287] used UTCHEM to simulate NAPL dissolution with very high cosolvent content in the flushing solution (%95 ethanol + %5 water) and demonstrated that the rate limited interphase mass transfer behavior vanishes at such high cosolvent contents. Kaye et al. [288] used the same model to investigate the occurrence of overand under-riding of the NAPL mass for different flushing solutions. Agaoglu et al. [54] simulated cosolvent flushing of multicomponent NAPL using UTCHEM and demonstrated that bypassing of the NAPL mass at intermediate cosolvent contents, which was observed in laboratory experiments conducted by the authors, could not be fully captured by the model. The development of preferential paths and NAPL bypassing was attributed to the mobilization of the NAPL mass residing within the larger pores due to the cosolvent-induced low interfacial tensions even for relatively homogeneous systems. Use of finer 3D meshes did not improve the simulation results because of the difficulty of accurately representing pore scale heterogeneity in the model. Lingieni and Dhir [289] simulated air sparging of NAPL contaminated soil (‘‘IWtoMN’’) by first solving the immiscible multiphase flow equations to determine air velocity which was then used to simulate the transport of the volatilized contaminant. Interphase mass transfer from NAPL to air modeled with both equilibrium and a rate limited approach, with the mass transfer coefficient related to the Reynolds and Schmidt numbers. Comparison of model simulations to column experimental results showed that better agreement was obtained with the rate limited approach. Lundegrad and Andersen [145] investigated the radius of influence and air flow path distribution during air sparging with a multiphase flow model initially developed for petroleum reservoir simulations (TETRAD). The development of radius of influence was also investigated by McCray and Falta [146] and Hein et al. [290] using the T2VOC code. These studies focusing on radius of influence were based on air–water system excluding NAPL. Mclure and Sleep [291] employed the compositional model developed by Sleep and Sykes [269] to investigate the bio-venting of soils to promote contaminant (toluene in this case) biodegradation. The parameters that have significant impact on this type of bio-remedial activity, namely: oxygen distribution in aqueous phase from the injected air, contaminant partitioning into the aqueous phase from the NAPL, contaminant evaporation and extraction through the injected air, were all modeled using a local equilibrium approach. Unger et al. [292] used the Compflow model to simulate air sparging of soil contaminated with a DNAPL. It was reported that at early stages recovery is high due to direct contact of NAPL with air flow-paths and direct vaporization. However, after this initial stage the recovery decreases and is controlled by the NAPL dissolution rate into the aqueous phase because contaminants must first dissolve into the aqueous phase and subsequently must travel to the air flow-paths for eventual recovery. Thomson and Johnson [137] showed that multiphase flow models may not be capable of accurately predicting mass recovery from air sparging activities due to the emergence of discrete capillaries leading to the formation of channels which are not accounted for in continuum-based models unless a very fine discretized grid is used which in practice is not possible. The presence of these channels is attributed to the formation of viscous fingering when air flows in a water-saturated medium and its sensitivity to pore scale heterogeneity. Yoon et al. [293] incorporated the rate limited mass transfer formulation of Miller et al. [8] to the STOMP model to investigate the impact of non-equilibrium conditions on soil vapor extraction applications. Jang and Aral [139] also used a rate limited mass transfer approach in a multiphase flow model to investigate the impact of different design parameters of air sparging such as injection rate, under various site conditions. Xu et al. [294] considered the simulation of CO2 into a deep sandstone aquifer and its immobilization through carbonate precipitation. In this flow configuration the CO2 is the mobile nonwetting phase that dissolves into the groundwater (i.e., ‘‘MNtoIW’’). The model was based on the computer code TOUGHREACT which is a part of the TOUGH2 suite of programs [295]. TOUGHREACT first solves the multiphase flow equation to determine phase velocities and uses this information in the reactive transport equation. The TOUGHREACT code assumes equilibrium partitioning of constituents between phases. To date nonequilibrium interphase mass transfer has not been used in CO2 simulations. CO2 injection in the subsurface involves high pressures and concentrations gradients resulting in significant convective mixing. Thiebeau and Dutin [296] used the computer code Eclipse to investigate the impact of convective mixing on a large field scale study (i.e., on the orders of kilometers). The computational grid used in this study was 10 m in the vertical direction. The authors stated that convective mixing (which is a result of density difference between two locations) and its influence on interphase mass transfer cannot be captured when such relatively large computational grids are employed. The authors attempted to overcome this problem by upscaling the diffusion coefficient of CO2 in aqueous phase to represent convective mixing. Enouy et al. [94] used the Compflow model to simulate column experiments involving injection of aqueous phase supersaturated with CO2 (i.e., ‘‘MWtoIN’’). Injection of aqueous phase supersaturated with CO2 has been recently investigated as a potential surrogate to air sparging activities to extract NAPL within exsolution gas or to promote biodegradation if oxygen is chosen as injected solute. The potential of this procedure was stated to lie in the more distributed gas formation in comparison to air sparging at least in the gas exsolution region because the flow paths in air sparging are highly dependent on the pore scale capillary pressure heterogeneity. The authors modeled the partitioning of CO2 from the supersaturated aqueous solution to the gas phase which emerges as nucleation, based on the mass transfer coefficient adopted from the formulations of Nambi and Powers [36] and Unger et al. [103]. Recently Niessner and Hassanizadeh [3,297] modeled the multiphase flow approach proposed by Hassanizadeh and Gray [169,298] that also incorporates the mass balance of specific interfacial area. The time evolution of the interfacial area was set equal to the interfacial area production rate, EWN, minus the advective transport of interfacial area. This latter term is set equal to the interfacial area permeability times the gradient of interfacial area which is related to the gradient of Gibbs free energy. The parameter EWN was determined from: Ewn ðSw ; Pc Þ ¼ ewn ðSw ; Pc Þ @Sw @t ð31Þ ewn is the strength of change in specific interfacial area determined by linear interpolation between its two limits which is determined based on the assumption that those limits occur at the beginning and end of drainage. The authors simulated a case where the non-wetting phase was injected into the porous media saturated with the wetting phase B. Agaoglu et al. / Advances in Water Resources 79 (2015) 162–194 (e.g., drainage). Interphase mass transfer between the wetting and non-wetting phases was also incorporated in the model. The interphase mass flux was expressed in terms of the mole fractions with the interphase mass transfer coefficient derived based on the film layer theory. However, this approach is not suitable for an REVbased model as film thickness (named diffusion length in this work) is not a known parameter for a computational node of REV-based model. It is a function of the various properties of the porous media as has been described in this Sections 2 and 3. As a result the authors defined this parameter as a fit parameter. The constitutive function between capillary pressure–saturation- interfacial area incorporated in the model was the one developed from the pore network model of Joeker-Niasar et al. [172]. The dissolution from the immobile wetting phase into the non-wetting phase (i.e., ‘‘IWtoMN’’) occurs mainly from film interfacial area as reported by Sahloul et al. [82] and Zhao and Ioannidis [88,241]. However this pore network model did not account for the film interfacial area. Recently multiscale modeling of multiphase flow in porous media was presented by Jackson et al. [299] which formulated the mass, momentum and energy balances over the volumes, interfacial areas and common curves with constraints arising from entropy inequality, the second law of thermodynamics and other equilibrium thermodynamics relations [300]. The multiphase system considered in this approach is immiscible [301]. 5. Summary and conclusions This review paper examines the factors influencing interphase mass transfer in porous media and the mathematical models developed to analyze this process. The large volume of studies reviewed in this manuscript attests to the significance of interphase mass transfer in porous media in various applications and the challenges encountered in characterizing this process. Because of the significance of scales on the apparent interphase mass transfer, we have framed this review according to the length scales of the reviewed studies, namely: the pore-scale, meso, and field scales. The characteristics of the multiphase flow regime in porous media is another important factor that influences interphase mass transfer. Therefore, to gain further insight, published studies were also discussed in terms of five flow configurations defined based on the mobility and wettability of the phases (Table 1). The review also surveys various mathematical models that have been developed for the analysis and simulation of interphase mass transfer in porous media, particularly those relating to the dissolution of NAPLs into the surrounding flow. Published studies have shown that the interphase mass transfer is influenced by a number of key parameters in a complex manner. Identifying these relationships is essential for further improvement in our understanding of and ability to model mass transfer. Foremost among these parameters is the interfacial area separating the interphases within the porous media. Recent technological advances in non-invasive imaging are promising in identifying the pore space geometry and phase distributions within. Such accurate determination of the interfacial area can potentially be incorporated into REV-based models and other modeling efforts where the lack of precise control on the interfacial area and how it varies with time is cited as a key factor contributing to uncertainties in model prediction. However, reported data clearly indicates that the relationship between the interfacial area and mass transfer rate is quite complex. For example; in case of residual NAPL, the porous medium with smaller median grain size generally leads to greater specific interfacial area. However, the generation of multipore blobs within finer medium induces stagnant regions around the interfacial domains, bringing about smaller mass transfer rates 187 in comparison to coarser media where the specific interfacial area is smaller but with less stagnation regions. In a related issue no REV-based model to date has attempted to incorporate separately the film and meniscus area data into the mass transfer term. The main reason for that is the difficulty in determining the distribution of the flow rates of the mobile phase within the film and meniscus areas of discrete pore throats. One promising avenue for future research is pore network modeling which can be used to efficiently evaluate the impact of different trapped phase configurations on interphase mass transfer for different realizations of pore sizes. The pore space and the variance between pore sizes/throats have been shown to have a significant impact on interphase mass flux through their influence on different parameters and mechanisms such as interfacial areas, entrapment, and flow field. Some authors have used the uniformity coefficient to account for the impact of pore size distribution in phenomenological models. However using the uniformity coefficient as a surrogate for the pore size distribution might not always be correct, as in some situations big pores may be filled with smaller grains meaning that pore space is mostly determined by the smaller grains. Non-invasive techniques for the observation of porous media with different uniformity coefficients and pore space may shed further light on the correlation between uniformity coefficients and pore size distribution. Wettability of the porous medium is another important parameter influencing interphase mass transfer as it has significant implications on a number of parameters, most notably the flow paths and interfacial domains. The impact of wettability on the pore-scale phase distribution especially in partially oil-wet media is still lacking. Moreover, the potential alteration of wettability with the factors influencing this mechanism needs further investigation, especially in relation to NAPL remediation activities. Most models ignore the impact of wettability because of the lack of reliable information on its influence even though the significance of wettability on multiphase flow and interphase mass flux is acknowledged. Besides wettability, the impact of aging, ionic strength and temperature are mostly ignored in existing models even though they may have significant impact on mass transfer. Flow instability is another important parameter that can have major implication on interphase mass transfer. The formation of instabilities has direct influence on interphase mass flux through their influence on interfacial areas and flow paths. Modeling efforts have shown that REV-based multiphase flow models are generally unable to simulate discrete capillaries that arise in viscous fingering due to the variance in pore sizes. The incorporation of spatially variable parameters defined on threedimensional, fine computational grids can partially address this issue but at a large computational cost and up to a limited extent because discrete pore throats are much smaller than the smallest possible mesh. There is the suggestion in the literature that the mobilization of NAPL entrapped in porous media may also influence the mass transfer characteristics. Recently mobilization of the trapped blobs which are reduced in size due to dissolution has been depicted by non-invasive methods. Micro-models may be used to capture this process and to provide some data to quantify this mass transfer mechanism. In groundwater remediation applications miscibilityenhancing agents can be used at intermediate contents whereby the interfacial tension is sufficiently reduced to cause some of the entrapped NAPL, particularly entrapped in larger pores, to be mobilized and lead to development of preferential flow paths similar to dissolution fingering. Most efforts to predict enhanced NAPL recovery have relied on equilibrium assumptions or have used mass transfer expressions that are originally developed for water–NAPL interphases. More broadly, there is a need to develop 188 B. Agaoglu et al. / Advances in Water Resources 79 (2015) 162–194 specific studies that address interphase mass transfer during enhanced NAPL recovery In this review we referred to a large number of phenomenological models that have attempted to relate the interphase mass transfer to key operational parameters. These models are generally deduced from meso-scale based experiments with limited information on pore scale processes. Although these models provide insight into the factors influencing interphase mass transfer and are useful as a predictive tool, it is important to note that the predictions of these different models are in some instances inconsistent with each other even when the same input parameters are used (see [83,84,90,98,102–108]). This is mostly attributed to differences in the specific experimental conditions used in developing these models. Further detailed studies that can isolate the impacts of each parameter are needed to develop more generally applicable interphase mass transfer correlations. In particular, the combination of meso-scale studies along with accurate pore-scale characterization of the porous media and the multiphase system via non-invasive methods and pore scale models may be an optimal avenue for further progress in the effort towards the development of generalized, accurate and robust Sherwood formulations that can be incorporated in REV-based models. Pore network models are especially useful in providing detailed information of the porous medium and flow that can be incorporated into these expressions, while non-invasive characterization techniques can provide testing of the assumptions used to develop pore network model. Attempts to incorporate more physics-based processes into Sherwood formulations, which have been mostly derived from empirical data, is one promising avenue for future research. A salient feature of field scale applications is the complex spatial structure of the permeability field and the NAPL distribution and, due to limited data, the lack of complete knowledge of these properties. The challenges observed at the field have led in some instances to the use of REV based numerical models to investigate field scale NAPL dissolution behavior. In a recent study it has been shown that the use of different mass transfer formulations for a well characterized aquifer results in different NAPL depletion times [84] whereas others stated that the difference is insignificant [38]. In many field cases, NAPL contamination are several decades old and to a large extent, the accessible residual NAPLs have already been depleted, with the remaining NAPL existing mostly in the form of pools. In modeling of such conditions, the use of local equilibrium in a grid block of an REV-based model that contain NAPL may be viable as they are likely to contain more NAPL than the threshold saturation that leads to achievement of solubility. However, further work is needed to evaluate the impact of dissolution of the inaccessible NAPLs on overall field scale dissolution performance, once they reduce below residual levels. The development of more precise mass transfer formulations that can be used in REV-based numerical models specifically applied to field applications will allow for the more accurate evaluation and prediction of NAPL behavior in the field. In summary, it is evident that interphase mass transfer in porous media is a result of the complex, dynamic multiphase flow occurring within porous media. It is equally evident that significant progress has been achieved in recent decades to better understand this critical process. This review shows that interphase mass transfer is governed by both pore scale characteristics of the multiphase system, such as wettability and pore size distributions, as well as large scale flow attributes such as NAPL bypassing due to porous media heterogeneities at the field scale. Interphase mass transfer is the results of interplay of all these factors. Therefore, further work must set sights on the characterization of the pore scale system and how it relates to the interphase mass transfer. Equally important, research efforts must also attempt to develop robust means of representing the larger scale dynamic nature of multiphase flow field within the porous medium and its influence on the interfacial fluid contact and mass transfer. References [1] Abriola LM. Modeling multiphase migration of organic chemicals in groundwater systems – a review and assessment. Environ Health Perspect 1989;83:117–43. http://dx.doi.org/10.1029/WR021i001p00011. [2] Yiotis AG, Stubos AK, Boudouvis AG, Yortsos YC. A 2-D pore-network model of the drying of single-component liquids in porous media. Adv Water Resour 2001;24:439–60. http://dx.doi.org/10.1016/S0309-1708(00)00066-X. [3] Niessner J, Hassanizadeh SM. Modeling kinetic interphase mass transfer for twophase flow in porous media including fluid–fluid interfacial area. Transp Porous Media 2009;80:329–44. http://dx.doi.org/10.1007/s11242-009-9358-5. [4] Joekar-Niasar V, Hassanizadeh SM. Analysis of fundamentals of two-phase flow in porous media using dynamic pore-network models: a review. Crit Rev http://dx.doi.org/10.1080/ Environ Sci Technol 2012;42:1895–976. 10643389.2011.574101. [5] Mercer J, Cohen R. A review of immiscible fluids in the subsurface: properties models characterization and remediation. J Contam Hydrol 1990;6. http:// dx.doi.org/10.1016/0169-7722(90)90043-G. [6] Dwivedi PN, Upadhyay SN. Particle–fluid mass transfer in fixed and fluidized beds. Ind Eng Chem 1977;16:157–65. http://dx.doi.org/10.1021/ i260062a001. [7] Lake WL. Enhanced oil recovery. Englewood Cliffs, NJ: Prentice-Hall; 1989. [8] Miller CT, Poirier-McNeill MM, Mayer AS. Dissolution of trapped nonaqueous phase liquids: mass transfer characteristics. Water Resour Res 1990;26:2783–96. http://dx.doi.org/10.1029/WR026i011p02783. [9] Powers SE, Abriola LM, Weber WJ. An experimental investigation of nonaqueous phase liquid dissolution in saturated subsurface systems: steady state mass transfer rates. Water Resour Res 1992;28:2691–705. http://dx.doi.org/10.1029/92WR00984. [10] Miller CT, Christakos G, Imhoff PT, McBride JF, Pedit JA, Trangenstein JA. Multiphase flow and transport modeling in heterogeneous porous media: challenges and approaches. Adv Water Resour 1998;21:77–120. http:// dx.doi.org/10.1016/S0309-1708(96)00036-X. [11] Christ JA, Ramsburg CA, Pennell KD, Abriola LM. Estimating mass discharge from dense nonaqueous phase liquid source zones using upscaled mass transfer coefficients: an evaluation using multiphase numerical simulations. Water Resour Res 2006;42. http://dx.doi.org/10.1029/2006WR004886. [12] Anslyn EV. Physical organic chemistry. Univ Sci Books; 2006. [13] Schwarzenbach RP, Gschwend PM, Imboden DM. Environmental organic chemistry. John Wiley Sons; 2005. [14] Israelachvilli J. Intermolecular and surface forces. San Diego: Academic Press Inc.; 1992. [15] Zumdahl SS, DeCoste DJ. Chemical principles. Cengage Learning; 2011. [16] Weber WJ, McGinley PM, Katz LE. Sorption phenomena in subsurface systems: concepts, models and effects on contaminant fate and transport. Water Res 1991;25:499–528. http://dx.doi.org/10.1016/00431354(91)90125-A. [17] Pfannkuch H. Determination of the contaminant source strength from mass exchange processes at the petroleum–ground–water interface in shallow aquifer systems. In: Water Well Assoc Am Pet Inst Houston, TX; 1984. [18] Baehr HD, Stephan K. Heat and mass transfer. Berlin, Heidelberg: Springer; 2006. http://dx.doi.org/10.1007/3-540-29527-5. [19] Whitman W. Gas transfer to and across an air–water interface. Chem Metall Eng 1923;29:146–8. http://dx.doi.org/10.1111/j.2153-3490.1977.tb00746.x. [20] Deacon EL. Gas transfer to and across an airt–water interface. Tellus 1977;29:363–84. http://dx.doi.org/10.1111/j.2153-3490.1977.tb00746.x. [21] Cussler EL. Diffusion: mass transfer in fluid systems. Cambridge University Press; 2009. [22] Clark MM. Transport modeling for environmental engineers and scientists. Wiley Interscience; 1996. [23] Frössling N. Über die Verdunstung fallender Tropfen. Gerlands Beitrage Zur Geophys; 1938. [24] Ranz W, Marshall W. Evaporation from drops. Chem Eng Prog 1952;48(3):141–6. [25] Friedlander SK. Mass and heat transfer to single spheres and cylinders at low Reynolds numbers. AIChE J 1957;3:43–8. http://dx.doi.org/10.1002/ aic.690030109. [26] Dobry R, Finn RK. Mass transfer to a cylinder at low Reynolds numbers. Ind Eng Chem 1956:540–3. http://dx.doi.org/10.1021/ie51400a044. [27] Kusik CL, Happel J. Boundary layer mass transport with heterogeneous catalysis. Ind Eng Chem Fundam 1962;1:163–72. http://dx.doi.org/10.1021/ i160003a002. [28] Williamson JE, Bazaire KE, Geankoplis CJ. Liquid-phase mass transfer at low Reynolds numbers. Ind Eng Chem Fundam 1963;249:3. [29] Pfeffer R, Happel J. An analytical study of heat and mass transfer in multiparticle systems at low Reynolds numbers. AIChE J 1964;10:605–11. http://dx.doi.org/10.1002/aic.690100507. [30] Wilson EJ, Geankoplis CJ. Liquid mass transfer at very low Reynolds numbers in packed beds. Ind Eng Chem Fundam 1966;5:9–14. http://dx.doi.org/ 10.1021/i160017a002. B. Agaoglu et al. / Advances in Water Resources 79 (2015) 162–194 [31] Van der Waarden M, Bridie ALAM, Groenewoud WM. Transport of mineral oil components to groundwater – I. Water Res 1971;5:213–26. http://dx.doi.org/ 10.1016/0043-1354(71)90054-6. [32] Imhoff PT, Jaffé PR, Pinder GF. An experimental study of complete dissolution of a nonaqueous phase liquid in saturated porous media. Water Resour Res 1994;30:307–20. http://dx.doi.org/10.1029/93WR02675. [33] Grant GP, Gerhard JI. Simulating the dissolution of a complex dense nonaqueous phase liquid source zone: 2. Experimental validation of an interfacial area-based mass transfer model. Water Resour Res 2007;43. http://dx.doi.org/10.1029/2007WR006039. [34] Genuchten M.Van., Alves W. Analytical solutions of the one-dimensional convective–dispersive solute transport equation, No. 157268, United States Department of Agriculture, Economic Research Service, 1982. [35] Powers SE, Loureiro CO, Abriola LM, Weber WJ. Theoretical study of the significance of nonequilibrium dissolution of nonaqueous phase liquids in subsurface systems. Water Resour Res 1991;27:463–77. http://dx.doi.org/ 10.1029/91WR00074. [36] Nambi IM, Powers SE. Mass transfer correlations for nonaqueous phase liquid dissolution from regions with high initial saturations. Water Resour Res 2003;39. http://dx.doi.org/10.1029/2001WR000667. [37] Schnaar G, Brusseau ML. Characterizing pore-scale dissolution of organic immiscible liquid in natural porous media using synchrotron X-ray microtomography. Environ Sci Technol 2006;40:6622–9. http://dx.doi.org/ 10.1021/es0508370. [38] Parker JC, Park E. Modeling field-scale dense nonaqueous phase liquid dissolution kinetics in heterogeneous aquifers. Water Resour Res 2004;40. http://dx.doi.org/10.1029/2003WR002807. [39] Held RJ, Celia MA. Pore-scale modeling and upscaling of nonaqueous phase liquid mass transfer. Water Resour Res 2001;37:539–49. http://dx.doi.org/ 10.1029/2000WR900274. [40] Setchenow M. Action de l’acide carbonique sur les solutions des sels a acides forts. Ann Chim Phys 1892;25. [41] Lesage S, Brown S. Observation of the dissolution of NAPL mixtures. J Contam Hydrol 1994;15:57–71. http://dx.doi.org/10.1016/0169-7722(94)90010-8. [42] Almeida MB, Alvarez AM, Miguel EMD, Hoyo ESD. Setchenow coefficients for naphthols by distribution method. Can J Chem 1983;61:244–8. http:// dx.doi.org/10.1139/v83-043. [43] Imhoff PT, Frizzell A, Miller CT. Evaluation of thermal effects on the dissolution of a nonaqueous phase liquid in porous media. Environ Sci Technol 1997;31:1615–22. http://dx.doi.org/10.1021/es960292x. [44] Horvath AL. Halogenated hydrocarbons: solubility-miscibility with water. CRC Press; 1982. [45] Tyn MT, Calus WF. Diffusion coefficients in dilute binary liquid mixtures. J Chem Eng Data 1975;20:106–9. http://dx.doi.org/10.1021/je60064a006. [46] Javaheri M, Abedi J, Hassanzadeh H. Linear stability analysis of doublediffusive convection in porous media, with application to geological storage of CO2. Transp Porous Media 2009;84:441–56. http://dx.doi.org/10.1007/ s11242-009-9513-z. [47] Aydin GA, Agaoglu B, Kocasoy G, Copty NK. Effect of temperature on cosolvent flooding for the enhanced solubilization and mobilization of NAPLs in porous media. J Hazard Mater 2011;186:636–44. http://dx.doi.org/10.1016/ j.jhazmat.2010.11.046. [48] Conrad SH, Wilson JL, Mason WR, Peplinski WJ. Visualization of residual organic liquid trapped in aquifers. Water Resour Res 1992;28:467–78. http:// dx.doi.org/10.1029/91WR02054. [49] Borden RC, Piwoni MD. Hydrocarbon dissolution and transport: a comparison of equilibrium and kinetic models. J Contam Hydrol 1992;10:309–23. http:// dx.doi.org/10.1016/0169-7722(92)90013-5. [50] Mukherji S, Peters CA, Weber WJ. Mass transfer of polynuclear aromatic hydrocarbons from complex DNAPL mixtures. Environ Sci Technol 1997;31:416–23. http://dx.doi.org/10.1021/es960227n. [51] Luthy RG, Ramaswami A, Ghoshai S, Merkel W. Interfacial films in coal tar nonaqueous-phase liquid–water systems. Environ Sci Technol 1993:2914–8. http://dx.doi.org/10.1021/es00049a035. [52] Denekas MO, Carlson FT, Moore JW, Dodd CG. Materials adsorbed at crude petroleum–water interfaces isolation and analysis of normal paraffins of high molecular weight. Ind Eng Chem 1951;43:1165–9. http://dx.doi.org/10.1021/ ie50497a048. [53] Oostrom M, Dane JH, Wietsma TW. A review of multidimensional, multifluid intermediate-scale experiments. Vadose Zone J 2006;5:570. http://dx.doi.org/ 10.2136/vzj2005.0125. [54] Agaoglu B, Scheytt T, Copty NK. Laboratory-scale experiments and numerical modeling of cosolvent flushing of multi-component NAPLs in saturated porous media. J Contam Hydrol 2012;140:80–94. http://dx.doi.org/10.1016/ j.jconhyd.2012.07.005. [55] Jawitz JW, Dai D, Rao PSC, Annable MD, Rhue RD. Rate-limited solubilization of multicomponent nonaqueous-phase liquids by flushing with cosolvents and surfactants: modeling data from laboratory and field experiments. Environ Sci Technol 2003;37:1983–91. http://dx.doi.org/10.1021/es0256921. [56] Taylor R, Krishna R. Multicomponent mass transfer. New York: Wiley; 1993. [57] Kooijman H, Taylor R. Estimation of diffusion coefficients in multicomponent liquid systems. Ind Eng Chem Res 1991;30:1217–22. http://dx.doi.org/ 10.1021/ie00054a023. [58] Imhoff PT, Gleyzer SN, Mcbride JF, Vancho LA, Okuda I, Miller CT. Cosolventenhanced remediation of residual dense nonaqueous phase liquids: [59] [60] [61] [62] [63] [64] [65] [66] [67] [68] [69] [70] [71] [72] [73] [74] [75] [76] [77] [78] [79] [80] [81] [82] [83] [84] [85] 189 experimental investigation. Environ Sci Technol 1995;29:1966–76. http:// dx.doi.org/10.1021/es00008a014. Abrams DS, Prausnitz JM. Statistical thermodynamics of liquid mixtures: a new expression for the excess Gibbs energy of partly or completely miscible systems. AIChE J 1975;21:116–28. http://dx.doi.org/10.1002/aic.690210115. Lee KY, Peters CA. UNIFAC modeling of cosolvent phase partitioning in nonaqueous phase liquid–water systems. J Environ Eng 2004;130:478–83. http://dx.doi.org/10.1061/(ASCE)0733-9372(2004)130:4(478). Hand D. Dineric distribution. J Phys Chem 1930;34(9):1961–2000. http:// dx.doi.org/10.1021/j150315a009. Delshad M. UTCHEM Version 6. 1 technical documentation. Austin, Texas Cent Pet Geosys-Tems Eng Univ Texas Austin; 1997. Reitsma S, Kueper BH. Non-equilibrium alcohol flooding model for immiscible phase remediation: 1. Equation development. Adv Water Resour 1998;21:649–62. http://dx.doi.org/10.1016/S0309-1708(97)00024-9. Chang YC, Moulton RW. Quternary liquid systems with two immiscible liquid pairs. Ind Eng Chem 1953:2350–61. http://dx.doi.org/10.1021/ ie50526a054. Yalkowsky SH. Solubility and solubilization in aqueous media. Oxford Univ Press; 1999. Brandes D, Farley KJ. Importance of phase behavior on the removal of residual DNAPLs from porous media by alcohol flooding. Water Environ Res 1993;65:869–78. Hayden NJ, Van der Hoven E. Alcohol flushing for enhanced removal of coal tar from contaminated soils. Water Environ Res 1996;68:1165–71. Ramakrishnan V, Ogram AV, Lindner AS. Impacts of co-solvent flushing on microbial populations capable of degrading trichloroethylene. Environ Health Perspect 2005;113:55–61. http://dx.doi.org/10.1289/ehp.6937. Abdul AS, Gibson TL, Devi NR. Selection of surfactants for the removal of petroleum products from shallow sandy aquifers. Ground Water 1990;28:920–6. http://dx.doi.org/10.1111/j.1745-6584.1990.tb01728.x. Fortin J, Jury WA, Anderson MA. Enhanced removal of trapped non-aqueous phase liquids from saturated soil using surfactant solutions. J Contam Hydrol 1997;24:247–67. http://dx.doi.org/10.1016/S0169-7722(96)00013-7. Carroll BJ. The kinetics of solubilization of nonpolar oils by nonionic surfactant solutions. J Colloid Interface Sci 1981;79:126–35. http:// dx.doi.org/10.1016/0021-9797(81)90055-2. Grimberg SJ, Nagel J, Aitken MD. Kinetics of phenanthrene dissolution into water in the presence of nonionic surfactants. Environ Sci Technol 1995;29:1480–7. http://dx.doi.org/10.1021/es00006a008. Zhong L, Mayer AS, Pope GA. The effects of surfactant formulation on nonequilibrium NAPL solubilization. J Contam Hydrol 2003;60:55–75. http:// dx.doi.org/10.1016/S0169-7722(02)00063-3. St-Pierre C, Martel R, Gabriel U, Lefebvre R, Robert T, Hawari J. TCE recovery mechanisms using micellar and alcohol solutions: phase diagrams and sand column experiments. J Contam Hydrol 2004;71:155–92. http://dx.doi.org/ 10.1016/j.jconhyd.2003.09.010. Abriola LM, Dekker TJ, Pennell KD. Surfactant-enhanced solubilization of residual dodecane in soil columns. 2. Mathematical modeling. Environ Sci Technol 1993;27:2332–40. http://dx.doi.org/10.1021/es00048a006. Johnson JC, Sun S, Jaffe PR. Surfactant enhanced perchloroethylene dissolution in porous media: the effect on mass transfer rate coefficients. Environ Sci Technol 1999;33:1286–92. http://dx.doi.org/10.1021/es980908d. Mayer AS, Zhong LR, Pope GA. Measurement of mass-transfer rates for surfactant-enhanced solubilization of nonaqueous phase liquids. Environ Sci Technol 1999;33:2965–72. http://dx.doi.org/10.1021/es9813515. Zimmerman JB, Kibbey TCG, Cowell MA, Hayes KF. Partitioning of ethoxylated nonionic surfactants into nonaqueous-phase organic liquids: influence on solubilization behavior. Environ Sci Technol 1999;33:169–76. http:// dx.doi.org/10.1021/es9802910. Sharmin R, Ioannidis MA, Legge RL. Effect of nonionic surfactant partitioning on the dissolution kinetics of residual perchloroethylene in a model porous medium. J Contam Hydrol 2006;82:145–64. http://dx.doi.org/10.1016/ j.jconhyd.2005.10.001. Lee DH, Cody RD, Hoyle EL. Laboratory evaluation of the use of surfactants for ground water remediation and the potential for recycling them. Ground Water Monit Remediat 2001;21(1):49–57. http://dx.doi.org/10.1111/j.17456592.2001.tb00630.x. Park S-K, Bielefeldt AR. Equilibrium partitioning of a non-ionic surfactant and pentachlorophenol between water and a non-aqueous phase liquid. Water Res 2003;37:3412–20. http://dx.doi.org/10.1016/S0043-1354(03)00237-9. Sahloul NA, Ioannidis MA, Chatzis I. Dissolution of residual non-aqueous phase liquids in porous media: pore-scale mechanisms and mass transfer rates. Adv Water Resour 2002;25:33–49. http://dx.doi.org/10.1016/S03091708(01)00025-2. Dillard LA, Blunt MJ. Development of a pore network simulation model to study nonaqueous phase liquid dissolution. Water Resour Res 2000;36:439–54. http://dx.doi.org/10.1029/1999WR900301. Maji R, Sudicky EA. Influence of mass transfer characteristics for DNAPL source depletion and contaminant flux in a highly characterized glaciofluvial aquifer. J Contam Hydrol 2008;102:105–19. http://dx.doi.org/10.1016/ j.jconhyd.2008.08.005. Kim H, Rao PSC, Annable MD. Gaseous tracer technique for estimating air– water interfacial areas and interface mobility. Soil Sci Soc. 1999:1554–60. http://dx.doi.org/10.2136/sssaj1999.6361554x. 190 B. Agaoglu et al. / Advances in Water Resources 79 (2015) 162–194 [86] Costanza-Robinson MS, Brusseau ML. Air–water interfacial areas in unsaturated soils: evaluation of interfacial domains. Water Resour Res 2002;38:131–7. http://dx.doi.org/10.1029/2001WR000738. [87] Costanza-Robinson MS, Harrold KH, Lieb-Lappen RM. X-ray microtomography determination of airwater interfacial area water saturation relationships in sandy porous media. Environ Sci Technol 2008:2949–56. http://dx.doi.org/10.1021/es072080d. [88] Zhao W, Ioannidis MA. Pore network simulation of the dissolution of a singlecomponent wetting nonaqueous phase liquid. Water Resour Res 2003;39. http://dx.doi.org/10.1029/2002WR001861. [89] Bradford SA, Vendlinski RA, Abriola LM. The entrapment and long-term dissolution of tetrachloroethylene in fractional wettability porous media. Water Resour Res 1999;35:2955–64. http://dx.doi.org/10.1029/ 1999WR900178. [90] Dillard LA, Essaid HI, Blunt MJ. A functional relation for field-scale nonaqueous phase liquid dissolution developed using a pore network model. J Contam Hydrol 2001;48:89–119. http://dx.doi.org/10.1016/S01697722(00)00171-6. [91] Leonenko Y, Keith DW. Reservoir engineering to accelerate the dissolution of CO2 stored in aquifers. Environ Sci Technol 2008;42:2742–7. http:// dx.doi.org/10.1021/es071578c. [92] Taku Ide S, Jessen K, Orr FM. Storage of CO2 in saline aquifers: effects of gravity, viscous, and capillary forces on amount and timing of trapping. Int J Greenhouse Gas Control 2007;1:481–91. http://dx.doi.org/10.1016/S17505836(07)00091-6. [93] Zhao W, Ioannidis MA. Gas exsolution and flow during supersaturated water injection in porous media: I. Pore network modeling. Adv Water Resour 2011;34:2–14. http://dx.doi.org/10.1016/j.advwatres.2010.09.010. [94] Enouy R, Li M, Ioannidis MA, Unger AJ. Gas exsolution and flow during supersaturated water injection in porous media: II. Column experiments and continuum modeling. Adv Water Resour 2011;34:15–25. http://dx.doi.org/ 10.1016/j.advwatres.2010.09.013. [95] Powers SE, Abriola LM, Dunkin JS, Weber WJ. Phenomenological models for transient NAPL–water mass-transfer processes. J Contam Hydrol 1994;16:1–33. http://dx.doi.org/10.1016/0169-7722(94)90070-1. [96] Hoffmann B. Über die ausbreitung gelöster kohlenwasserstoffe im grundwasserleiter. Inst F Wasserwirtschaft U Landwirtschaftl Wasserbau D Techn Univ; 1969. [97] Kokkinaki A, O’Carroll DM, Werth CJ, Sleep BE. Coupled simulation of DNAPL infiltration and dissolution in three-dimensional heterogeneous domains: process model validation. Water Resour Res 2013;49:7023–36. http:// dx.doi.org/10.1002/wrcr.20503. [98] Kokkinaki A, O’Carroll DM, Werth CJ, Sleep BE. An evaluation of Sherwood– Gilland models for NAPL dissolution and their relationship to soil properties. J http://dx.doi.org/10.1016/ Contam Hydrol 2013;155:87–98. j.jconhyd.2013.09.007. [99] Mayer AS, Miller CT. The influence of mass transfer characteristics and porous media heterogeneity on nonaqueous phase dissolution. Water Resour Res 1996;32:1551–67. http://dx.doi.org/10.1029/96WR00291. [100] Zhou D, Dillard LA, Blunt MJ. A physically based model of dissolution of nonaqueous phase liquids in the saturated zone. Transp Porous Media 2000:227–55. http://dx.doi.org/10.1023/A:1006693126316. [101] Chrysikopoulos CV, Kim TJ. Local mass transfer correlations for nonaqueous phase liquid pool dissolution in saturated porous media. Transp Porous Media 2000:167–87. [102] Baldwin CA, Gladden LF. NMR Imaging of nonaqueous-phase liquid dissolution in a porous medium. AIChE J 1996;42:1341–9. http://dx.doi.org/ 10.1002/aic.690420514. [103] Unger AJ, Forsyth P, Sudicky E. Influence of alternative dissolution models and subsurface heterogeneity on DNAPL disappearance times. J Contam Hydrol 1998;30:217–42. http://dx.doi.org/10.1016/S0169-7722(97)00041-7. [104] Bradford SA, Phelan TJ, Abriola LM. Dissolution of residual tetrachloroethylene in fractional wettability porous media: correlation development and application. J Contam Hydrol 2000;45:35–61. http:// dx.doi.org/10.1016/S0169-7722(00)00118-2. [105] Saba T, Illangasekare TH. Effect of groundwater flow dimensionality on mass transfer from entrapped nonaqueous phase liquid contaminants. Water Resour Res 2000;36:971–9. http://dx.doi.org/10.1029/1999WR900322. [106] Zhu J, Sykes J. The influence of NAPL dissolution characteristics on field-scale contaminant transport in subsurface. J Contam Hydrol 2000;41:133–54. http://dx.doi.org/10.1016/S0169-7722(99)00064-9. [107] Marble JC, DiFilippo EL, Zhang Z. Application of a lumped-process mathematical model to dissolution of non-uniformly distributed immiscible liquid in heterogeneous porous media. J Contam Hydrol 2008;100:1–10. http://dx.doi.org/10.1016/j.jconhyd.2008.04.003. [108] Zhang C, Yoon H, Werth CJ, Valocchi AJ, Basu NB, Jawitz JW. Evaluation of simplified mass transfer models to simulate the impacts of source zone architecture on nonaqueous phase liquid dissolution in heterogeneous porous. J Contam Hydrol 2008;102:49–60. http://dx.doi.org/10.1016/ j.jconhyd.2008.05.007. [109] Zilliox L, Muntzer P, Menanteau J. Probleme de l’echange entre un produit petrolier immobile et l’eau en mouvement dans un milieu poreux. Rev IFP 1973;28:185–200. [110] Powers SE, Abriola LM, Weber WJ. An experimental investigation of nonaqueous phase liquid dissolution in saturated subsurface systems: transient mass transfer rates. Water Resour Res 1994;30:321–32. [111] Culligan KA, Wildenschild D, Christensen BSB, Gray WG, Rivers ML. Porescale characteristics of multiphase flow in porous media: a comparison of air–water and oil–water experiments. Adv Water Resour 2006;29:227–38. http://dx.doi.org/10.1016/j.advwatres.2005.03.021. [112] Brusseau ML, DiFilippo EL, Marble JC, Oostrom M. Mass-removal and massflux-reduction behavior for idealized source zones with hydraulically poorlyaccessible immiscible liquid. Chemosphere 2008;71:1511–21. http:// dx.doi.org/10.1016/j.chemosphere.2007.11.064. [113] Voudrias EA, Yeh MF. Dissolution of a toluene pool under constant and variable hydraulic gradients with implications for aquifer remediation. Ground Water 1994, http://dx.doi.org/10.1111/j.1745-6584.1994. tb00645.x. [114] Whelan MP, Voudrias EA, Pearce A. DNAPL pool dissolution in saturated porous media: procedure development and preliminary results. J Contam Hydrol 1994;15:223–37. http://dx.doi.org/10.1016/0169-7722(94)90026-4. [115] Anderson M, RL J, Pankow J. Dissolution of dense chlorinated solvents into groundwater. 3. Modeling contaminant plumes from fingers and pools of solvent. Environ Sci Technol 1992, http://dx.doi.org/10.1111/j.1745-6584. 1992.tb01797.x. [116] Imhoff PT, Miller CT. Dissolution Fingering During the Solubilization of Nonaqueous Phase Liquids in Saturated Porous Media: 1. Model Predictions. Water Resour Res 1996;32:1919–28. http://dx.doi.org/10.1029/96WR00602. [117] Imhoff PT, Thyrum GP, Miller CT. Dissolution Fingering During the Solubilization of Nonaqueous Phase Liquids in Saturated Porous Media: 2. Experimental Observations. Water Resour Res 1996;32:1929–42. http:// dx.doi.org/10.1029/96WR00601. [118] Imhoff PT, Farthing MW, Gleyzer SN, Miller CT. Evolving interface between clean and nonaqueous phase liquid (NAPL)-contaminated regions in twodimensional porous media. Water Resour Res 2002;38. http://dx.doi.org/10. 1029/2001WR000290. [119] Seyedabbasi MA, Farthing MW, Imhoff PT, Miller CS. The influence of wettability on NAPL dissolution fingering. Adv Water Resour 2008;31:1687–96. http://dx.doi.org/10.1016/j.advwatres.2008.08.003. [120] Brusseau ML, Nelson NT, Oostrom M, Zhang Z, Johnson GR, Wietsma TW. Influence of heterogeneity and sampling method on aqueous concentrations associated with NAPL dissolution. Environ Sci Technol 2000:3657–64. http:// dx.doi.org/10.1021/es9909677. [121] Nambi IM, Powers SE. NAPL dissolution in heterogeneous systems: an experimental investigation in a simple heterogeneous system. J Contam Hydrol 2000;44:161–84. http://dx.doi.org/10.1016/S0169-7722(00)00095-4. [122] Zhang C, Werth CJ, Webb AG. Characterization of NAPL Source Zone Architecture and Dissolution Kinetics in Heterogeneous Porous Media Using Magnetic Resonance Imaging. Environ Sci Technol 2007;41:3672–8. http://dx.doi.org/10.1021/es061675q. [123] Page JWE, Soga K, Illangasekare T. The significance of heterogeneity on mass flux from DNAPL source zones: an experimental investigation. J Contam Hydrol 2007;94:215–34. http://dx.doi.org/10.1016/j.jconhyd.2007.06.004. [124] Frind EO, Molson JW, Schirmer M, Guiger N. Dissolution and mass transfer of multiple organics under field conditions: The Borden emplaced source. Water Resour Res 1999;35:683–94. http://dx.doi.org/10.1029/1998WR900064. [125] Zhu J, Sykes JF. Simple screening models of NAPL dissolution in the subsurface. J Contam Hydrol 2004;72:245–58. http://dx.doi.org/10.1016/ j.jconhyd.2003.11.002. [126] Jawitz JW, Annable MD, Rao PS. Miscible fluid displacement stability in unconfined porous media. J Contam Hydrol 1998;31:211–30. http:// dx.doi.org/10.1016/S0169-7722(97)00062-4. [127] Ramsburg CA, Pennell KD. Density-modified displacement for DNAPL source zone remediation: density conversion and recovery in heterogeneous aquifer cells. Environ Sci Technol 2002;36:3176–87. http://dx.doi.org/10.1021/ es011403h. [128] Lunn SRD, Kueper BH. Manipulation of density and viscosity for the optimization of DNAPL recovery by alcohol flooding. J Contam Hydrol 1999;38:427–45. http://dx.doi.org/10.1016/S0169-7722(99)00008-X. [129] Martel R, Gelinas PJ, Desnoyers JE. Aquifer washing by micellar solutions: 1. Optimization of alcohol-surfactant-solvent solutions. J Contam Hydrol 1998;29:319–46. http://dx.doi.org/10.1016/S0169-7722(97)00029-6. [130] Van Valkenburg ME, Annable MD. Mobilization and entry of DNAPL pools into finer sand media by cosolvents: two-dimensional chamber studies. J Contam Hydrol 2002;59:211–30. http://dx.doi.org/10.1016/S01697722(02)00058-X. [131] Zhong L, Mayer A, Glass RJ. Visualization of surfactant-enhanced nonaqueous phase liquid mobilization and solubilization in a two-dimensional micromodel. Water Resour Res 2001;37:523. http://dx.doi.org/10.1029/ 2000WR900300. [132] Brooks MC, Wise WR, Annable MD. Fundamental changes in in situ air sparging how patterns. Ground Water Monit Remediat 1999, http://dx.doi. org/10.1111/j.1745-6592.1999.tb00211.x. [133] Peteron JW, Murray KS, Tulu YI, Peuler BD, Wilkens DA. Air-flow geometry in air sparging of fine-grained sands. Hydrogeol J 2001;9:168–76. http:// dx.doi.org/10.1007/s100400000104. [134] Selker JS, Niemet M, Mcduffie NG, Gorelick SM, Parlange J-Y. The Local Geometry of Gas Injection into Saturated Homogeneous Porous Media. Transp Porous Media 2006;68:107–27. http://dx.doi.org/10.1007/s11242006-0005-0. [135] Heron G, Gierke JS, Faulkner B, Mravik S, Wood L, Enfield CG. Pulsed air sparging in aquifers contaminated with dense nonaqueous phase liquids. B. Agaoglu et al. / Advances in Water Resources 79 (2015) 162–194 [136] [137] [138] [139] [140] [141] [142] [143] [144] [145] [146] [147] [148] [149] [150] [151] [152] [153] [154] [155] [156] [157] [158] [159] [160] [161] [162] Ground Water Monit Remediat 2002, http://dx.doi.org/10.1111/j.1745-6592. 2002.tb00773.x. Berkey JS, Lachmar TE, Doucette WJ, Ryan Dupont R. Tracer studies for evaluation of in situ air sparging and in-well aeration system performance at a gasoline-contaminated site. J Hazard Mater 2003;98:127–44. http:// dx.doi.org/10.1016/S0304-3894(02)00293-5. Thomson NR, Johnson RL. Air distribution during in situ air sparging: an overview of mathematical modeling. J Hazard Mater 2000;72:265–82. http:// dx.doi.org/10.1016/S0304-3894(99)00143-0. Waduge WAP, Soga K, Kawabata J. Effect of NAPL entrapment conditions on air sparging remediation efficiency. J Hazard Mater 2004;110:173–83. http:// dx.doi.org/10.1016/j.jhazmat.2004.02.050. Jang W, Aral MM. Multiphase flow fields in in-situ air sparging and its effect on remediation. Transp Porous Media 2008;76:99–119. http://dx.doi.org/ 10.1007/s11242-008-9238-4. Hunt JR, Sitar N, Udell KS. Nonaqueous phase liquid transport and cleanup: 2. Experimental studies. Water Resour Res 1988;24:1259–69. http://dx.doi.org/ 10.1029/WR024i008p01259. Schmidt R, Betz C, Faerber A. LNAPL and DNAPL behaviour during steam injection into the unsaturated zone. IAHS Publ Assoc Hydrol Sci 1998:111–7. She HY, Sleep BE. Removal of perchloroethylene from a layered soil system by steam flushing. Ground Water Monit Rem 1999;19:70–7. http://dx.doi.org/ 10.1111/j.1745-6592.1999.tb00207.x. Kaslusky SF, Udell KS. A theoretical model of air and steam co-injection to prevent the downward migration of DNAPLs during steam-enhanced extraction. J Contam Hydrol 2002;55:213–32. http://dx.doi.org/10.1016/ S0169-7722(01)00191-7. Ji W, Dahmani A, Ahlfeld DP, Lin JD, Hill E. Laboratory study of air sparging: air flow visualization. Groundw Monit Remediat 1993;13:115–36. http:// dx.doi.org/10.1111/j.1745-6592.1993.tb00455.x. Lundegard PD, Andersen G. Multiphase numerical simulation of air sparging performance. Groundwater 1996;34(3):451–60. http://dx.doi.org/10.1111/ j.1745-6584.1996.tb02026.x. McCray JE, Falta RW. Defining the air sparging radius of influence for groundwater remediation. J Contam Hydrol 1996;24:25–52. http:// dx.doi.org/10.1016/0169-7722(96)00005-8. FaisalAnwar AHM, Bettahar M, Matsubayashi U. A method for determining air–water interfacial area in variably saturated porous media. J Contam Hydrol 2000;43:129–46. http://dx.doi.org/10.1016/S0169-7722(99)00103-5. Peng S, Brusseau ML. Impact of soil texture on air–water interfacial areas in unsaturated sandy porous media. Water Resour Res 2005;41:89. http:// dx.doi.org/10.1029/2004WR003233. Rao PSC, Annable MD, Kim H. NAPL source zone characterization and remediation technology performance assessment: recent developments and applications of tracer techniques. J Contam Hydrol 2000;45:63–78. http:// dx.doi.org/10.1016/S0169-7722(00)00119-4. Kim H, Rao PSC, Annable MD. Determination of effective air–water interfacial area in partially saturated porous media using surfactant adsorption. Water Resour Res 1997;33:2705–11. http://dx.doi.org/10.1029/97WR02227. Schaefer CE, DiCarlo DA, Blunt MJ. Determination of water–oil interfacial area during 3-phase gravity drainage in porous media. J Colloid Interface Sci 2000;221:308–12. http://dx.doi.org/10.1006/jcis.1999.6604. Cho J, Annable MD. Characterization of pore scale NAPL morphology in homogeneous sands as a function of grain size and NAPL dissolution. http://dx.doi.org/10.1016/ Chemosphere 2005;61:899–908. j.chemosphere.2005.04.042. Brusseau ML, Peng S, Schnaar G, Murao A. Measuring air–water interfacial areas with X-ray microtomography and interfacial partitioning tracer tests. http://dx.doi.org/10.1021/ Environ Sci Technol 2007;41:1956–61. es061474m. Kim H, Choi K-M, Moon J-W, Annable MD. Changes in air saturation and air– water interfacial area during surfactant-enhanced air sparging in saturated sand. J Contam Hydrol 2006;88:23–35. http://dx.doi.org/10.1016/ j.jconhyd.2006.05.009. Schaefer CE, Callaghan AV, King JD, Mccray JE. Dense nonaqueous phase liquid architecture and dissolution in discretely fractured sandstone blocks. Environ Sci Technol 2009;43:1877–83. http://dx.doi.org/10.1021/es8011172. Johns ML, Gladden LF. Magnetic resonance imaging study of the dissolution kinetics of octanol in porous media. J Colloid Interface Sci 1999;210:261–70. http://dx.doi.org/10.1006/jcis.1998.5950. Johns ML, Gladden LF. Probing ganglia dissolution and mobilization in a water-saturated porous medium using MRI. J Colloid Interface Sci 2000;225:119–27. http://dx.doi.org/10.1006/jcis.2000.6742. Leverett MC. Capillary behaviour in porous solids. Trans AIME 1941:152–70. Morrow NR. Physics and thermodynamics of capillary action in porous media. Ind Eng Chem Res 1970;62:32–56. Bradford SA, Leij FJ. Estimating interfacial areas for multi-fluid soil systems. J http://dx.doi.org/10.1016/S0169Contam Hydrol 1997;27:83–105. 7722(96)00048-4. Dalla E, Hilpert M, Miller CT. Computation of the interfacial area for two-fluid porous medium systems. J Contam Hydrol 2002;56:25–48. http://dx.doi.org/ 10.1016/S0169-7722(01)00202-9. Hilpert M, Miller CT. Pore-morphology-based simulation of drainage in totally wetting porous media. Adv Water Resour 2001;24:243–55. http:// dx.doi.org/10.1016/S0309-1708(00)00056-7. 191 [163] Dobson R, Schroth MH, Oostrom M, Zeyer J. Determination of NAPL–water interfacial areas in well-characterized porous media. Environ Sci Technol 2006;40:815–22. http://dx.doi.org/10.1021/es050037p. [164] Oostrom M, White MD, Brusseau ML. Theoretical estimation of free and entrapped nonwetting–wetting fluid interfacial areas in porous media. Adv Water Resour 2001;24. http://dx.doi.org/10.1016/S0309-1708(01)00017-3. [165] Grant GP, Gerhard JI. Simulating the dissolution of a complex dense nonaqueous phase liquid source zone: 1 Model to predict interfacial area. Water Resour Res 2007;43. http://dx.doi.org/10.1029/2007WR006038. [166] Brusseau ML, Peng S, Schnaar G, Costanza-Robinson MS. Relationships among air–water interfacial area, capillary pressure, and water saturation for a sandy porous medium. Water Resour Res 2006;42. http://dx.doi.org/10.1029/ 2005WR004058. [167] Fenwick DH, Blunt MJ. Three-dimensional modeling of three phase imbibition and drainage. Adv Water Resour 1998;21:121–43. http:// dx.doi.org/10.1016/S0309-1708(96)00037-1. [168] Keller AA, Blunt MJ, Roberts APV. Micromodel observation of the role of oil layers in three-phase flow. Transp Porous Media 1997:277–97. http:// dx.doi.org/10.1023/A:1006589611884. [169] Hassanizadeh SM, Gray WG. Thermodynamic basis of capillary pressure in porous media. Water Resour Res 1993;29(10):3389–405. http://dx.doi.org/ 10.1029/93WR01495. [170] Reeves PC, Celia MA. A functional relationship between capillary pressure, saturation, and interfacial area as revealed by a pore-scale network model. Water Resour Res 1996;32:2345–58. http://dx.doi.org/10.1029/96WR01105. [171] Held RJ, Celia MA. Modeling support of functional relationships between capillary pressure, saturation, interfacial area and common lines. Adv Water Resour 2001;24:325–43. http://dx.doi.org/10.1016/S0309-1708(00)00060-9. [172] Joekar-Niasar V, Hassanizdeh SM, Leijnse A. Insights into the relationships among capillary pressure, saturation, interfacial area and relative permeability using pore-network modeling. Transp Porous Media 2008;74:201–19. http://dx.doi.org/10.1007/s11242-007-9191-7. [173] Cheng JT, Pyrak-Nolte LJ, Nolte DD, Giordano NJ. Linking pressure and saturation through interfacial areas in porous media. Geophys Res Lett 2004;31. http://dx.doi.org/10.1029/2003GL019282. [174] Chen D, Pyrak-Nolte LJ, Griffin J, Giordano NJ. Measurement of interfacial area per volume for drainage and imbibition. Water Resour Res 2007;43(12). http://dx.doi.org/10.1029/2007WR006021. [175] Chen L, Kibbey TCG, Street WB. Measurement of air–water interfacial area for multiple hysteretic drainage curves in an unsaturated fine sand. 2006:6874–80. http://dx.doi.org/10.1029/2004WR003278. (1). [176] Joekar-Niasar V, Hassanizadeh SM, Pyrak-Nolte LJ, Berentsen C. Simulating drainage and imbibition experiments in a high-porosity micromodel using an unstructured pore network model. Water Resour Res 2009;45:1–15. http:// dx.doi.org/10.1029/2007WR006641. [177] Porter ML, Schaap MG, Wildenschild D. Lattice-Boltzmann simulations of the capillary pressure–saturation-interfacial area relationship for porous media. Adv Water Resour 2009;32:1632–40. http://dx.doi.org/10.1016/j.advwatres. 2009.08.009. [178] Porter ML, Wildenschild D, Grant G, Gerhard JI. Measurement and prediction of the relationship between capillary pressure, saturation, and interfacial area in a NAPL–water–glass bead system. Water Resour Res 2010;46:1–10. http://dx.doi.org/10.1029/2009WR007786. [179] Lenordmand R, Toulboul E, Zarcone C. Numerical models and experiments on immiscible displacements in porous media. J Fluid Mech 1988;189. http:// dx.doi.org/10.1017/S0022112088000953. [180] Lenordmand R, Zarcone C, Sarr A. Mechanisms of the displacement of one fluid by another in a network of capillary ducts. J Fluid Mech 1983;135:337–53. http://dx.doi.org/10.1017/S0022112083003110. [181] Zhang C, Werth CJ, Webb AG. A magnetic resonance imaging study of dense nonaqueous phase liquid dissolution from angular porous media. Environ Sci Technol 2002;36:3310–7. http://dx.doi.org/10.1021/es011497v. [182] Blunt MJ, Scher H. Pore-level modeling of wetting. Phys Rev E 1995;52:6387–403. http://dx.doi.org/10.1103/PhysRevE.52.6387. [183] Hughes RG, Blunt MJ. Pore scale modeling of rate effects in imbibition. Transp http://dx.doi.org/10.1023/ Porous Media 2000;40:295–322. A:1006629019153. [184] Chen JD, Wilkinson D. Pore-scale viscous fingering in porous media. Phys Rev Lett 1985;55. http://dx.doi.org/10.1103/PhysRevLett.55.1892. [185] Sinha PK, Wang C. Pore-network modeling of liquid water transport in gas diffusion layer of a polymer electrolyte fuel cell. Electrochim Acta 2007;52:7936–45. http://dx.doi.org/10.1016/j.electacta.2007.06.061. [186] Powers SE, Anckner WH, Seacord TF. Wettability of NAPL-contaminated sands. J Environ Eng 1996:889–96. http://dx.doi.org/10.1061/(ASCE)07339372(1996)122:10(889). [187] Kovscek AR, Wong H, Radke CJ. A pore level scenario for the development of mixed-wettability in oil reservoirs. US Dep Energy; 1992. [188] Valvatne PH, Blunt MJ. Predictive pore-scale modeling of two-phase flow in mixed wet media. Water Resour Res 2004;40. http://dx.doi.org/10.1029/ 2003WR002627. [189] Bradford SA, Abriola LM, Leij FJ. Wettability effects on two- and three-fluid relative permeabilities. J Contam Hydrol 1997;28:171–91. http://dx.doi.org/ 10.1016/S0169-7722(97)00024-7. [190] Hwang II S, Lee KP, Lee DS, Powers SE. Effects of fractional wettability on capillary pressure–saturation–relative permeability relations of two-fluid 192 [191] [192] [193] [194] [195] [196] [197] [198] [199] [200] [201] [202] [203] [204] [205] [206] [207] [208] [209] [210] [211] [212] [213] [214] [215] [216] B. Agaoglu et al. / Advances in Water Resources 79 (2015) 162–194 systems. Adv Water Resour 2006;29:212–26. http://dx.doi.org/10.1016/ j.advwatres.2005.03.020. Mayer AS, Miller CT. The influence of porous medium characteristics and measurement scale on pore-scale distributions of residual nonaqueous-phase liquids. J Contam Hydrol 1992;11. http://dx.doi.org/10.1016/01697722(92)90017-9. Schnaar G, Brusseau ML. Pore-scale characterization of organic immiscibleliquid morphology in natural porous media using synchrotron X-ray microtomography. Environ Sci Technol 2005;39:8403–10. http://dx.doi.org/ 10.1021/es0508370. Knutson CE, Werth CJ, Valocchi AJ. Pore-scale modeling of dissolution from variably distributed nonaqueous phase liquid blobs. Water Resour Res 2001;37:2951–63. http://dx.doi.org/10.1029/2001WR000587. Chomsurin C, Werth CJ. Analysis of pore-scale nonaqueous phase liquid dissolution in etched silicon pore networks. Water Resour Res 2003;39. http://dx.doi.org/10.1029/2002WR001643. Reddy KR, Adams JA. Effects of soil heterogeneity on airflow patterns and hydrocarbon removal during in situ air sparging. J Geotech 2001:234–47. http://dx.doi.org/10.1061/(ASCE)1090-0241(2001)127:3(234). Peterson JW, Lepczyk PA, Lake KL. Effect of sediment size on area of influence during groundwater remediation by air sparging: a laboratory approach. Environ Geol 1999;38:1–6. http://dx.doi.org/10.1007/s002540050394. Braida W, Ong SK. Influence of porous media and airflow rate on the fate of NAPLs under air sparging. Transp Porous Media 2000:29–42. http:// dx.doi.org/10.1023/A:1006603031438. Rogers SW, Ong SK. Influence of porous media, airflow rate, and air channel spacing on benzene NAPL removal during air sparging. Environ Sci Technol 2000;34:764–70. http://dx.doi.org/10.1021/es9901112. Russo AE, Narter M, Brusseau ML. Characterizing pore-scale dissolution of organic immiscible liquid in a poorly-sorted natural porous medium. Environ Sci Technol 2009:5671–8. http://dx.doi.org/10.1021/es803158x. Hunt JR, Sitar N, Udell KS. Nonaqueous phase liquid transport and cleanup: 1. Analysis of mechanisms. Water Resour Res 1988;24:1247–58. http:// dx.doi.org/10.1029/WR024i008p01247. Schwille F. Petroleum contamination of the subsoil – a hydrological problem. In: Hepple P, editor. The joint problems of the oil and water industries, proceedings of a symposium. London: Inst Pet; 1967. Illangasekare TH, Armbruster III EJ, Yates DN. Non-aqueous-phase fluids in heterogeneous aquifers-experimental study. J Environ Eng http://dx.doi.org/10.1061/(ASCE)07331995;121(8):571–9. 9372(1995)121:8(571). Imhoff PT, Arthur MH, Miller CT. Complete dissolution of trichloroethylene in saturated porous media. Environ Sci Technol 1998;32:2417–24. http:// dx.doi.org/10.1021/es970965r. Mackay DM, Cherry JA. Groundwater contamination: pump-and-treat remediation. Environ Sci Technol 1989;23(6):630–6. http://dx.doi.org/ 10.1021/es00064a001. Ball WP, Roberts PV. Diffusive rate limitations in the sorption of organic chemicals. In: Baker R, editor. Organic substances and sediments in water. Process anal, vol. 2. Chelsea, MI: Lewis Publishers Inc; 1991. p. 273–310 [chapter 13]. Luthy RG, Aiken GR, Brusseau ML, Cunningham SD, Gschwend PM, Pignatello JJ, et al. Sequestration of hydrophobic organic contaminants by geosorbents. Environ Sci Technol 1997:3341–7. http://dx.doi.org/10.1021/es970512m. Grathwohl P. Diffusion in natural porous media: contaminant transport, sorption/desorption and dissolution kinetics. Kluwer Academic Publishers; 1998. Kueper BH, Frind EO. An overview of immiscible fingering in porous media. J Contam Hydrol 1988;2:95–110. http://dx.doi.org/10.1016/01697722(88)90001-0. Mackay DM, Roberts PV, Cherry JA. Transport of organic contaminants in groundwater. Environ Sci Technol 1985;19:384–92. http://dx.doi.org/ 10.1021/es00135a001. Palomino AM, Grubb DG. Recovery of dodecane, octane and toluene spills in sandpacks using ethanol. J Hazard Mater 2004;110:39–51. http://dx.doi.org/ 10.1016/j.jhazmat.2004.02.035. Grubb DG, Sitar N. Mobilization of trichloroethene (TCE) during ethanol flooding in uniform and layered sand packs under confined conditions. Water Resour Res 1999;35:3275–89. http://dx.doi.org/10.1029/1999WR900222. Grubb DG, Sitar N. Horizontal ethanol floods in clean, uniform, and layered sand packs under confined conditions. Water Resour Res 1999;35:3291–302. http://dx.doi.org/10.1029/1999WR900223. Xu X, Chen S, Zhang D. Convective stability analysis of the long-term storage of carbon dioxide in deep saline aquifers. Adv Water Resour 2006;29:397–407. http://dx.doi.org/10.1016/j.advwatres.2005.05.008. Kneafsey TJ, Pruess K. Laboratory flow experiments for visualizing carbon dioxide-induced. density-driven brine convection. Transp Porous Media 2009;82:123–39. http://dx.doi.org/10.1007/s11242-009-9482-2. Hatfield K, Stauffer TB. Transport in porous media containing residual hydrocarbon. I: Model. J Environ Eng 1993;119:540–58. http://dx.doi.org/ 10.1061/(ASCE)0733-9372(1993)119:3(540). Soerens TS, Sabatini DA, Harwell JH. Effects of flow bypassing and nonuniform NAPL distribution on the mass transfer characteristics of NAPL dissolution. Water Resour Res 1998;34:1657–73. http://dx.doi.org/10.1029/ 98WR00554. [217] Chrysikopoulos CV, Voudrias EA, Fyrillas MM. Modeling of contaminant transport resulting from dissolution of nonaqueous phase liquid pools in saturated porous media. Transp Porous Media 1994;16:125–45. http:// dx.doi.org/10.1007/BF00617548. [218] Chrysikopoulos CV. Three-dimensional analytical models of contaminant transport from nonaqueous phase liquid pool dissolution in saturated subsurface formations. Water Resour Res 1995;31:1137–45. http:// dx.doi.org/10.1029/94WR02780. [219] Holman HYN, Javandel I. Evaluation of transient dissolution of slightly watersoluble compounds from a light nonaqueous phase liquid pool. Water Resour Res 1996;32:915–23. http://dx.doi.org/10.1029/96WR00075. [220] Sale TC, McWhorter DB. Steady state mass transfer from single-component dense nonaqueous phase liquids in uniform flow fields. Water Resour Res 2001;37:393–404. http://dx.doi.org/10.1029/2000WR900236. [221] Hunt B. Dispersion sources in uniform groundwater flow. Hydraul Div Am Soc Civ Eng 1960;104:75–85. [222] Falta RW, Suresh Rao P, Basu N. Assessing the impacts of partial mass depletion in DNAPL source zones. I. Analytical modeling of source strength functions and plume response. J Contam Hydrol 2005;78:259–80. http:// dx.doi.org/10.1016/j.jconhyd.2005.05.010. [223] Lemke LD, Abriola LM, Lang JR. Influence of hydraulic property correlation on predicted dense nonaqueous phase liquid source zone architecture, mass recovery and contaminant flux. Water Resour Res 2004;40. http://dx.doi.org/ 10.1029/2004WR003061. [224] Christ JA, Ramsburg CA, Pennell KD, Abriola LM. Predicting DNAPL mass discharge from pool-dominated source zones. J Contam Hydrol 2010;114:18–34. http://dx.doi.org/10.1016/j.jconhyd.2010.02.005. [225] Rao PSC, Jawitz JW, Enfield CG, Falta RW, Annable MD, Wood AL. Technology integration for contaminated site remediation: clean-up goals and performance criteria. Groundwater Qual 2001:571–8. [226] DiFilippo EL, Carroll KC, Brusseau ML. Impact of organic-liquid distribution and flow-field heterogeneity on reductions in mass flux. J Contam Hydrol 2010;115:14–25. http://dx.doi.org/10.1016/j.jconhyd.2010.03.002. [227] Falta RW. REMChlor user’s manual beta version 1.0. n.d. [228] Jawitz JW, Fure AD, Demmy GG, Berglund S, Rao PSC. Groundwater contaminant flux reduction resulting from nonaqueous phase liquid mass reduction. Water Resour Res 2005;41. http://dx.doi.org/10.1029/ 2004WR003825. [229] Fure AD, Jawitz JW, Annable MD. DNAPL source depletion: linking architecture and flux response. J Contam Hydrol 2006;85:118–40. http:// dx.doi.org/10.1016/j.jconhyd.2006.01.002. [230] Chen X, Jawitz JW. Reactive tracer tests to predict dense nonaqueous phase liquid dissolution dynamics in laboratory flow chambers. Environ Sci Technol 2008;42:5285–91. [231] Wood AL, Enfield CG, Espinoza FP, Annable M, Brooks MC, Rao PSC, et al. Design of aquifer remediation systems: (2) estimating site-specific performance and benefits of partial source removal. J Contam Hydrol 2005;81:148–66. http://dx.doi.org/10.1016/j.conhyd.2005.08.004. [232] Braida W, Ong SK. Modeling of air sparging of VOC-contaminated soil columns. J Contam Hydrol 2000;41:385–402. http://dx.doi.org/10.1016/ S0169-7722(99)00075-3. [233] Chao KP, Ong SK, Protopapas A. Water-to-air mass transfer of VOCs: laboratory-scale air sparging system. J Environ Eng 1998;124(11):1054–60. http://dx.doi.org/10.1061/(ASCE)0733-9372(1998)124:11(1054). [234] Jia C, Shing K, Yortsos YC. Visualization and simulation of non-aqueous phase liquids solubilization in pore networks. J Contam Hydrol 1999;35:363–87. http://dx.doi.org/10.1016/S0169-7722(98)00102-8. [235] Celia MA, Reeves PC, Ferrand LA. Recent advances in pore scale models for multiphase flow in porous media. Rev Geophys 1995:1049–57. http:// dx.doi.org/10.1029/95RG00248. [236] Joekar-Niasar V, Prodanović M, Wildenshild D, Hassanizadeh SM. Network model investigation of interfacial area, capillary pressure and saturation relationships in granular porous media. Water Resour Res 2010;46:1–18. http://dx.doi.org/10.1029/2009WR008585. [237] Al-Gharbi M, Blunt MJ. Dynamic network modeling of two-phase drainage in porous media. Phys Rev E 2005;71:016308. http://dx.doi.org/10.1103/ PhysRevE.71.016308. [238] Blunt MJ. Flow in porous media – pore-network models and multiphase flow. Curr Opin Colloid Interface Sci 2001;6(3):197–207. http://dx.doi.org/ 10.1016/S1359-0294(01)00084-X. [239] Wilkinson D, Willemsen JF. Invasion percolation: a new form of percolation theory. J Phys A Math Gen 1983;16:3365–76. http://dx.doi.org/10.1088/ 0305-4470/16/14/028. [240] Parker JC, Katyal AK, Kaluarachchi JJ, Lenhard RJ, Johnson TJ, Jayaraman K, et al. Modeling multiphase organic chemical transport in soils and ground water. EPA/600/2-91/042, US Environ Prot Agency, Washington, DC; 1991. [241] Zhao W, Ioannidis MA. Effect of NAPL film stability on the dissolution of residual wetting NAPL in porous media: a pore-scale modeling study. Adv Water Resour 2007;30:171–81. http://dx.doi.org/10.1016/ j.advwatres.2005.03.025. [242] Zhao W, Ioannidis MA. Convective mass transfer across fluid interfaces in straight angular pores. Transp Porous Media 2007:495–509. http:// dx.doi.org/10.1007/s11242-006-0025-9. B. Agaoglu et al. / Advances in Water Resources 79 (2015) 162–194 [243] Martys NS, Chen H. Simulation of multicomponent fluids in complex threedimensional geometries by the lattice Boltzmann method. Phys Rev E 1996;53. [244] Swift MR, Orlandini E, Osborn WR, Yeomans JM. Lattice Boltzmann simulations of liquid–gas and binary fluid systems. Phys Rev E Stat Phys Plasmas Fluids Relat Interdisciplin Top 1996;54:5041–52. http://dx.doi.org/ 10.1103/PhysRevE.54.5041. [245] Pan C, Hilpert M, Miller CT. Lattice-Boltzmann simulation of two-phase flow in porous media. Water Resour Res 2004;40. http://dx.doi.org/10.1029/ 2003WR002120. [246] Schaap MG, Porter ML. Comparison of pressure–saturation characteristics derived from computed tomography and lattice Boltzmann simulations. Water Resour Res 2007;43:1–15. http://dx.doi.org/10.1029/ 2006WR005730. [247] Meakin P, Tartakovsky AM. Modeling and simulation of pore-scale multiphase fluid flow and reactive transport in fractured and porous media. Rev Geophys 2009:1–47. http://dx.doi.org/10.1029/ 2008RG000263.1.INTRODUCTION. [248] Yuan P, Schaefer L. Equations of state in a lattice Boltzmann model. Phys Fluids 2006;18:042101. http://dx.doi.org/10.1063/1.2187070. [249] Bao J, Schaefer L. Lattice Boltzmann equation model for multi-component multi-phase flow with high density ratios. Appl Math Model 2013;37:1860–71. http://dx.doi.org/10.1016/j.apm.2012.04.048. [250] Blunt MJ, Bijeljic B, Dong H, Gharbi O, Iglauer S, Mostaghimi P, et al. Advances in water resources pore-scale imaging and modelling. Adv Water Resour 2013;51:197–216. http://dx.doi.org/10.1016/j.advwatres.2012.03.003. [251] Raeini AQ, Blunt MJ, Bijeljic B. Modelling two-phase flow in porous media at the pore scale using the volume-of-fluid method. J Comput Phys 2012;231:5653–68. http://dx.doi.org/10.1016/j.jcp.2012.04.011. [252] Raeini AQ. Modelling multiphase flow through micro-CT images of the pore space [Ph.D. Diss.]. Imperial College London; 2013:1–152. [253] Kim TJ, Chrysikopoulos CV. Mass transfer correlations for nonaqueous phase liquid pool dissolution in saturated porous media. Water Resour Res 1999;35:449–59. http://dx.doi.org/10.1029/1998WR900053. [254] Schaerlaekens J, Vanderborght J, Merckx R, Feyen J. Surfactant enhanced solubilization of residual trichloroethene: an experimental and numerical analysis. J Contam Hydrol 2000;46:1–16. [255] Ji W, Brusseau ML. A general mathematical model for chemical-enhanced flushing of soil contaminated by organic compounds. Water Resour Res 1998;34:1635. http://dx.doi.org/10.1029/98WR01040. [256] Geller JT, Hunt JR. Mass transfer from nonaqueous phase organic liquids in water-saturated porous media. Water Resour Res 1993;29:833–45. http:// dx.doi.org/10.1029/92WR02581. [257] Powers SE, Nambi IM, Curry GW. Non-aqueous phase liquid dissolution in heterogeneous systems: mechanisms and a local equilibrium modeling approach. Water Resour Res 1998;34:3293–302. http://dx.doi.org/10.1029/ 98WR02471. [258] Corey AT. The interrelation between gas and oil relative permeabilities. Prod Mon 1954;19:38–41. [259] Brusseau ML, Zhang Z, Nelson NT, Cain RB, Tick GR, Oostrom M. Dissolution of nonuniformly distributed immiscible liquid: intermediate-scale experiments and mathematical modeling. Environ Sci Technol 2002;36:1033–41. http:// dx.doi.org/10.1021/es010609f. [260] Miller CT, Gleyzer SN, Imhoff PT. Numerical modeling of NAPL dissolution fingering in porous media. In: Selim HM, Ma L, editors. Physical nonequilibrium in soils: modeling and application. Mich: Ann Arbor Press Ann Arbor; 1998. p. 389–415. [261] Mason AR, Kueper BH. Numerical simulation of surfactant-enhanced solubilization of pooled DNAPL. Environ Sci Technol 1996;30:3205–15. http://dx.doi.org/10.1021/es9507372. [262] Rathfelder KM, Abriola LM, Taylor TP, Pennell KD. Surfactant enhanced recovery of tetrachloroethylene from a porous medium containing low permeability lenses. 2. Numerical simulation. J Contam Hydrol 2001;48:351–74. http://dx.doi.org/10.1016/S0169-7722(00)00186-8. [263] Abriola LM, Lang JR, Rathfelder KM. Michigan soil vapor extraction remediation MISER. Model: a computer program to model soil vapor extraction and bioventing of organic chemicals in unsaturated geologic materials. US Environ Prot Agency, EPAr600rR-97r099, Natl Risk Manag Res Lab Ada, OK; 1997. [264] Taylor TP, Pennell KD, Abriola LM, Dane JH. Surfactant enhanced recovery of tetrachloroethylene from a porous medium containing low permeability lenses. 1. Experimental studies. J Contam Hydrol 2001;48:325–50. http:// dx.doi.org/10.1016/S0169-7722(00)00185-6. [265] Abriola LM, Pinder GF. A multiphase approach to the modeling of porous media contamination by organic compounds: 1. Equation development. Water Resour Res 1985;21:11–8. http://dx.doi.org/10.1029/ WR021i001p00011. [266] Kueper BH, Frind EO. Two-phase flow in heterogeneous porous media: 1. Model development. Water Resour Res 1991;27:1049–57. http://dx.doi.org/ 10.1029/91WR00266. [267] Faust CR. Transport of immiscible fluids within and below the unsaturated zone: a numerical model. Water Resour Res 1985;21:587–96. http:// dx.doi.org/10.1029/WR021i004p00587. [268] Osborne M, Sykes J. Numerical modeling of immiscible organic transport at the Hyde Park landfill. Water Resour Res 1986;22:25–33. http://dx.doi.org/ 10.1029/WR022i001p00025. 193 [269] Sleep BE, Sykes JF. Compositional simulation of groundwater contamination by organic compounds: 1. Model development and verification. Water Resour Res 1993;29:1697–708. http://dx.doi.org/10.1029/93WR00283. [270] White MD, Oostrom M, Imhard J. Modeling fluid flow and transport in variably saturated porous media with the STOMP simulator. 1. Nonvolatile three-phase model description. Adv Water Resour 1995;18:353–64. http:// dx.doi.org/10.1016/0309-1708(95)00018-E. [271] Rathfelder K, Abriola LM. Mass conservative numerical solutions of the headbased Richards equation. Water Resour Res 1994;30:2579–86. http:// dx.doi.org/10.1029/94WR01302. [272] Guiguer N, Frind EO. Dissolution and mass transfer processes for residual organics in the saturated groundwater zone. In: Dracos Stauffer, editor. Transport and reactive processes in aquifers. Rotterdam: Balkema; 1994. p. 475–80. [273] Dekker TJ, Abriola LM. The influence of field-scale heterogeneity on the infiltration and entrapment of dense nonaqueous phase liquids in saturated formations. J Contam Hydrol 2000;42:187–218. http://dx.doi.org/10.1016/ S0169-7722(99)00092-3. [274] Hinkelmann R, Breiting TR, Helmig R. Modeling of Hydrosystems with MUFTE-UG: multiphase flow and transport processes in the subsurface. In: Fourth ınt conf hydroinformatics, Iowa, USA; 2000. [275] Hinkelmann R, Sheta H, Class H, Sauter EJ, Helmig R, Schlüter M. Numerical simulation of freshwater, salt water and methane interaction processes in a coastal aquifer. In: Resource recovery, confinement, and remediation of environmental hazards. New York: Springer; 2002. p. 263–81. [276] Helmig R, Class H, Huber R, Sheta H, Ewing R, Hinkelmann R, et al. Architecture of the modular program system MUFTE-UG for simulating multiphase flow and transport processes in heterogeneous porous media. Math Geol 1998;64:123–31. [277] Hinkelmann R. Efficient numerical methods and information-processing techniques for modeling hydro-and environmental systems. Berlin: Springer; 2005. p. 21. [278] Falta RW. Modeling sub-grid-block-scale dense nonaqueous phase liquid (DNAPL) pool dissolution using a dual-domain approach. Water Resour Res 2003;39. http://dx.doi.org/10.1029/2003WR002351. [279] Falta RW, Pruess K, Finsterle S, Battistelli A. T2VOC user’s guide. Rep LBL36400, Lawrence Berkeley Natl Lab, Berkeley, CA; March 1995. [280] Basu NB, Fure AD, Jawitz JW. Predicting dense nonaqueous phase liquid dissolution using a simplified source depletion model parameterized with partitioning tracers. Water Resour Res 2008;44:1–13. http://dx.doi.org/ 10.1029/2007WR006008. [281] Basu N, Rao P, Falta R, Annable M, Jawitz J, Hatfield K. Temporal evolution of DNAPL source and contaminant flux distribution: Impacts of source mass depletion. J Contam Hydrol 2008;95:93–109. http://dx.doi.org/10.1016/ j.jconhyd.2007.08.001. [282] O’Carroll DM, Sleep BE. Role of NAPL thermal properties in the effectiveness of hot water flooding. Transp Porous Media 2009;79:393–405. http:// dx.doi.org/10.1007/s11242-008-9329-2. [283] Reitsma S, Kueper BH. Non-equilibrium alcohol flooding model for immiscible phase remediation: 2. Model development and application. Adv http://dx.doi.org/10.1016/S0309Water Resour 1998;21:663–78. 1708(97)00025-0. [284] Gray WG, Celia MA, Reeves PC. Incorporation of interfacial areas in models of two-phase flow. Int keamey found soil sci conf vadose zo hydrol – cut across discip. Univ California, Davis; Sept 6–8 1995. [285] Li B, Fu J. Interfacial tensions of two-liquid-phase ternary systems. J Chem Eng Data 1992;37:172–4. http://dx.doi.org/10.1021/je00006a009. [286] Roeder E, Falta RW. Modeling unstable alcohol flooding of DNAPLcontaminated columns. Adv Water Resour 2001;24. http://dx.doi.org/ 10.1016/S0309-1708(00)00072-5. [287] Liang H, Falta RW. Modeling field-scale cosolvent flooding for DNAPL source zone remediation. J Contam Hydrol 2008;96:1–16. http://dx.doi.org/10.1016/ j.jconhyd.2007.09.005. [288] Kaye AJ, Cho J, Basu NB, Chen X. Laboratory investigation of flux reduction from dense non-aqueous phase liquid (DNAPL) partial source zone remediation by enhanced dissolution. J Contam Hydrol 2008;102:17–28. http://dx.doi.org/10.1016/j.jconhyd.2008.01.006. [289] Lingineni S, Dhir VK. Controlling transport processes during NAPL removal by soil venting. Adv Water Resour 1997;20:157–69. http://dx.doi.org/10.1016/ S0309-1708(96)00028-0. [290] Hein GL, Gierke JS, Hutzler NJ, Falta RW. Three-dimensional experimental testing of a two-phase flow-modeling approach for air sparging. Ground Water Monit Rem 1997;17(3):222–30. http://dx.doi.org/10.1111/j.17456592.1997.tb00597.x. [291] McClure PD, Sleep BE. Simulation of bioventing for soil and ground-water remediation. J Environ Eng 1996:1003–12. http://dx.doi.org/10.1061/ (ASCE)0733-9372(1996)122:11(1003). [292] Unger AJA, Sudicky EA, Forsyth PA. Mechanisms controlling vacuum extraction coupled with air sparging for remediation of heterogeneous formations contaminated by dense nonaqueous phase liquids. Water Resour Res 1995;31:1913–25. http://dx.doi.org/10.1029/95WR00172. [293] Yoon H, Werth CJ, Valocchi AJ, Oostrom M. Impact of nonaqueous phase liquid (NAPL) source zone architecture on mass removal mechanisms in strongly layered heterogeneous porous media during soil vapor extraction. J Contam Hydrol 2008;100:58–71. http://dx.doi.org/10.1016/ j.jconhyd.2008.05.006. 194 B. Agaoglu et al. / Advances in Water Resources 79 (2015) 162–194 [294] Xu T, Apps JA, Pruess K. Mineral sequestration of carbon dioxide in a sandstone–shale system. Chem Geol 2005;217(3):295–318. http:// dx.doi.org/10.1016/j.chemgeo.2004.12.015. [295] Xu T, Sonnenthal E, Spycher N, Pruess K. TOUGHREACT user’s guide: a simulation program for non-isothermal multiphase reactive geochemical transport in variably saturated geologic media, V1; 2008. [296] Thibeau S, Dutin A. Large scale CO2 storage in unstructured aquifers: Modeling study of the ultimate CO2 migration distance. Energy Proc 2011;4:4230–7. http://dx.doi.org/10.1016/j.egypro.2011.02.371. [297] Niessner J, Hassanizadeh SM. A model for two-phase flow in porous media including fluid–fluid interfacial area. Water Resour Res 2008;44(8). http:// dx.doi.org/10.1029/2007WR006721. [298] Hassanizadeh SM, Gray WG. Mechanics and thermodynamics of multiphase flow in porous media including interphase boundaries. Adv Water Resour 1990;13(4):169–86. http://dx.doi.org/10.1016/0309-1708(90)90040-B. [299] Jackson AS, Miller CT, Gray WG. Thermodynamically constrained averaging theory approach for modeling flow and transport phenomena in porous medium systems: 6. Two-fluid-phase flow. Adv Water Resour 2009;32:779–95. http://dx.doi.org/10.1016/j.advwatres.2008.11.010. [300] Miller CT, Dawson CN, Farthing MW, Hou TY, Huang J, Kees CE, et al. Numerical simulation of water resources problems: models, methods, and trends. Adv Water Resour 2013;51:405–37. http://dx.doi.org/10.1016/ j.advwatres.2012.05.008. [301] Jackson ABS. Multiscale modeling of multiphase flow in porous media using the thermodynamically constrained averaging theory approach [Ph.D. diss]. Univ North Carolina Chapel Hill; 2011.

RELATED PAPERS

RELATED TOPICS

Log In

Interphase mass transfer between fluids in subsurface formations: A review

Interphase mass transfer between fluids in subsurface formations: A review

Related Papers

RELATED PAPERS

RELATED TOPICS