Effects of Pore‐Scale Heterogeneity on Macroscopic NAPL Dissolution Efficiency: A Two‐Scale Numerical Simulation Study

Interphase mass transfer or dissolution coefficient of nonaqueous phase liquids (NAPL) is an important parameter in predicting the transport of contaminant species in porous media. While the literature offers valuable insights into the dependence of this coefficient on different parameters at the continuum scale (e.g., contaminant saturation and Darcy velocity), effects of pore‐scale heterogeneity on macroscopic dissolution coefficient have received little attention. In this work a three‐dimensional pore‐scale model is developed to simulate interphase mass transfer over different synthetic pore network structures with various pore radii correlation lengths. The pore network modeling simulates dissolution of immobile NAPL into water (single phase) through diffusive throats for the water‐NAPL interface. The impacts of pore network spatially correlated heterogeneities, NAPL saturation/distribution, and aqueous phase velocity on NAPL mass transfer coefficient and water‐NAPL interfacial surface area are studied. These macroscopic properties are then employed in two‐dimensional continuum‐scale domains formed by concatenating 20 by 20 pore networks in x and y directions. The results highlight the impact of pore‐scale heterogeneity on the distribution of NAPL and subsequently on the dissolution rate (i.e., dissolution coefficient). An uncorrelated distribution of pore radii consistently leads to higher NAPL dissolution coefficient than spatially correlated heterogeneity. The results of continuum modeling show that NAPL dissolution rates are only different between domains formed by correlated and uncorrelated pore networks at very high flow rates and Darcy velocities. However, for typical values of Darcy velocity in groundwater systems, variation in mass transfer coefficient due to pore‐scale heterogeneity is minimal for efficient mass removal.


Introduction
Interphase mass transfer between fluids in porous media occurs in numerous industrial and natural applications (Abriola, 1989). Such processes involve a series of phenomena, including (i) diffusion and convection through the bulk phase of one fluid toward the interface between phases; (ii) accumulation, adsorption/desorption, convection, diffusion, or chemical reaction at the interface; and (iii) diffusion and convection away from the interface into the bulk phase of the second fluid (Giavedoni & Deiber, 1986). Some of these processes may be sufficiently fast that they may be neglected when considering overall transport rates. Others may be rate limiting due to chemical reactions or physical impedances (Powers, 1992). Significant work has been conducted in recent years especially in the field of nonaqueous phase liquid (NAPL) remediation to elucidate the factors that affect this process (e.g., Aydin-Sarikurt et al., 2017;Brusseau et al., 2002;Held & Celia, 2001;Imhoff et al., 1994;Keller & Sirivithiyapakorn, 2000;Kokkinaki et al., 2013;Nambi & Powers, 2003;Pan et al., 2007;Saba & Illangasekare, 2000), and the review papers of Agaoglu et al. (2015) and Essaid et al. (2015) have been dedicated to this topic. The characteristic lengths of the conducted studies range from the pore scale to the mesoscale and further on to the field scale.
Numerous studies have shown that as a result of soil heterogeneity and the nonuniform spatial distribution of the NAPL mass, NAPL bypassing by water can occur at all length scales resulting in complex relations between interphase mass transfer and flow and porous media properties (Agaoglu et al., 2012;Held & Celia, 2001;Maji & Sudicky, 2008;Marble et al., 2008;Unger et al., 1998;Zhou et al., 2000). Among these works, however, pore-scale modeling studies that focus on the impacts of pore-scale heterogeneities on interphase mass transfer are either limited or only focused on random heterogeneity in the lattices. For example, Jia et al. (1999) modeled NAPL dissolution from etched-glass micromodel-based pore networks (regular lattice of 29 × 15 pores of volume 0.5 to 2.0 cm3); Knutson et al. (2001) used lattice Boltzmann modeling to simulate water flow and solute transport from randomly distributed circular NAPL blobs in a twodimensional porous media; Zhao and Ioannidis (2003) modeled NAPL dissolution in homogeneous and heterogeneous lattice networks. More recently, Agaoglu et al. (2016) modeled NAPL dissolution in a lattice pore network of two-dimensional square domains for uniform and heterogeneous media representative of the porous medium used in the experiments of Powers (1992). The cumulative distribution function of pore throat diameters and pore throat lengths was defined with a van-Genuchten-type equation. Bahar et al. (2016) used a two-dimensional (2-D) pore-scale heterogeneous architecture to model bioenhanced NAPL dissolution. Subsurface hydrogeologic heterogeneity generally produces complex NAPL source zones that are different from uniform dispersed zones of residual NAPL blobs and ganglia. Therefore, different forms of heterogeneity can significantly complicate interphase mass transfer processes Essaid et al. (2015).
At the pore scale, intensive research has been conducted to elucidate the impact of the fluid-fluid interfacial area on dissolution mechanism and evaluate the fluid-fluid interfacial area as a function of capillary pressure and saturation (Godinez-Brizuela et al., 2017;Joekar-Niasar et al., 2010;Reeves & Celia, 1996). However, less attention was dedicated to determine the specific effects that pore-scale heterogeneity has on dissolution processes. One naturally occurring form of heterogeneity that exists at the pore scale is that of spatially correlated pore radii and pore throat distribution. Several comprehensive studies have emphasized the importance of accounting for the effects of spatial correlations of macroscopic-level parameters such as permeability or hydraulic conductivity, flow channeling, and connectivity on transport coefficients and solute transport (e.g., Dentz et al., 2002;Knudby & Carrera, 2005;Knudby & Carrera, 2006;Renard & Allard, 2013;Vogel, 2000). However, to our best knowledge, only a limited number of studies have been specifically dedicated to understand the effect of pore-level spatial correlation on macroscopic transport properties (see, e.g., Bernabé & Bruderer, 1998;Bruderer & Bernabé, 2001;Bruderer-Weng et al., 2004;le Borgne et al., 2011;Liu & Mostaghimi, 2017;Makse et al., 2000). More recently, Babaei and Joekar-Niasar (2016) and Liu and Mostaghimi (2017) incorporated pore-scale correlated heterogeneity into pore network models and studied their impacts on macroscopic diffusivity and reaction rates, respectively.

Objectives of This Study
The introduction above specifically highlighted a significant knowledge gap related to the impacts of porescale correlated heterogeneity on dissolution rate and macroscopic interphase mass transfer coefficient. Therefore, this study will present the results of single-phase NAPL dissolution predictions from spatially correlated and uncorrelated pore network models and will present the impacts of pore-level correlated heterogeneity on dissolution rate and interphase mass transfer. Two main macroscopic properties that are derived from pore-scale modeling, effective permeability and interphase mass transfer coefficient, will be employed in a single-phase continuum-scale framework where each grid block is represented by one pore network. Our hypothesis is that the combination of effects on both properties may introduce changes in NAPL dissolution rate and remediation efficiency.

NAPL Dissolution Modeling
In this study, NAPL dissolution is numerically simulated at the pore scale within three-dimensional unstructured pore networks. The macroscopic upscaled properties from pore network simulations are applied in a two-dimensional continuum-scale model. While the governing equations for flow and NAPL dissolution at the continuum scale are described in the Supporting Information (Part A), here we describe pore-scale governing equations and our simulation assumptions. The simulation code uses a commonly applied simplification of NAPL dissolution in water based on the work of Jia et al. (1999) that was adapted later by others (e.g., Agaoglu et al., 2016;Aydin-Sarikurt et al., 2017;Held & Celia, 2001). No aqueous flow is assumed in the pores with NAPL, and NAPL is assumed as an immobile phase. A steady state mass balance equation for each of the pores is used to calculate the pressure distribution and aqueous fluid flow over the pore network: where Q ij is the flow rate of aqueous phase between pore i and j and N i is the number of pores that are connected to pore i. The flow rate within the pore throat that connects pore i to j, Q ij ,is calculated using Hagen-Poiseuille equation: where r ij is the radius of the throat, which connects pore i to j, l ij is the length of that throat, μ is the aqueous phase viscosity, and ΔP ij is pressure difference between two sides of throat. It is assumed that there is no pressure drop in the pore bodies. Substituting equation (2) into equation (1) gives a mass balance equation for each of the pores as a function of pressure: Applying equation (3) for all pores results in a linear system of equation for pressure of pores. Pressures are defined for the inlet and outlet face of the pore network to solve the linear system of equation and to obtain the pressure distribution over the pore network. By applying the Hagen-Poiseuille expression of equation (2), the flow rates within the throats are calculated.
After simulation of steady state aqueous phase flow, the dissolved NAPL transport in the aqueous phase is simulated using the advection-diffusion equation: where C tþ1 i is the concentration of pore i at the next time step; V i is the volume of pore i; C *t ij is the concentration of pore j, which is connected to pore i if the flow direction is from pore j to i and it is the concentration of pore i if the direction is from pore i to j; ΔC t ij is the difference of concentration between pore i and j; A ij is the throat cross-section area between pore i and j; and D is the diffusion coefficient of NAPL in the aqueous phase. The third term on the right-hand side of equation (4) is applied when pore i is connected to pore j, which is occupied by NAPL. NAPL dissolution occurs only at the throats that connect a pore filled with the aqueous phase with a pore filled with the NAPL phase. The concentration of NAPL at the interface of water and NAPL is assumed to be the solubility of NAPL in the aqueous phase (C s ); that is, it applies the thermodynamic driving force for the mass transfer between aqueous phase and NAPL phase. The fourth term of equation (4) is only applied for those pores at the outlet, as a sink term, and is defined as flow rate for those pores at the outlet face of the pore network. For the upstream boundary condition, it is assumed that no NAPL enters the pore network.
The term ∑ Water Resources Research species from its own phase, along the boundary layer (in our case diffusive throats) and into the aqueous phase. This simplification has been used elsewhere Held & Celia, 2001;Jia et al., 1999). The three steps of interphase mass transfer are schematically shown in Figure 1. In the formulation above, only the third step is considered and it is assumed that the other two steps are not rate limited. The interface of NAPL-water phase forms in the stagnant throats where only diffusion takes place. We refer to these throats as stagnant throats because they have a zero flow rate and are not part of the advective regime whereby water flow occurs. Also, such stagnant throats are not located inside the NAPL clusters where neither diffusion nor advection takes place. Recent investigations have shed light into the importance of stagnant regions in two-phase porous media (for example, NAPL-water systems) on hydrodynamics of solute transport (e.g., Aziz et al., 2018;Hasan et al., 2019;Karadimitriou et al., 2016); therefore, we perform an analysis of the number of stagnant throats for different networks to discover possible links between this parameter and the effective dissolution coefficient of NAPL.
A steady state condition of solute transport is required to calculate the upscaled NAPL mass transfer rate coefficient (Powers et al., 1994). There are two approaches to achieve the steady state conditions for solute transport in the pore network: (1) solve solute transport model (equation (4)) explicitly in time until a fully developed solute transport is reached (solute concentration of each of pores does not change over the time) and (2) solve the steady state form of the solute transport model (equation (4)); that is, the left-hand side of equation (4) is set to zero. Same distributions of NAPL concentration through the pore network are achieved using both approaches in the code developed. Figure 1. Schematic of NAPL dissolution into water from a pore filled with NAPL phase to a pore filled with the aqueous phase through a diffusive stagnant water-filled throat. Three steps of dissolution include (1) diffusion of NAPL species in its own phase: D n ∇C n NAPL where D n is the self-diffusivity of NAPL, (2) partitioning through thermodynamic equilibrium: where K p is the partitioning coefficient and C int NAPL is the NAPL concentration at the interface, and (3) diffusion of NAPL species in the aqueous phase: D w ∇C w NAPL where D w = D is the diffusion coefficient of species in the aqueous phase.

Water Resources Research
Basically, there are two approaches for calculation of NAPL mass transfer coefficient. In the practical approach for field application or in the lab scale (e.g., in a one-dimension column experiment), it is difficult to measure the average residence concentration of NAPL. Therefore, NAPL mass transfer can be calculated using effluent concentration of NAPL based on the steady state transport equation (Powers et al., 1994): where q is the aqueous velocity, L is the length of the column, a is the specific interfacial area, and C eff is the NAPL effluent concentration. However, a theoretically sounder approach-that reflects more dynamics of the system, flow conditions, and stagnant zones' effects for NAPL dissolution analysis-is based on the residence concentration of NAPL in an REV (Representative Elementary Volume) scale. We use this approach in our work. Dissolution rate coefficient of NAPL as an upscaled mass transfer coefficient is calculated based on the simulated NAPL concentration of pores in the pore network model (Agaoglu et al., 2016): where A is the overall interfacial area (interphase mass transfer area) in the pore network, that is, the summation of cross-sectional areas of the throats, which connect NAPL-filled pores with water-filled pores. The concentration, C, is defined as the average resident concentration of NAPL in the aqueous phase over the entire pore network model. The scale of domain is representative of pore-scale heterogeneity and correlation lengths chosen so that the average resident concentration is representative . The NAPL dissolution rate into the aqueous phase, ∂M/∂t, is equal to the rate of NAPL leaving the pore network at steady state condition; that is, it is calculated by multiplying the NAPL concentration of pores at the outlet with its flow rate.
According to the terminology commonly used in soil remediation, K f is also referred to as explicit mass transfer coefficient as opposed to lumped mass transfer coefficient ( . This terminology also leads to the definition of explicit Sherwood number and modified Sherwood number as Sh = K f d m /D and Sh′ ¼ K L d 2 m =D, respectively, where d m is the mean grain diameter. The simulation code is first validated against the experimental results of Aydin-Sarikurt et al. (2017) for single-phase constant-flux pool dissolution of immobile 1,2-dichlorobenzene NAPL. The rationale to use these experimental data is that chlorinated solvents are major sources of groundwater contamination and they have been used in numerous studies varying from pore to field scale, and thus, understanding such processes may help for remediation effectiveness. Whereas the experiment is about pool dissolution, since similar formulation as above has been used by the authors to validate their code, we also use the experiment for validation of our code.
Similar to the original paper, a two-dimensional lattice pore network with cubical pore bodies and properties summarized in Aydin-Sarikurt et al. (2017; Table 1) is set up. The authors (Aydin-Sarikurt et al., 2017) conducted glass bead experiments, which can be represented sufficiently accurately by the formulation presented above. In fact, given the uniform porous media used in the experiments, the pore chambers can be assumed to have a cubical shape, which is commonly used in pore network models . The domain is set to the size of an experimental tank with length of 100 mm, height of 30 mm, filled with uniform porous media with glass bead diameters of L c = 0.7 mm, square cross-section area of A= 0.25 mm 2 , and length of l= 0.3 mm. The inflow concentration is assumed to be equal to zero. Along the bottom boundary, NAPL dissolution (∑ In the analysis of results in section 3 we use Péclet number defined as Pe = Lu/D, where L is the characteristic length (here, it is the average throat length of the pore network), u is the interstitial velocity of the network, and D is diffusion coefficient. The Péclet number is a dimensionless number to describe the ratio of advective transport rate to diffusive transport rate and as such can provide an indication of transport regime and diffusion versus advection dominance. Commonly, mass transfer coefficient or Sherwood number is represented against Péclet number.

Simulation Parameters 2.2.1. Properties of Pore Networks
First, we generate a series of pore networks based on the algorithm described in Babaei and Joekar-Niasar (2016) and further information provided in the Supporting Information (Part B). Table 1 summarizes the geometric properties and correlation lengths of 14 networks generated by the algorithm used in this study. Two geometric network sizes of 4.6 mm (=L y ) × 4.6 mm (=L z ) × 7 mm (=L x ), referred to as PN-1, and 4.6 mm (=L y ) × 4.6 mm (=L z ) × 14 mm (=L x ), referred to as PN-2, are used to investigate the effect of physical length scale on dissolution. In terms of pore radius distribution, the networks have μ R = 16.65 μm, σ R = 9.57 μm, and R ≤ 50 μm, where μ R is the mean of pore radii, σ R is the standard deviation of pore radii, and R is the pore radii. The mean and variance of R (μ R and σ R ) are related to the mean and variance of the . The minimum distance between each pair of points (body centres) is set to 100 μm, which allows for a maximum pore radius of 50 μm. The average coordination number of the generated network is 5.5. After randomly populating the space with pores, the mean of the distances between the centers of the pore bodies is on average equal to 119 μm. This distance is often used as the characteristic length of the pore network models (Bijeljic et al., 2004). It is in good agreement with the Berea sandstone with a characteristic length of 131 μm and also is with the measurement range of 100 to 150 μm determined by other modeling studies (Bijeljic et al., 2004;Mostaghimi et al., 2012;Øren & Bakke, 2003). In total, as reported in Table 1, we assign the pore centers of PN-1 and PN-2, such that we have two uncorrelated pore networks and 12 correlated pore networks with different correlation lengths. The correlation lengths are either isotropic or anisotropic in the x direction, the y direction, or the diagonal direction. The pore throat size distributions for the network sizes where r th is the radius of throat and r i and r j are the radii of two neighbouring pore bodies.

NAPL Initial Distribution
Tetrachloroethylene (also referred to as Perchloroethene, Perchloroethylene, or PCE) is selected as the entrapped NAPL phase in the networks with a diffusion coefficient of 8.2 × 10−6 cm2/s and solubility of 0.0012 mole/L in water. The residual saturation and distribution of NAPL in the pore networks are functions of fluid properties, pore geometry (topology and pore/throat size distribution of pore networks), and the wettability of pore network (Mercer & Cohen, 1990). Here, it is assumed that the NAPL is trapped in the large pores. This can be justified based on the fact that NAPL generally occupies large pore spaces in a water-wet medium (Powers, 1992), and most of the studies have shown that nonwetting phase tends to be trapped in the large pores. For example, in a set of two-dimensional micromodel experiments, Chatzis et al. (1983) and Wardlaw (1982) showed that NAPL preferentially distributes in large clusters of pores with a high pore body to throat size ratio. Additionally, a direct visualization method using magnetic resonance imaging in a three-dimensional packed system with a heterogeneous distribution of sand fraction showed that most NAPL is trapped in the coarsest sand as pools and not as residual ganglia (Zhang & Guo, 2008). Furthermore, as described in the Supporting Information (Part C), we conducted a simple pore network modeling of drainage and imbibition to confirm that the large pores are occupied by NAPL.

Continuum-Scale Properties
Four hundred of PN-2 networks (4.6 × 4.6 × 14 mm 3 ) are put next to each other in a 20 (=N y ) × 1 (=N z ) × 20 (=N x ) mesh forming a 9.2 cm (=L Y ) × 0.46 cm (=L Z ) × 28 cm (=L X ) mesoscopic domain, where each network is treated as one continuum scale grid block. We do this PN-2 L0, iso. L2, aniso-x.L2, aniso-y.L2, and aniso-d. L2. Therefore, five continuum models are formed with homogeneous permeability equal to that of each constituent pore network listed in Table 1. The continuum domains are assumed to be initially 75% saturated with NAPL.
Next, the macroscopic properties including effective permeability as a function of NAPL saturation, NAPL dissolution coefficient as a function of NAPL saturation and flow rate, and specific water-NAPL interfacial area, are used to model single phase NAPL dissolution at the continuum scale using the governing equations presented in the Supporting Information (Part A). Pore network modeling provides information on NAPL mass transfer coefficient, specific water-NAPL interfacial area, and effective permeability at different NAPL saturations and aqueous velocities for different constitutive pore networks. Therefore, pore-scale physics associated with different NAPL distribution at various NAPL saturations (represented in the pore-scale model) is incorporated into the continuum-scale simulation. That is for each saturation, an interpolation is made between pore networks used with S n = 19%, S n = 38%, S n = 55%, S n = 75%, and S n = 85%. Other properties used in the continuum scale are listed in Table 2. A constant water flow rate in the x direction from left-hand side is prescribed, and the extraction of effluent concentration takes place from right-hand side boundary. The initial concentration of NAPL at aqueous phase is assumed to be the solubility of the NAPL species.
Three different flow rates are tested (q= 0.7, 2.2, and 11 cm 3 /min) to evaluate effects of flushing rate (or residence time) on dissolution rate. Given that the porosity from our pore networks is too low (ϕ = 2.5%), the flow rates chosen lead to very high interstitial Darcy velocities (95, 300, and 1,500 m/day) that can only occur in the vicinity of the wellbore. However, at such very large velocities, mobilization of NAPL can occur and the dissolution can no longer be stated as the linear driving force model presented in the Supporting  (Abriola et al., 2004), occurring again only in the vicinity of injection or extraction wells and under two-phase conditions. The reason to choose broadly unrealistically high flow rates is to illustrate and highlight the effects of pore-scale heterogeneity (if any) at the continuum scale. This will be further detailed in the next section.

Effective Permeability, Interfacial Area, and Diffusive Throats
First, we calculate the aqueous phase effective permeability of the pore networks at different NAPL saturations, as well as interfacial areas between NAPL and aqueous phase (A) in the pore networks. The results are shown in Figure 3. While absolute permeability (as listed in Table 1) reduces as correlation length increases, effective permeability follows the same trend except for the pore networks correlated in the y direction and diagonally at higher NAPL saturation. For fixed NAPL saturation (S n ), the effective permeabilities of the uncorrelated pore networks (L0) are lower than the correlated ones (isotropic and anisotropic in the x direction), clearly pointing out the importance of correlated pore-scale heterogeneity, which exhibits higher permeability. Highest effective permeability for fixed S n is obtained for the anisotropic pore networks in the direction of flow (aniso-x.L2). Increase in S n decreases the effective permeability as the connectivity of pore networks is reduced by immobile NAPL phase. Similarly, the impact of correlation lengths is reduced for higher S n due to a reduction in pore connectivity.
At lower S n , the water-NAPL interfacial area in the uncorrelated pore networks (L0) is lower than the correlated pore networks. This is because at lower S n , NAPL is not distributed as clusters in the correlated pore networks, and the large pores, which are filled with NAPL, are neighbors through correlation. Therefore, there is a greater chance that a large pore filled with NAPL is connected to a large pore filled with water in the correlated network. In contrast, a large pore filled with NAPL is more likely to be connected to a random small pore in the uncorrelated network. At higher S n , the water-NAPL interfacial area in the uncorrelated pore networks (L0) is higher than correlated pore networks, because NAPL is distributed as clusters in the correlated pore networks and this results in lower interfacial area between NAPL and aqueous phase.
The subsequent analysis is conducted to calculate the number of diffusive/stagnant throats in the domains (throats similar to the one shown in Figure 1). As shown in Figure 4a, the number of stagnant pore throats increases as NAPL saturation increases. That is, as NAPL saturation increases, the throat discontinuity and therefore number of stagnant pore throat increase. However, interestingly, the differences in number of stagnant throats between the uncorrelated and correlated networks also increase as NAPL saturation increases, so that the uncorrelated network has the highest number and the isotropic. L2 (correlated) has the lowest number of stagnant throats (as shown in Figures 4b and 4c). The reason for higher number of stagnant throats in the uncorrelated compared to the correlated pore networks is that while NAPL is distributed as clusters in the correlated pore networks (Figure 4c), NAPL is distributed through the entire uncorrelated pore networks (Figure 4b). This results in higher discontinuity and stagnant pore throats.

Mass Transfer Coefficient
K f and Sh as functions of aqueous velocity and Péclet number for all pore networks for various NAPL saturations are presented in Figure 5 for S n = 38%, Figure 6 for S n = 55%, Figure 7 for S n = 75%, and Figure 8 for S n = 85% (the results for S n = 19% are very similar to S n = 38% and therefore are not shown).
When comparing NAPL mass transfer coefficients for different correlation lengths, the results are similar for both PN-1 and PN-2 ( Figures 5-8). For example, the highest mass transfer coefficient is observed in the uncorrelated pore networks for both PN-1 and PN-2. Moreover, the range of mass transfer coefficient for a specific correlation length is very similar for both PN-1 and PN-2. For example, the mass transfer coefficient ranged between 2.6 × 10 −6 to 6.9 × 10 −6 m/s for uncorrelated PN-1 and 2.9 × 10 −6 to 6.8 × 10 −6 m/s for uncorrelated PN-2. This result indicates that the mass transfer coefficient is independent of the scale and the length of pore network in this study.

Water Resources Research
The results also demonstrate that as expected, the aqueous velocity plays a major factor in the NAPL mass transfer coefficient. The rate of dissolved NAPL removal from the pore network domain increases with the increase of the aqueous pore-scale velocity. As aqueous velocity increases, the average concentration of NAPL in the aqueous phase decreases; hence, larger driving force (C s − C) between NAPL and aqueous phase remains. Consequently, increasing NAPL dissolution flux leads to an increase in macroscopic NAPL dissolution coefficient. The NAPL mass transfer coefficient does not change significantly at a very high aqueous velocity and increases to a maximum as velocity is increased.
While the minimum values of the NAPL mass transfer coefficient for different pore network structures are different from each other, the maximum values are close to each other as shown in Figure 5. The maximum value of NAPL mass transfer coefficient depends on pore network structure, NAPL distribution, and molecular diffusion coefficient of NAPL component in the aqueous phase. We propose that the maximum NAPL mass transfer coefficient can be directly calculated from D/ < l>, in which <l> is the average length for those throats that connect pores filled with NAPL to the pores filled with the aqueous phase. For example, the average throat length (connecting of aqueous and NAPL pores) for the 4.6 × 4.6 × 7.aniso-L1 pore networks with 38% NAPL saturation is 119.7 μm; this results in D/ < l > = 6.85 × 10 −6 m/s, which is very close the  Among the various pore networks used in this work, uncorrelated pore networks (L0) result in the highest NAPL mass transfer coefficients and the isotropically correlated pore networks (iso. L2) have the lowest NAPL mass transfer coefficient. These trends are consistent for various NAPL saturations. To discuss these observations, first we focus on high NAPL saturation values. For these cases, both water-NAPL interfacial area and number of stagnant throats are the highest for L0 and lowest for iso. L2; therefore, the results for NAPL mass transfer coefficient at high NAPL saturations are justified and are as expected. It is more likely to distribute and dissolve NAPL through the entire pore network with the higher water-NAPL interfacial area and more stagnant throats as they are in the uncorrelated pore networks, which result in an increase in NAPL dissolution rate and coefficient. Also as shown in Figure 3, the water-NAPL interfacial area decreases in isotropic pore networks especially at higher NAPL saturation (indicating that the NAPL is not well distributed but rather more consistent to that of a pooled distribution), which results in a lower NAPL dissolution rate.
For the lower NAPL saturation simulations, although the water-NAPL interfacial area for uncorrelated pore network is lower than the correlated pore networks and the number of stagnant throats are equal between various networks at low NAPL saturation, NAPL mass transfer coefficients of L0 are still mostly higher compared to those for the correlated pore networks. While neither NAPL-water interfacial area (NAPL mass transfer coefficients are nearly independent of specific NAPL-water interfacial area; Cho et al., 2005) nor the number of stagnant throats are directly responsible for this observation (Mobile et al., 2016), it seems that

Water Resources Research
the pore-scale correlation leading to distribution of NAPL is the main factor for this observation. That is, the well-distributed NAPL in the uncorrelated pore networks overcomes the lower interfacial area leading to the higher K f . That is, while the black points (NAPL-filled pores) are well distributed through the entire uncorrelated pore network (Figure 4b), they are distributed as clusters in isotropic pore network (Figure 4c). This positive distribution cannot be attributed to water-NAPL interfacial area or number of stagnant throats. This finding clearly highlights the importance of considering the pore-scale correlation lengths for low NAPL saturations (which are more likely to occur than high NAPL saturations). Also in Figure 9, the average K f over all velocities are reported for different networks and NAPL saturations (including S n = 19%) for better comparison. For instance, quantitative analysis demonstrates that at the highest S n saturation of 85%, for a dimensionless correlation length of 0.1, a reduction of mass transfer coefficient from 5 × 10 −6 to 3 × 10 −6 m/s is observed.
As shown in Figures 10 and 11, while NAPL is well distributed through the uncorrelated pore network (black points), NAPL is poorly distributed as clusters in the correlated pore networks. The well-distributed NAPL results in its dissolution through the entire pore network, which increases the NAPL removal rate from the pore network and then results in an increase of K f . As shown in Figure 11i, NAPL is not well distributed (trapped as a cluster) for the isotropic high NAPL saturation network, and average residence concentration of dissolved NAPL in the aqueous phase is less than the other pore networks. These reasons cause a reduction in NAPL mass transfer coefficient compared to other pore networks.
Furthermore, the results also indicate that the orientation of the NAPL clusters has a consistently slightly positive impact on NAPL mass transfer coefficient, especially at higher NAPL saturations. The orientation of the NAPL clusters relative to the flow direction has an impact on NAPL bypassing by the aqueous phase, surface contact between NAPL and fresh water, and average effluent concentration (Agaoglu et al., 2016). Among the anisotropic pore networks, although there is a similar water-NAPL interfacial area for certain correlation lengths (as shown in Figure 3), those pore networks with a diagonal correlation length show better performance in terms of NAPL dissolution and dissolved phase removal. This is because that diagonal correlation length leads to less NAPL bypassing by aqueous phase, which results in higher dissolution rate and effluent concentration. For example, for the 85% NAPL saturation scenario, the average concentration of NAPL through the diagonal-anisotropic pore network (5.83 × 10 −4 mole/L, Figure 11g) was greater than those of anisotropic in the x and y directions (which are 5.53 × 10 −4 and 5.44 × 10 −4 mole/L, respectively, Figures 11e and 11c), indicating higher NAPL dissolution rate and mass transfer coefficient in diagonal pore networks. Generally, the diagonal correlation length orientations in the pore networks help dissolve the NAPL through the entire pore network in both x and y directions.

Continuum-Scale Modeling
NAPL recovery as a function of PV flushed is simulated for five continuum models (described in section 2.2.3) under various pore network regimes with explicit representation of effective permeability, water-NAPL interfacial area (Figure 3), and NAPL dissolution coefficient (Figures 5-7). Figures 12a Results indicate that the effect of pore-scale correlation length (and its impact on NAPL dissolution coefficient) gradually increases as water flushing rate increases. For example, although the largest difference in NAPL dissolution coefficient is observed between the uncorrelated pore network (L0) and isotropic (iso. L2) pore network (as shown in Figure 7), their NAPL recoveries are quite similar at lower flushing rate of 0.7 cm 3 /min. The reason for this behavior is that the residence time for the aqueous phase increases at low flushing rate. For example, the residence times of the aqueous phase in one grid block in these continuum models-defined as the ratio of the PV of grid block over aqueous flow rate in that grid block-are

10.1029/2019WR026035
Water Resources Research 0.21, 0.07, and 0.013 min for flushing rates of 0.7, 2.2, and 11 cm 3 /min, respectively. As the residence time increases, regardless of NAPL dissolution coefficient, there is sufficient time (i.e., mass transfer time) for NAPL dissolution so that NAPL concentration reaches its solubility limit. Therefore, NAPL recovery is

Water Resources Research
not noticeably affected by dissolution coefficient or water-NAPL interfacial area at low flushing rates. Twodimensional NAPL saturation maps inset in Figure 12 indicates that NAPL removal from the domain at the low flushing rate is more piston-like than the scenarios at the high flushing rates. This means that a

Water Resources Research
sweeping dissolution regime occurs at low flushing rates (i.e., we observe more effective uniformly swept zone throughout the entire system, which allows for steady state equilibrium dissolution to occur). Therefore, if the flushing rate is sufficiently low (high residence time), NAPL phase is dissolved to its solubility limit in the aqueous phase.
In contrast to the low flushing rate, pore-scale correlation length of the porous network affects NAPL removal under higher flushing rate conditions. For the higher flushing rate scenario (11 cm 3 /min), NAPL recovery is highest for the uncorrelated porous media (L0) and lowest for the isotropic case at early stages of the simulation. These results are consistent with pore network modeling results where Figure 12. NAPL recovery as a function of aqueous pore volume (PV) flushed with three different flushing rates of (a) 0.7 and 2.2 cm 3 /min and (b) and 11 cm 3 /min for five continuum models described in section 2.2.3. Insets: twodimensional NAPL saturation profiles at specific PVs injected and average NAPL dissolution rate as a function of PV flushed for the case of (a) 2.2-and (b) 11-cm 3 /min flushing rates.
the highest NAPL mass transfer coefficient ( Figure 9) and water-NAPL interfacial areas ( Figure 3) were observed for the uncorrelated pore networks (L0) under higher NAPL saturation conditions. This results in higher NAPL dissolution rates for uncorrelated porous media (L0) at early stages of water flushing as shown in the inset of Figure 12b for the average NAPL dissolution rate as a function of PV flushed for 11-cm 3 /min flushing. However, NAPL dissolution rate decreases as NAPL is dissolved and NAPL saturation decreases. This reduction of NAPL dissolution rate is due to the decrease of NAPL dissolution coefficient and water-NAPL interfacial area at lower NAPL saturations. The reduction of NAPL dissolution rate for uncorrelated porous media (L0) is more than correlated structures. This is because while the NAPL dissolution coefficients of all pore networks are almost similar at a high flow rate (as shown in Figure 5), water-NAPL interfacial area for uncorrelated pore network (L0) is lower than correlated pore networks at lower NAPL saturation (as shown in Figure 3). Therefore, the time taken for complete NAPL removal from correlated pore networks is less for the uncorrelated pore network. Thus, the efficiency of NAPL dissolution and removal at lower NAPL saturation (after a long period of water flushing) for continuum domains formed by correlated pore networks is higher than that of formed by uncorrelated pore networks.
We chose very high flow rates that led to large interstitial velocities. These flow rates highlighted the trend of difference in NAPL recovery between continuum-scale domains formed by correlated and uncorrelated networks, albeit under unrealistic settings. For example, the same trend of difference observed for 11 cm 3 /cm is observed for 2.2 cm 3 /cm (in smaller magnitudes). Nevertheless, as discussed before, the Darcy velocity in field-scale groundwater and remediation processes is significantly lower than the ones tested in this work. Therefore, we expect negligible effects of pore-scale heterogeneity for the dissolution rate of NAPL at the continuum scale. Of course, this conclusion is also true for any other literature study that has obtained an order of magnitude variation in mass transfer coefficient rate of NAPL at pore scale. In other words, since we used a continuum model, we could evaluate NAPL recovery for the range of variations for this parameter.

Conclusions and Future Work
The results of this research can be summarized as follows: 1. As the pore-scale heterogeneity shifts from a random uncorrelated distribution to correlated distribution, the initial distribution of NAPL, water-NAPL interfacial area, effective permeability, and interphase mass transfer changed as functions of NAPL saturation and Darcy velocity (i.e., water flow rate). 2. NAPL mass transfer coefficient was the highest for uncorrelated pore networks as NAPL was well distributed in this structure. This was observed for different initial NAPL saturations, and the trend could only be attributed to the nature of heterogeneity rather than water-NAPL interfacial area and number of stagnant/diffusive throats. 3. NAPL mass transfer coefficient increased with water velocity and reached similar maximum values regardless of pore network structure (correlated and uncorrelated). The maximum NAPL mass transfer coefficient is dependent upon the average throat length (the throats that connect NAPL-filled pores with water-filled pores) and the diffusion coefficient of NAPL in the aqueous phase through the D/ < l> relationship. 4. Using macroscopic properties obtained from pore network modeling for two-dimensional homogeneous continuum-scale domains, it was determined that although NAPL dissolution coefficients and water-NAPL interfacial area were different in various pore networks at a flushing rate of 0.7 cm 3 /min, NAPL dissolution/removal rates were the same. Such behavior is attributed to the fact that the residence time was sufficiently high for NAPL dissolution at typically low flow rates and Darcy velocities occurring in groundwater regardless of the porous media structure. 5. As the flushing rate and Darcy velocity increased well beyond typical magnitudes occurring in groundwater (e.g., in the vicinity of well), NAPL removal rates became different for correlated and uncorrelated porous media structures, whereby the continuum models constructed with the uncorrelated pore networks had lower NAPL dissolution and removal than those formed by the correlated networks. The difference is such that the NAPL recovery was higher for continuum domains constituted by correlated pore networks. The reason for this behavior is that water-NAPL interfacial area for uncorrelated pore 10.1029/2019WR026035

Water Resources Research
network is lower than correlated pore networks at lower NAPL saturation, while their NAPL dissolution coefficients are almost similar at high flow rates.
While we studied the interphase mass transfer coefficient in the context of NAPL dissolution in groundwater, this parameter is also an important for other applications of subsurface engineering such as geological CO 2 storage (Babaei & Islam, 2018;Emami-Meybodi et al., 2015), acidification processes in petroleum engineering Panga et al., 2005), and miscible displacement in porous media (Bretz & Orr, 1987). Therefore, the findings of this research are applicable to these areas of research as well.
Recommendations for future work consist of expanding the methodology for SEAR (surfactant-enhanced aquifer remediation), in situ chemical oxidation, and other remediation applications in which this two-scale numerical modeling approach can be used to determine mesoscopic properties via pore-scale modeling. One important and often missed research challenge for NAPL remediation is the consideration of concurrent dissolution and mobilization of NAPL that can be systematically investigated for such heterogeneous pore-scale models. Some experimental examples include the works of Detwiler et al. (2001), Ghosh and Tick (2013), and Tsakiroglou et al. (2013). Pore-scale heterogeneity may also play a role on mobilization extent given that the Darcy velocity involved in mobilization is higher than dissolution processes.