Upscaling dissolution and remobilization of NAPL in surfactant-enhanced aquifer remediation from microscopic scale simulations

The dissolution and mobilization of non-aqueous phase liquids (NAPL) blobs in Surfactant-Enhanced Aquifer Remediation (SEAR) processes are upscaled using dynamic pore network modelling of three-dimensional and unstructured networks. We considered corner ﬂow and micro-ﬂow mechanisms including snap-oﬀ and piston-like movement for two-phase ﬂow. Moreover, NAPL entrapment and remobilization were evaluated using force analysis to develop capillary desaturation curve (CDC) and predict the onset of remobilization and complete removal of entrapped NAPL blobs. The corner diﬀusion mechanism was also applied in the modeling of interphase mass transfer to represent NAPL dissolution as the dominant mass transfer process. Our model showed that although surfactants enhance NAPL recovery during two-phase ﬂow, surfactant-enhanced remediation of


Introduction
Non-aqueous phase liquids (NAPLs) are liquid solution contaminants, such as chlorinated solvents and petroleum hydrocarbons, which can be long-term threats to subsurface water resources because of relatively low water solubility and high toxicity (Farthing et al., 2012;Mayer & Hassanizadeh, 2005).NAPLs migrate to the subsurface and can become trapped as disconnected blobs, ganglia, and pools due to gravitational, viscous, and capillary forces ( Abriola, 1984).The entrapped NAPLs in the water saturated zone could be removed during displacement, dissolution, and remobilization processes.Surfactant-enhanced aquifer remediation (SEAR) is one of the most effective techniques that improves remediation of entrapped contaminants by increasing NAPL solubility in water and reduction of interfacial tension (IFT) (Baran et al., 1994;Huo et al., 2020;Javanbakht et al., 2017;Jeong & Corapcioglu, 2003;McGuire et al., 2006;Ramsburg et al., 2005;Roote, 1997).The presence of surfactants in the aqueous phase can change interfacial area which could have a significant effect on the dissolution process, i.e., the accumulation of surfactant monomers and micelles near the interface can retard interphase mass transfer processes (Aydin-Sarikurt et al., 2016;Babaei & Copty, 2019;Grimberg et al., 1995;Javanbakht & Goual, 2016;Lee et al., 1998;Ramezanzadeh et al., 2019).Hence, a practical design of remediation strategies needs reliable numerical simulations in which the complex behaviors of NAPLs and surfactants, as well as mechanisms of multiphase flow and transport should be considered.
As natural porous media are intrinsically heterogeneous, pore-scale heterogeneities can result in variabilities in flow and transport properties at various length scales in a porous medium (Li et al., 2006).Indeed, transport modeling at large scales is not capable of taking into account porescale heterogeneities.Hence, pore-scale modeling and upscaling of subsurface processes are important to predict flow and transport properties with high accuracy in heterogeneous porous media (Farthing et al., 2012;Li et al., 2006;Raoof et al., 2010).As pore network modelling (PNM) has substantial saving in computational cost compared with direct methods (Hammond & Unsal, 2012), this numerical approach was used to simulate SEAR processes at pore-scale in this work.PNM is a conventional method, which represents porous media by a network of interconnected pores and throats with idealized geometries.This approach was proposed for the first time by Fatt (1956) and has been applied extensively by numerous researchers for studying a wide variety of physical and chemical processes in porous media (Blunt et al., 2002;Blunt & Scher, 1995;Dias & Payatakes, 1986;Joekar-Niasar et al., 2007;Yin et al., 2019).
Comprehensive reviews of pore-network models can be found in Blunt et al. (2013) and Joekar-Niasar and Hassanizadeh (2012).Pore-network can be classified into quasi-static and dynamic network models.The quasi-static approach has been successfully applied when capillary force is dominant (Piri & Blunt, 2005), so this approach cannot be used to study cases in which capillarydominated assumptions do not apply.Viscous forces play significant roles in various subsurface flow processes, such as SEAR and EOR methods of surfactant flooding, and high velocity flow regimes happened in naturally and hydraulically induced fractures as well as near well-bore areas (Aghaei & Piri, 2015).In these cases, the combination of capillary and viscous forces determines the flow and transport behavior in porous media.However, in dynamic pore-network approach, viscous forces are taken into account (Khasi et al., 2019).The Dynamic approach can be categorized into two types of algorithms.The first is single-pressure algorithm, which is computationally efficient.Aker et al. (1998) developed a two-dimensional network and investigated the dynamics of drainage process.They calculated pressure using solving mass conservation at nodes, which were assumed to be occupied by only one fluid.Although main terminal fluid-fluid interfaces were modeled, corner flow and snap-off mechanism were not considered.Al-Gharbi and Blunt (2005) implemented a dynamic pore-network model with single-pressure algorithm for modeling of two-phase drainage.They considered wetting layer flow, meniscus oscillation, and snap-off displacement mechanism in their network models.The other algorithm of dynamic pore network modeling is two-pressure algorithm, in which pressure distribution, saturations, and mass flux can be determined for both fluids occupying pore spaces.This algorithm was suggested for the first time by Thompson (2002) for studying imbibition process in fibrous materials.Later, this algorithm was developed and implemented in structured network with coordination number of six for studying two-phase flow in porous media (Joekar-Niasar et al., 2010;Joekar-Niasar & Hassanizadeh, 2011).
Interphase mass transfer can be simulated and upscaled using pore network modelling of twophase flow that occurs when NAPL ceases to be displaced (with an initial residual saturation).
Retraction of NAPL could occur through some processes such as dissolution and volatilization, but as existing phases in saturated media are water and NAPL, the principal mass transfer mechanism is dissolution.A common approach for evaluating NAPL dissolution in pore networks is mechanistic modeling.Zhou et al. (2000) developed mechanistic modeling based on simplified geometric properties and the physical mechanisms controlling NAPL dissolution at pore-scale including pore diffusion, corner diffusion, and mixing streams depending on Péclet number ().Pore diffusion mechanism is defined as the diffusion of solute through NAPLaqueous phase interfaces in the opposite direction of aqueous phase flow, however, corner diffusion represents the diffusion of solute from NAPL-aqueous phase interface into the aqueous phase surrounding the NAPL.The mixing of streams of aqueous phase is a type of dissolution mechanism in which streams of aqueous phase with different dissolved NAPL concentrations are mixed together from different pores (Zhou et al., 2000).Dillard and Blunt (2000) applied the corner diffusion process for NAPL dissolution in a water-wet pore network.The corner diffusion model did not account for a velocity distribution in the aqueous phase.It was also appropriate only for short contact times as infinite boundary condition was employed for concentration.Sahloul et al. (2002) proposed the corner transport model in which the corner diffusion modeling was modified by considering velocity profiles for water flowing around entrapped NAPL ganglia.Zhao and Ioannidis (2003) applied the corner transport model to simulate dissolution process in a NAPL-wet pore network.The latter study used a structured two-dimensional network that was not able to represent pore structure of real porous media.Besides, none of the mentioned studies can model the enhancements by surfactant during dissolution process.
In SEAR and surfactant enhanced oil recovery (EOR), flow and transport processes in porous media has been investigated extensively using numerical techniques at multiple scales.
Numerous researchers developed large-scale continuum models to simulate surfactant flushing in the subsurface (Abriola et al., 1993;Delshad et al., 1996;Fujioka, 2013;Muradoglu & Tryggvason, 2014;Qin et al., 2007;Rathfelder et al., 2003;Sulaiman & Soo Lee, 2012;Zhang et al., 2018;Zhang et al., 2014).Pore-scale modeling is often used to understand fundamental phenomena as well as flow and transport mechanisms at pore-scale (Ghosh et al., 2019;Xiong et al., 2016).Although various studies have investigated pore-scale modeling of SEAR and surfactant EOR either using pore-network modeling or by direct numerical modeling (Furtado & Skartlien, 2010;Landry et al., 2014;Li et al., 2017;Liu et al., 2018;Skartlien et al., 2011;Unsal et al., 2011), dynamic pore-scale modeling of two-phase flow, dissolution, and remobilization processes in SEAR has not received adequate attention.An existing work by Tsakiroglou et al. (2013) simplifies two-phase flow of SDS solution and n-C 10 (as contaminant) with a single phase flow where oil (NAPL) mobilisation is replaced with a dispersion-like phenomenological equation.Therefore this work is also not a two-phase modelling study.
In this work, we present dynamic pore-network modeling and upscaling of surfactant-enhanced aquifer remediation using a three-dimensional, unstructured pore network.The main focus of this work as presented in Section 2 is to investigate and upscale interphase mass transfer in SEAR processes by simulating two-phase flow, dissolution, remobilization of NAPL blobs, and surfactant effects on these processes.Section 3 aims to validate the pore network modelling code through comparison of computed results to previously published experimental data.In Section 4, the results of modeling of SEAR process for two types of surfactants, SDS and Triton X-100, are presented.Finally, we end the paper with summary and conclusions in Section 5.

Flow model
A single-pressure dynamic flow model is developed by incorporating both viscous force and capillary force into the pore network model.The pressures of water-invaded pores are calculated at each time step and the next throat filled by water is determined based on pore pressures and local capillary pressures.In this model, piston-like and snap-off fluid displacement mechanisms are considered according to the local aqueous phase velocity.While, trapping of NAPL does not occur during displacing phase in the piston-like mechanism, a large amount of NAPL could be trapped during the snap-off micro-flow mechanism (i.e., the aqueous phase flows through layers and occupies corners of the pore spaces).To investigate two-phase flow and determine initial NAPL saturation and distribution through the pore network, the network system is initially saturated with NAPL, and water is injected from the inlet under constant flow rate condition.
Pressure distribution through the pore network is calculated by implementing conservation of mass at water-invaded pores and solving a system of equations.Eq. 1 is a typical mass conservation equation for a single water-invaded pore  th .
∑  , = 0   =1 (1) where  , shows the local flow rate between two neighbouring pore bodies  ℎ and  ℎ , and   is the number of throats connected to a single pore  th .The local flow rate is calculated using Hagen-Poiseuille's law, Eq. 2 for invaded pore bodies  ℎ and  ℎ .
where   is flow conductivity of a throat, which connects pore  ℎ to pore  ℎ ,   is pressure difference between two sides of a throat, and    is the length of the throat between pore bodies  ℎ and  ℎ .The flow conductivity depends on the type of displacement mechanisms (e.g.snapoff and piston-like).In piston-like displacement,   is calculated using Eq. 3 (Al-Gharbi, 2004).
where   is the viscosity of water and  eff is the effective radius of throat which can be defined as (Bryant & Blunt, 1992): where   is the radii of the inscribed circle of each throat cross-section and   is the radii of a circle with the same cross-sectional area which is calculated using Eq. 5 and   the crosssectional area of throats: When the snap-off micro-flow is the fluid displacement mechanism, the flow conductivity   is calculated using Eq.7 (Wang & Dong, 2011).
where   is the cross-sectional area of the aqueous phase in the corners of throats which is approximated using Eq. 8 (Lago & Araujo, 2001).
where   is the interior angle of regular polygon with   corners, and   is the curvature radii of fluid-fluid interface in corners of throats which is calculated using Eq. 9 (Ramezanzadeh et al., 2019).
where  is a dimensionless shape factor which can be determined using area of throat (  ) and throat cross section perimeter (  ) as  =     2 ⁄ (Mason & Morrow, 1991).Finally,   in Eq.( 7) is a dimensionless parameter for flow of the aqueous phase in corners of throats, which can be calculated using (Zhou et al., 1997): Prior to calculation of pressure distribution, the type of displacement mechanism is necessary to determine at each time step.To predict snap-off mechanism, the force balance between viscous, gravitational, and capillary forces should be examined as follows (Brutsaert, 2005;Pennell et al., 1996): where   and   are the pressure force on the receding side and advancing side of a NAPL blob,   and  ℎ are radii of pore and throat, respectively, and  is the contact angle.It should be noted that the advancing and receding contact angles are assumed to be the same.Held and Celia (2001a) suggested that snap-off mechanism occurs when   =  ℎ   ⁄ is less than 0.5.The combination of Eq. 11 with the trapping condition proposed by Held and Celia (2001a)

results in
Eq. 12 which could be used as a new trapping condition to determine snap-off mechanism (if the condition is satisfied) and piston-like displacement mechanism (if the condition is not satisfied).
In Eq. 12,   is the entry capillary pressure in the throats.The new trapping condition enables us to consider the effects of pore structure, viscous force, and capillary force on trapping and remobilization of fluids.
The next invaded pore is decided by evaluating the maximum driving force, which is sum of viscous pressure difference in connecting throats and entry capillary pressure.Among the throats connected to an invaded pore body, the next invaded throat is the throat with maximum driving force.The entry capillary pressure in the throats can be calculated using Eq. 13 (Lago & Araujo, 2001;Mason & Morrow, 1991).
During the surfactant injection in the SEAR process, the interfacial tension is reduced by adsorbing surfactants to the interface, i.e., the van der Waals forces of attraction between surfactant and solvent molecules at the interface increases (Stebe & Lin, 2001).Therefore, the capillary pressure needs to be updated due to the effect of surfactants on the interfacial tension.
A thermodynamic analysis of the effect of surfactant at interface was investigated by Gibbs based on the increase in chemical potential of an adsorbed surfactant species at the interface, and the relationship between interfacial tension, , and surfactant concentration,  surfactant , is given by Szyszkowskie-Langmuir equation as (Karakashev et al., 2008): where  0 is NAPL-water interfacial tension,  is the gas constant,  is the absolute temperature, Γ  is the maximum surfactant concentration at the NAPL-aqueous phase interface,   is the Langmuir equilibrium adsorption constant, and CMC is the critical micelle concentration.
In our pore network model, two-phase displacement continues until NAPL ganglia are trapped in the pore spaces.NAPL entrapment and remobilization of entrapped NAPL blobs can be predicted using Eq.14.

Solute transport model
The concentration distribution for the invaded pores can be determined based on the pressure distribution through the network using applying the solute mass conservation.A typical pore  th with pressure of   can be connected either to water-invaded pores or to a NAPL saturated ones.
A typical pore  th , which is connected to pore  th can have a pressure higher than, lower than, or equal to   .Based on the pressure difference between each two neighboring pores and the type of displacement mechanisms, four different solute transport conditions can occur.A typical solute mass conservation equation for a single water-invaded pore  ℎ is as follows: where  , is flux of NAPL from pore  ℎ to pore  ℎ .The first type of NAPL flux occurs due to convection and diffusion when pressure of the adjoining pore  th is higher than that of the pore  ℎ .Hence,  , , can be calculated as (Zhao & Ioannidis, 2003): where  and   are the concentration and the aqueous solubility of NAPL, respectively.  is the normalized average solute concentration at the exit of each throat.In order to evaluate   , the diffusion of solute from NAPL-aqueous phase interface into the flowing wetting phase in corners of throats should be modeled in a simplified geometry as explained in Section 2.3.The second type of NAPL flux occurs by means of diffusion only when the adjoining pore body is saturated with NAPL species. , in this case can be obtained as: In the third type of NAPL flux, when pressure difference between the pore body  ℎ and the adjacent pore body  th is zero, solute flux occurs due to concentration difference between the two neighbouring pores.Hence,  , is calculated as: Finally, when the adjacent pore pressure   is greater than the pressure in pore  ℎ ,  , can be estimated using the following equation: Hence, based on the classification of NAPL flux, the concentration field can be determined for water-invaded pores.
where   0 is the equilibrium solubility of NAPL in water,   is the dissolution enhancement factor, and   is the micelle partition coefficient obtained based on experimental results.

Dissolution model
NAPL species is dissolved into the aqueous phase through two principal mechanisms: pore diffusion and corner diffusion (Zhou et al., 2000).As dissolution process is often dominated by corner diffusion mechanism in water-wet porous media, Dillard and Blunt, (2000) considered this mechanism for modelling of NAPL dissolution process.They assumed a simplified interface in which the corner diffusion is a one-dimensional diffusion from a flat NAPL surface into the aqueous phase surrounding it in corners of a water-wet throat.The pore-scale transport of NAPL species into aqueous phase across the simplified interface is evaluated using the following equation (Zhao & Ioannidis, 2003): here   = (, ) is the velocity profile of the aqueous phase flowing in the corner of the throat, and   is the molecular diffusion coefficient.To obtain NAPL flux in corners of the throat,   should be calculated using the following equation (Zhao & Ioannidis, 2003): where   is the normalized average solute concentration (the average solute concentration is normalized by aqueous solubility of NAPL) at the exit of the throat.To solve Eq. ( 22), uniform flow of the aqueous phase and 2D coordinate geometry of  are considered.The initial and boundary conditions are presented in Eqs. ( 23), ( 24), and (25).Eq. ( 26) can be derived from Eq. ( 22) based on the work done by Dillard and Blunt (2000).
(2).   indicates the dimensionless time   which is defined as   = (2  /3) to describe the duration of contact between NAPL and aqueous phase in corners of throats.
As the complementary error function in Eq. ( 28) is negligible,   can be estimated using the following equation at small values of   : The effect of surfactant on the level of solute concentration at the exit of each throat can be evaluated using thermodynamics of flow and reactive transport.
The adsorption of surfactants on NAPL-aqueous phase interface changes interface chemical potentials, and as a result, it lowers the interphase mass transfer from NAPL into aqueous phase.
For example, Mayer et al. (1999) showed experimentally that the mass transfer rate coefficients for surfactant-enhanced solubilization for TCE are at least an order of magnitude lower than those obtained for freshwater dissolution.Mayer et al. (1999) provided a potential explanation for this phenomenon by attributing the accumulation of surfactants near the NAPL-aqueous phase interface to retardation of the diffusion of either NAPL dissolved in the aqueous phase or micelles containing NAPL.However, to the best of our knowledge, there is no study to underpin this phenomenon.
Here the influence of adsorbed surfactant on interphase mass transfer was modeled using a thermodynamic procedure in accordance with Langmuir's energy barrier model (Bothe, 2015;Langmuir & Langmuir, 2002).Bothe (2015) has shown how interface concentrations and interface chemical potentials affect interphase mass transfer by changes in interface energy due to the adsorbed surfactants.Assuming incompressible liquid-liquid interface, reduction in dissolution rate due to the presence of surfactant can be predicted as (Bothe, 2015): where   * is the normalized average solute (NAPL) concentration in the presence of surfactants at the exit of each throat.Experimental data of mass transfer rate coefficients for different surfactant concentrations (Mayer et al., 1999;Ramezanzadeh et al., 2019) show that  changes with  surfactant until surfactant concentration reaches CMC.After CMC,  is a constant parameter with estimated value of −0.5 based on experimental data of Mayer et al. (1999).
Hence, the reduction in interphase mass transfer due to the adsorption of surfactants at NAPLaqueous phase interface can be evaluated as: The dissolution model can be applied to investigate and upscale interphase mass transfer in porous media.

Mass transfer rate coefficient
Upscaled mass transfer rate coefficient   can be evaluated using the definition of mass flux from entrapped NAPL into the aqueous phase as (Held & Celia, 2001b): here  ̅ is the average NAPL concentration in the pore network.If NAPL density   is constant, Eq. 31 can be rewritten as: where   and   are NAPL volume and total pore spaces in the pore network, respectively.  can be estimated as: in which   is total number of throats in the pore network.The mass transfer rate coefficient is obtained by insertion of Eq. ( 34) in Eq. ( 32) as follows: Mass transfer rate coefficients obtained by PNM can be upscaled according to Eq. ( 34) to predict interphase mass transfer with high accuracy at larger scales.

Properties of pore network and fluids
We develop a three-dimensional unstructured pore network in which each pore is represented as a cubic chamber with a number of square cross section throats connecting neighbouring pores based on their coordination number.The length (  ) of pore throat is obtained using locations of the two neighbouring pore bodies  ℎ and  ℎ .Also, the radii (  ) of pore throats are calculated based on the radii of the two neighbouring pore bodies  ℎ and  ℎ as follows (Joekar-Niasar et al., 2008): =     (  −0.9 +   −0.9 ) −0.9 where   and   are the radii of  th and  th pore bodies.
The size of the network is 1.15 mm (  ) × 1.15 mm (  ) × 7 mm (  ) with 5,552 pores.The average coordination number of the pore network is 5.38, and pore bodies in the network have average pore radii ( ̅  ) of 15.16 μm and standard deviation (  ) of 9.04 μm.
In our simulation, surfactant solutions were injected at different surfactant concentrations into an initially NAPL saturated pore network.Tetrachloroethylene (PCE) was selected as a dense nonaqueous phase liquid.Sodium dodecyl sulfate (SDS) and Triton X-100 (TX100) were used as surfactants.The physiochemical properties of fluids and the parameters required for evaluating adsorption of SDS and TX100 surfactants are listed in Table 1.

PNM Code Validation
The experimental data of Corapcioglu et al. (2009) and Pennell et al. (1996) are used to validate the PNM code.(1996) for a water-wetting porous medium (100-120 mesh Ottawa sand) during SDS flushing of PCE.This figure shows the predicted CDC is consistent with the experimental data of Pennell et al. (1996).
In addition to the PNM validation using CDCs, the interfacial area between NAPL and water after the end of two-phase flow could be validated using specific interfacial area (  : ratio of NAPL-water interfacial area to pore volume of the network).

Two-phase flow and NAPL entrapment
Surfactant can reduce the residual NAPL saturation as a result of IFT reduction (Alzahid et al., 2019), however, this depends on the capillary number.Here, the aqueous phase with surfactant (SDS) concentration of 0.0070 mol/L is injected at a constant flow rate of 0.00005 cc/min from the right-hand side to the left-hand side of the pore network (the boundary and initial conditions of this modelling were described in Section 2.3).
The pressure distribution of aqueous phase and distribution of NAPL through the pore network at different pore volume injection (PVI) are presented in Figs. 2 and 3, respectively.The breakthrough of the injected aqueous phase occurs at PVI of 0.24.The results show that the end of two-phase displacement is at PVI of about 1.15 (phase displacement is defined as the NAPL phase movement), and as shown in Fig. 2, there are still some pore bodies which are not invaded by the aqueous phase (marked as yellow points) before the end of two-phase displacement.The aqueous phase invades the pore spaces either using piston-like displacement or by snap-off mechanism until two-phase displacement ends.Pressure distribution does not change significantly after the end of NAPL phase displacement as a result of trapping of NAPL blobs in the pore bodies, and dissolution is the dominant mechanism in NAPL removal at this situation.
The simulation marks the residual NAPL saturation of 0.268 (for the case of SDS injection) after the end of two-phase displacement which is almost the same as the residual NAPL saturation of 0.273 for the surfactant-free case.Therefore, surfactant could not reduce residual NAPL saturation significantly at low flow rate.
Although the aqueous phase invades all pore bodies after the end of two-phase displacement due to dissolution (as shown in Fig. 2-f

Effect of flow rate on NAPL recovery and interfacial area
The effect of flow rate on NAPL recovery before and after the end of two-phase displacement at a constant surfactant (SDS) concentration of 0.0035 mol/L was investigated.Fig. 4a shows related recovery (the ratio of recovered NAPL after the end of two-phase flow to the initial saturation of NAPL) at different flow rates of 0.00005, 0.0001, 0.0005, 0.001, 0.005, 0.01, and 0.05 cm 3 /min.Residual NAPL saturations (NAPL saturations after the end of two-phase displacement) for different flow rates are also presented in Fig. 4b.The figure shows that increase in flow rate enhances NAPL recovery before the end of NAPL displacement due to increase in viscous force, i.e., the residual NAPL saturation after the end of two-phase displacement (at PVI of 0) decreases with increase of flow rate (capillary number).Fig. 4b indicates that surfactant can reduce residual NAPL saturation at the end of two-phase displacement at high flow rates.The residual NAPL saturations could be almost the same with and without surfactant at low capillary number (low flow rate).However, in terms of pore volume injection, NAPL recovery decreases with flow rate after the end of two-phase displacement as shown in Fig. 4a.This is because while NAPL mass transfer rate coefficient increases with flow rate (Aminnaji et al., 2019), the residence time decreases with flow rate which results in NAPL recovery reduction as a pore volume injection.
Figs. 4c and 4d show intrinsic interfacial area (  the ratio of NAPL-water interfacial area to NAPL volume), and specific interfacial area (  ) as a function of NAPL saturation at different flow rates from 0.00001 to 0.001 cm3/min.Intrinsic interfacial area is an index which could indicate the NAPL distribution through the porous media.Intrinsic interfacial area for a porous medium which consists of many NAPL droplets is higher than for a porous medium that consists of fewer NAPL droplets but with the same NAPL volume.As shown in Fig. 4c, the pore network model shows that   is independent of flow rate except for NAPL saturations of less than 1.5%.
The latter result indicates that the size and distribution of the trapped blobs are comparable at similar values of NAPL saturation when different flow rates are applied.This is consistent with the experimental results of Corapcioglu et al. (2009) which investigated the relationship between   and flow rate for different NAPL saturations.Hence,   can be applied in different porous media and flow rates in order to evaluate distribution of entrapped NAPL blobs and ganglia at similar NAPL saturations (Corapcioglu et al., 2009;Saripalli et al., 1997).Furthermore,   represents interfacial area in porous media which has a linear relation with NAPL saturation as shown in Fig. 4d.The linear relationship between   and NAPL saturation has been shown by other researchers (Corapcioglu et al., 2009;Joekar-Niasar et al., 2008).It should be mentioned that the small differences in the values of   for higher flow rates in Fig. 4c

Effect of surfactant type on NAPL recovery and interface surface area
The impact of SDS and TX100 on two-phase flow and NAPL entrapment is investigated using pore network modeling of two-phase flow at constant flow rate of 0.005 cm 3 /min with different surfactant injection concentrations of 0, 0.0035, 0.0070, 0.0087, and 0.0139 mol/L for SDS and concentrations of 0, 0.0001, 0.0002, 0.0005, and 0.0008 mol/L for TX100.The results of NAPL recovery during two-phase displacement are shown in concentrations after the end of two-phase displacement when the dissolution is the principal recovery process.The reason for this behaviour is that while SDS lowers interphase mass transfer (as it will be discussed in Section 4.4), increase in NAPL solubility due to SDS injection could not enhance NAPL dissolution, thereby reducing relative recovery.However, for SDS concentrations greater than CMC, interphase mass transfer remains constant, and NAPL recovery slightly increases due to enhancement of NAPL solubility.
In contrast to the SDS, simulations for TX100 injection above the CMC (0.0005 and 0.0008 mol/L) show a significant enhancement of NAPL remediation after the end of two-phase displacement.When TX100 was injected into the NAPL contaminated medium, NAPL recovery increased through dissolution due to a considerable enhancement in solubility as shown in Fig. 6c and Fig. 6d.In this case, due to high micelle partition coefficient and low CMC of TX100 (Harendra & Vipulanandan, 2011), the enhancement of solubility overcompensates the reduction in mass transfer rate coefficient.Therefore, as shown in Fig. 6c and Fig. 6d, for the TX100 concentrations above the CMC (0.0005 and 0.0008 mol/L), the recoveries are greater than the untreated case.However, for the TX100 concentrations below the CMC (0.0001 mol/L) and at CMC (0.0002 mol/L), the recovery curves fall below than that of the untreated simulation.This will be further explained in the Section 4.4.

Upscaled NAPL mass transfer rate coefficient
In this section, mass transfer rate coefficient (  ) is upscaled using the pore network modelling and Eq.35 in the presence of SDS and TX100.& Werth, 2003;Corapcioglu et al., 2009;Geller & Hunt, 1993;Mehdi Ramezanzadeh et al., 2019).This results confirms the hypothesis of Mayer et al. (1999) that the accumulation of surfactant at interface lowers interphase mass transfer resulting in reduction in NAPL effluent concentration.
However, at surfactant concentrations above the CMC,   and NAPL effluent concentration remain almost constant with the increase of SDS concentration.
In contrast to the SDS, when TX100 was injected through the pore network for NAPL remediation, mass transfer rate coefficient increased with TX100 concentration above the CMC (as shown in Figs.9c and 9d).Fig. 9c clearly indicates the variations (first decrease and then increase) of   with TX100 concentration at different NAPL saturations of 2, 5, and 10%.Fig.
9d also shows the similar variations in normalized effluent NAPL concentration with increase of TX100 concentration.
These results indicate SDS and TX100 have two different behaviours in mass transfer rate coefficient and NAPL recovery; and two reasons are suggested as being the root of these behaviours.Firstly, the micelle partition coefficient of SDS,   is too small to increase solubility of PCE in water.Harendra and Vipulanandan (2011) showed that solubility enhancement in the presence of PCE and SDS is limited by small value of   = 0.17 which is very lower than for TX100 (  =7.92).Secondly, CMC of SDS is 0.008 mol/L which is higher than that of TX100 (0.0002 mol/L).According to Eq. 20, the difference between surfactant concentration and CMC determines solubility enhancement, so the difference for SDS is far smaller than that for TX100 even at high surfactant concentration of 0.0139 mol/L (Harendra & Vipulanandan, 2011).

Remobilization
In most of the studies of interphase mass transfer processes, it is usually assumed that the shrinking NAPL blobs, which are trapped during two-phase displacement, is not remobilized during the dissolution process (Aminnaji et al., 2019;Corapcioglu et al., 2009;Khasi et al., 2019).However, the analysis of forces acting on a trapped blob in Section 2.1 indicates that the trapped NAPL blobs are not always stagnant.The trapped NAPL blobs remobilize when viscous forces are greater than capillary forces which hold NAPL blobs in pore bodies.NAPL blobs shrink due to dissolution resulting in reduction in radii of NAPL blobs, so both capillary and viscous forces are changed over time.When the viscous force exceeds the capillary force, the shrinking NAPL blobs can be remobilized.In this section, the onset of remobilization due to dissolution is studied by running simulations at different surfactant concentrations of 0, 0.0035, 0.007, 0.0087, and 0.0139 mol/L.The ratio of remobilized NAPL volume to trapped NAPL volume,   , versus capillary number are presented in Fig. 10.The onset capillary number of   = 3 × 10 −5 is estimated above which viscous forces significantly mobilized the shrinking NAPL blobs.The high ratios of remobilized NAPL, after the onset capillary number, suggest that remobilization due to dissolution should be considered in modeling with high capillary number.The results of the remobilization process are consistent with different experimental data (Datta et al., 2014;Duffield et al., 2003;Pennell et al., 1996).Various researchers have shown that remobilization occurs when   is close to 10 −5 , and the achievement of complete remediation occurs (regardless of gravity force) when   exceeds 10 −3 (Abriola et al., 2008;Boving & Brusseau, 2000;Li et al., 2007).

Pore network structure
Four different pore networks are used to study the effects of network structure (pore radii distribution) on upscaled mass transfer rate coefficient as well as NAPL recovery of SEAR.Fig. 11 illustrates pore radius distribution of the adopted networks.Porosities of the generated networks with the same bulk volumes are 1.99, 0.9, 2.72, and 3.72 %, respectively.The results of NAPL recovery, effluent NAPL concentration, mass transfer rate coefficient, and effective specific interfacial area (which is defined as    ⁄ ) are compared for the networks at injection velocity of 0.001 cm 3 /min and surfactant concentration of 0.0035 kg/m 3 in Fig. 12. Fig. 12(d) shows that the heterogeneity of the network can have an important influence on NAPL distribution as effective specific interfacial area changes significantly in the networks.As network number 3 and 4 in Fig. 12 have more similar pore radius distribution, the effective specific interfacial area is almost the same in these cases.

Summary, Conclusions and Future Works
In this work, we have studied the effect of surfactant solutions on immiscible two-phase flow, enhanced dissolution process, and remobilization of entrapped NAPL blobs in SEAR using dynamic pore network model in a three-dimensional unstructured pore network.To the best of our knowledge, the pore scale modelling of SEAR has received inadequate attention in literature.
To understand how surfactants influence SEAR process, we simulated the injection of surfactant solutions with different surfactant concentrations into an initially NAPL saturated pore network.
Corner flow and micro-flow mechanisms including snap-off and piston-like were considered in the simulation of two-phase flow.Furthermore, NAPL entrapment and remobilization were evaluated using force analysis in order to predict the onset of remobilization for different surfactant concentrations.Besides, corner diffusion mechanism was applied in the modeling of interphase mass transfer to evaluate NAPL dissolution as the principal mass transfer process in saturated media.We also applied PNM to investigate and upscale interphase mass transfer.The evaluated mass transfer rate coefficients and thermodynamic analysis showed that although surfactants enhance NAPL recovery during two-phase flow, surfactant-enhanced remediation of NAPL contaminated media is highly dependent on surfactant type.Based on the results obtained in this work, the following conclusions can be drawn:  NAPL recovery and mass transfer rate coefficients may either decrease or increase with surfactant concentrations depending on surfactant type after the end of two-phase displacement where dissolution is principal mechanism.Surfactants with high CMC and low micelle partition coefficients significantly reduce NAPL recovery after the end of two-phase flow because mass transfer rate coefficients decrease due to considerable changes in interface chemical potentials.However, surfactants with low CMC and high micelle partition coefficients improve NAPL recovery, because an enhancement of solubility at surfactant concentrations greater than CMC offsets the mass transfer reductions.
 We determined the onset of remobilization and the capillary number, where complete NAPL removal due to remobilization occurs, using force analysis on trapped NAPL blobs.The results of remobilization are consistent with the experimental results in the literature.
 A capillary desaturation curve (CDC) was developed using pore-scale modeling of twophase flow and NAPL entrapment at different flow rates.The results of initial residual saturations are consistent with the experimental results of Pennell et al. (1996).
 We found that intrinsic interfacial area, as an indicator of NAPL distribution in porous media, is independent of flow rate, except for NAPL saturations of less than 1.5%.This parameter increases significantly in the presence of surfactant at lower NAPL saturations due to changes in NAPL distribution in porous media.The results are in agreement with the experimental work done by Corapcioglu et al. (2009)  The results of pore network modeling of SEAR for different networks show that effective specific interfacial area, mass transfer rate coefficient, and NAPL recovery increase as average pore size decreases in the networks.
We used SDS surfactant in this work to highlight the role of surfactants on interphase mass transfer.To investigate the effect of surfactant type on interphase mass transfer and solubilization after the end of two phase flow, the pore network modeling of SEAR in the presence of TX100 was also simulated.The results indicate that the role of surfactants in aquifer remediation is complex, and one should consider all of surfactant effects to find the optimal surfactant concentration and composition in SEAR.Furthermore, SDS alone may not enhance aquifer remediation considerably, while the presence of TX100 can accelerate removal of NAPL through dissolution.Experimental studies have also indicated that the removal efficiencies by TX100 and SDS-TX100 are much larger than those by SDS (Mazen and Radzuan 2009).Mazen and Radzuan (2009) stated that SDS (anionic surfactant) adsorption was not detected to any of the adsorbent samples in their experiments, and surfactant adsorbed/Kg adsorbent was lower for TX100-SDS mixtures in comparison to TX100.Therefore, an optimal mixture of SDS-TX100 due to different adsorption behaviour (which is not considered in our simulation results) might result in maximum removal efficiency in SEAR and surfactant enhanced oil recovery (EOR) applications.

Acknowledgements
Fig. 1a shows the capillary desaturation curves (CDCs) obtained by pore network model and experimental data.The CDC was predicted by pore network modeling of two-phase flow and NAPL entrapment at different flow rates during surfactant flushing at a constant SDS concentration of 0.0035 mol/L.The experimental data used in Fig. 1a is from Pennell et al.
Fig. 1.(a) Capillary desaturation curves showing the ratio of residual saturation to initial residual

Fig. 2 .Fig. 3 .
Fig. 2. Pressure distribution of the aqueous phase in Pascal through the pore network at (a) PVI of 0.24, Fig. 4. (a) Relative NAPL recovery (the ratio of recovered NAPL after the end of two-phase flow to initial

Fig. 5 .
Fig. 5. NAPL recovery as a function of PVI with the injection flow rate of 0.005 cm 3 /min during the two-

Fig. 7a and
Fig. 7a and 7b show   and   as functions of NAPL saturation for different surfactant (SDS) concentrations, respectively.As   is an indicator of NAPL distribution in porous media, the significant changes in   at low NAPL saturations show that surfactant flooding could change Fig. 8a shows the NAPL mass transfer rate coefficients versus NAPL saturation at different injection flow rates at a constant SDS concentration of 0.0035 mol/L.The results indicate the increase of mass transfer rate coefficients with flow rate.The reason for this behaviour is that the increases of flow rate causes reduction in residence NAPL concentration which results in increase of the driving force (  − ) and NAPL dissolution rate.In addition to the increase of NAPL dissolution rate, the convective mass transfer in the pore throats increases with flow rate resulting in increase of   .The increase of mass transfer rate coefficient with flow rate has been indicated by other researchers (Chomsurin

Fig. 8b represents
Fig. 8b represents the effluent NAPL concentration normalized by solubility of NAPL (at SDS

Figs. 9a and
Figs. 9a and 9b show the effect of surfactant (SDS) concentration on   and normalized NAPL effluent concentration at constant flow rate of 0.005 cm 3 /min.The results show the reduction in Fig. 9. (a) Mass transfer rate coefficients versus SDS concentrations at different NAPL saturations of 2, 5,

Fig. 10 .
Fig. 10.The ratio of remobilized NAPL volume to trapped NAPL volume versus   for different surfactant concentrations.
Fig. 12(c)  indicates that mass transfer rate coefficient increases as average pore size decreases in the networks.The latter result is consistent with experimental data inCho et al. (2005) in which interphase mass transfer increases as grain size decreases.

Table 1 :
Physicochemical properties of fluids Morteza Aminnaji and Masoud Babaei contribution into this research was provided by the Engineering and Physical Science Research Council, Grant No. EP/R009678/1.MA and MB thank this organization.The data presented in this paper (simulation results of pore networks) are freely available to download at https://doi.org/10.6084/m9.figshare.12776249.