Vol.:(0123456789) Transport in Porous Media (2023) 146:639–667 https://doi.org/10.1007/s11242-022-01879-0 1 3 Effect of Pore Space Stagnant Zones on Interphase Mass Transfer in Porous Media, for Two‑Phase Flow Conditions H. Gao1  · A. B. Tatomir1,4 · N. K. Karadimitriou2 · H. Steeb2,3 · M. Sauter1 Received: 8 June 2022 / Accepted: 30 October 2022 / Published online: 10 November 2022 © The Author(s) 2022 Abstract Interphase mass transfer is an important solute transport process in two-phase flow in porous media. During two-phase flow, hydrodynamically stagnant and flowing zones are formed, with the stagnant ones being adjacent to the interfaces through which the inter- phase mass transfer happens. Due to the existence of these stagnant zones in the vicinity of the interface, the mass transfer coefficient decreases to a certain extent. There seems to be a phenomenological correlation between the mass transfer coefficient and the extent of the stagnant zone which, however, is not yet fully understood. In this study, the phase- field method-based continuous species transfer model is applied to simulate the interphase mass transfer of a dissolved species from the immobile, residual, non-aqueous phase liquid (NAPL) to the flowing aqueous phase. Both scenarios, this of a simple cavity and this of a porous medium, are investigated. The effects of flow rates on the mass transfer coefficient are significantly reduced when the stagnant zone and the diffusion length are larger. It is found that the stagnant zone saturation can be a proxy of the overall diffusion length of the terminal menisci in the porous medium system. The early-stage mass transfer coefficient continuously decreases due to the depletion of the solute in the small NAPL clusters that are in direct contact with the flowing water. The long-term mass transfer mainly happens on the interfaces associated with large NAPL clusters with larger diffusion lengths, and the mass transfer coefficient is mainly determined by the stagnant zone saturation. Article highlights • Flow rates can hardly affect the mass transfer coefficient when the stagnant zone and the diffusion length are large. • The stagnant zone saturation can be a proxy of the overall diffusion length of the menisci in the porous media system. * H. Gao huhao.gao@uni-goettingen.de 1 Department of Applied Geology, University of Göttingen, 37073 Göttingen, Germany 2 Institute of Applied Mechanics (CE), University of Stuttgart, 70569 Stuttgart, Germany 3 Stuttgart Center for Simulation Science (SC SimTech), Pfaffenwaldring 5a, 70569 Stuttgart, Germany 4 Independent researcher, Göttingen, Germany http://orcid.org/0000-0003-0123-9155 http://crossmark.crossref.org/dialog/?doi=10.1007/s11242-022-01879-0&domain=pdf 640 H. Gao et al. 1 3 • The correlations between the mass transfer coefficient and the stagnant zone saturations are built, based on simulations at different flow rates. Keywords Interphase solute mass transfer · Phase field method · Two-phase flow · Porous media List of symbols Awn Interfacial area [m2] Adif wn Diffusive interface [m3] awn Specific interfacial area [m−1] bwn Diffuse interface thickness [m] cwater Average aqueous phase concentration [kg·s−1  m−3] cNAPL Average non-aqueous phase concentration [kg·s−1  m−3] ck Solute concentration [kg·s−1  m−3] cs Equilibrium solute concentration [kg·s−1  m−3] Dk Molecular diffusion coefficient [m2/s] fmix Free energy density [J/m3] G Chemical potential [J/m3] h Slide length of mesh triangle [m] J Solute diffusion flux [mol·s−1  m−2]] kf Mass transfer coefficient L Diffusion length [m] Lc Characteristic length M Interphase transferred solute mass [kg] p Pressure [pa] Pow Partitioning coefficient Pe Péclet number Re Reynolds number Sh Sherwood number SNAPL NAPL saturation Ssw Saturation of stagnant water t Time [s] u Velocity [m/s] uin Inflow velocity [m/s] Vpm Volume of the porous media [m3] Vf Volume fraction of fluid Vwater Region in domain occupied by water [m3] VNAPL Region in domain occupied by NAPL [m3] Vinterface Region in domain occupied by diffusive interface [m3] Δt Time interval [s] � Phase variable λ Magnitude of the mixing energy [N] ε Capillary thickness [m] σ Surface tension [N/m] γ Mobility [m3s/kg] χ Characteristic mobility [m·s/kg] ρ Density [kg/m3] 641Effect of Pore Space Stagnant Zones on Interphase Mass Transfer… 1 3 μ Viscosity [pa·s] �w Contact angle 1 Introduction Interphase mass transfer in porous media two-phase flow systems is an important process relevant in a wide range of scientific and engineering applications including the reme- diation of aquifers contaminated by non-aqueous phase liquids (NAPL) (Brusseau, 1992; Miller et  al., 1990), geological storage of carbon dioxide (Schaffer et  al., 2013; Tatomir et al., 2016), enhanced oil recovery (Graveleau et al., 2017), etc. In some of these applica- tions, the dissolved chemical species (solute) in one phase of fluid may cross the interface and enter the other phase. The interphase mass transfer can be classified into different con- ditions depending on the fluid phase as the source of the species is mobile or immobile, wetting or non-wetting, according to Agaoglu et al. (2015). This study focuses on the inter- phase mass transfer process of the dissolved species from the residual immobile wetting liquid into the mobile non-wetting liquid in porous media. 1.1 Interphase Mass Transfer The interphase mass transfer in porous media has been investigated both at the pore scale and at the macro-scale. At the pore scale, the pore structures and fluid–fluid interfaces are explicitly described in the studies. The interphase mass transfer occurs when the compound of interest has two individual chemical potentials with respect to each phase on both sides of the interface, and the thermodynamic equilibrium is violated (Weber et al., 1991). The upscaling techniques (e.g., volume averaging) have been widely used to derive the macro- scale equations from the pore-scale physics (Coutelieris et al., 2006; Kechagia et al., 2002; Quintard and Whitaker, 1994; Soulaine et al., 2011). Physical parameters are averaged over the representative element volume (REV) in the macro-scale studies. At the macro-scale, a common approach to determine the interphase mass transfer rate is to define it as the product of a mass transfer coefficient that is unique for the system, the interfacial area, and the driving force, which is determined by the difference between the current solute con- centration in the solvent phase and the solute concentration in the solvent phase at thermo- dynamic equilibrium (Agaoglu et al., 2015; Miller et al., 1990). Thus, the interphase mass transfer rate can be formulated as: where M (kg) is interphase transferred solute mass, t (s) is time, Awn (m2) is the interfacial area, c (kg/m3) is the current solute concentration, cs (kg/m3) is the equilibrium solute con- centration, and kf (m/s) is the mass transfer coefficient. The mass transfer coefficient can be expressed in a dimensionless form as the Sherwood number: where Lc (m) is characteristic length equal to the mean grain diameter, and D (m2/s) is the solute molecular diffusion coefficient. In some studies, the interfacial area is not available, (1) dM dt = kf Awn(cs − c) (2)Sh = kf Lc D 642 H. Gao et al. 1 3 and thus a lumped mass transfer coefficient is introduced, which incorporates the above mass transfer coefficient kf and the interfacial area (Brusseau, 1992; Imhoff et  al., 1994; Miller et  al., 1990). Since the interfacial area can be directly measured/estimated with pore-scale numerical simulation approaches, this study mostly focuses on the explicitly formulated mass transfer coefficient rather than the lumped mass transfer coefficient. The mass transfer coefficient for a given system will approach a constant value after a certain period of time, due to the establishment of equilibrium in the driving forces and mass transfer rates (Miller et al., 1990; Soulaine et al., 2011). This mass transfer coefficient can be affected by many factors in the two-phase flow process, including interfacial area, saturation, flow rates, phase entrapments, NAPL distribution, NAPL mobilization, fluid density and viscosity, the physical properties of the porous medium, e.g., wettability, mean grain diameter, pore size distribution (Agaoglu et al., 2015), and the chemical properties of the NAPL and dissolved solute (Essaid et al., 2015). Flow velocity is considered an impor- tant parameter for interphase mass transfer since it determines the renewal of the com- pound near the interface (Graveleau et al., 2017). Numerous researches have been reported to build the correlations between the flow velocity (expressed as the Reynolds number or the Péclet number) and the mass transfer coefficient (expressed as the Sherwood number) (Aydin Sarikurt et al., 2017; Imhoff et al., 1994; Kennedy and Lennox, 1997; Miller et al., 1990; Powers et al., 1994b, 1994a, 1992). The proposed correlations are different depend- ing on the properties of the porous material, the fluid saturations, the type of the NAPL, etc. 1.2 Effects of Stagnant Zones The stagnant zones are commonly defined as the regions with a very low flow veloc- ity, casting diffusive mass transport (Karadimitriou et  al., 2016). The non-wetting phase belonging to the stagnant zones and the flowing zones, and the residual wetting phase are distinguished in Fig. 1. It was found that the presence of stagnant zones can influence the tracer transport in unsaturated systems in the column experiments, i.e., larger longitudinal dispersion, early breakthrough and long tailings (Bond and Wierenga, 1990; Smedt and Wierenga, 1984; van Genuchten and Wierenga, 1976). The recent pore-scale micro-model experiments and numerical simulations further elaborate the critical roles the stagnant zones play in solute transport for two-phase fluids systems, e.g., affecting the dispersion Fig. 1 Sketch of the distribu- tion of the two phases of fluids involved in the solute interphase mass transfer in the porous media (black arrows show the flow direction), where flowing zones and stagnant zones of the non- wetting fluid are distinguished with the frame of solid lines and dashed lines 643Effect of Pore Space Stagnant Zones on Interphase Mass Transfer… 1 3 coefficient, in correlation to the saturation topology and the relative permeability, etc. (Aziz et al., 2018; Dou et al., 2022; Hasan et al., 2020, 2019; Karadimitriou et al., 2017, 2016). Besides, the existence of the stagnant zones is also an important factor affecting the interphase mass transfer rate and the mass transfer coefficient. With the use of pore-scale numerical simulation with the lattice Boltzmann method, Knutson et al. (2001) found that the interphase mass transfer is limited by the diffusive transport process through a thin zone of stagnant water at NAPL/water interface. Zhang et al. (2002) used magnetic resonance imaging to study NAPL dissolution in water-saturated flowthrough columns packed with angular silica gel of two different sizes. They found that the mass transfer coefficient was smaller in experiments with finer grains. This was attributed to the fact that finer grains led to the formation of a larger number of NAPL ganglia and an increment of the quantity of stagnant water adjacent to the interfaces. Thus, the extent of flowing water bypass (without contact to the interface) was increased in their system. These results are in agreement with the micromodel, fluorescence microscopy observations of Chomsurin and Werth (2003). They also found that the Sherwood number is highly dependent on the average diffusion length scale associated with the stagnant zones. Aminnaji et al., (2019) presented a pore- network model to study the NAPL dissolution, where links between the number of stagnant throats for different networks and the corresponding mass transfer coefficients are studied. They concluded that the NAPL mass transfer coefficient mainly depends on the distribution of the NAPL, rather than the interfacial area or the number of stagnant throats. Despite the important findings in the above research studies, they were mostly focused on NAPL dissolution problems, where the species of interest is the NAPL itself rather than another species dissolved in it. One difference between these two conditions is the equilib- rium solute concentration (cs) in Eq. (1). For NAPL dissolution, the equilibrium concentra- tion (cs) is a constant equal to the solubility limit, while for interphase transfer of a dilute species (other than NAPL itself) the equilibrium concentration (cs) is a variable dependent on its concentration in the NAPL. Thus, the mass transfer coefficients obtained from these two conditions can be different. To the authors’ knowledge, the effects of stagnant zones on the interphase transfer of dilute species in porous media were not yet studied. Furthermore, the stagnant zone is not directly quantified in these studies, e.g., by a stagnant zone satura- tion (the ratio of the volume of the stagnant zone to the total volume of pore space). Thus, the parameters representing the extent of the stagnant zones (e.g., stagnant zone saturation) have not yet been used to build the correlations for the characterization of interphase mass transfer. It is worth noting that the effects of stagnant zones cannot be represented with the fluid saturation in the correlations, because it has been shown that the volume of the stag- nant zone depends in a non-monotonic way on the saturation for a fixed saturation topol- ogy (Dou et al., 2022; Hasan et al., 2020; Karadimitriou et al., 2017, 2016). How the mass transfer coefficients change with respect to different volumes of the stagnant zone at differ- ent saturations is not fully understood. Therefore, the correlation between the stagnant zone saturation and the mass transfer coefficient will be studied here by performing a pore-scale investigation using direct numerical simulation. 1.3 Numerical Methods at Pore Scale The pore-scale numerical simulations of the interphase mass transfer have been done with the pore network models (PNM) (Dillard et al., 2001; Dillard and Blunt, 2000; Held and Celia, 2001; Jia et al., 1999; Zhao and Ioannidis, 2003). However, the previous studies using PNM are mostly regarding to the NAPL dissolution, and the simplified representations of the pore 644 H. Gao et al. 1 3 space may make the calculation of the stagnant zones saturation inaccurate. As an alternative, computational fluid dynamics (CFD) approaches can be used for simulation of the interphase mass transfer in complex (natural) pore structures. The CFD approaches, e.g., the volume-of- fluid method (VOF), the level-set method (LSM), the phase-field method (PFM), solve the Navier–Stokes equations in a discretized domain, and the interface between the two fluids is represented by an indicator function (Graveleau et al., 2017; Maes and Soulaine, 2018). One major challenge for the simulation of species interphase mass transfer with the CFD method is the handling of the concentration jump at the fluid–fluid interface (Graveleau et al., 2017). This problem is solved by Haroun et al., (2010a) with a new single-field Continuous Species Transfer (CST) formulation. The CST formulation, developed in the framework of the VOF method, introduced a constant partitioning coefficient in the transport equation to describe the thermodynamic equilibrium of the species at the fluid–fluid interface (Haroun et al., 2010a). The VOF-CST model was successfully applied to simulate the interphase mass transfer in liq- uid film flow along corrugated surfaces (Haroun et al., 2012, 2010b), in bubble flow (Deis- ing et  al., 2016; Marschall et  al., 2012), and in subsurface two-phase flow problems with moving contact lines (Graveleau et al., 2017). Yang et al. (2017) found that the solute mass partitioning between two phases is not accurately simulated by the VOF-CST model at high Péclet numbers, due to the numerical diffusion in the phase concentration. By adding a com- pressive term in the transport equation, Maes and Soulaine (2018) managed to improve the VOF-CST model by significantly diminishing the numerical errors. Gao et al. (2021a) further derived the CST formulations based on the PFM, and the model was labelled as PFM-CST. Both approaches have their pro and con. The PFM with a smooth indicator function can better simulate the curvature and smooth physical quantities at the interface (Akhlaghi Amiri and Hamouda, 2013). However, the PFM is generally more computationally expensive due to the necessity to resolve the diffusive interface on a fine grid (Maes and Soulaine, 2018). 1.4 Objectives In this study, we employ the previously developed PFM-CST models to simulate the inter- phase mass transfer of a dissolved dilute species in the residual immobile NAPL, into the flowing water, in a cavity and in a porous medium. We focus on the understanding of how the stagnant zones affect the mass transfer coefficient at the pore scale. The novelty of the paper is that the correlations between the mass transfer coefficient (for interphase transfer of dilute species) and the parameters representing the extent of the stagnant zones (the stagnant zone saturation and the average diffusion length) are first time investigated, at different water flow rates and different saturation topologies. For the porous medium domain, different phase dis- tribution patterns (with different volumes of stagnant zones) are obtained from the displace- ment processes at different flow rates in the preparatory step. The paper is organized as follows: Sect. 2 reviews the mathematical model. In Sect. 3, we explain the numerical implementation and the details of the numerical setup and data process- ing. The simulation results and discussions are demonstrated in Sect. 4, and in Sect. 5, finally, we list the main conclusions. 645Effect of Pore Space Stagnant Zones on Interphase Mass Transfer… 1 3 2 Mathematical Model 2.1 Phase‑Field Method The basic theory of the PFM is first reviewed briefly here according to the studies, e.g., (Alpak et al., 2016; Basirat et al., 2017; Jacqmin, 1999; Yue et al., 2004). In the PFM, sharp interface between two fluid phases is replaced by a diffusive interfacial layer, with a small but finite thickness (Basirat et al., 2017). Following Yue et al. (2004), a phase variable ( � ), changing smoothly between -1 and 1, is introduced to describe the composition of the fluid mixture at the interface. Thus, the volume fractions of two fluids (Vf) are defined as 1−� 2 and 1+� 2 , respec- tively. The � profile across the interface is determined by the free energy density (Alpak et al., 2016). The free energy density (mixing energy) (J/m3) for isothermal mixing of two fluids can be expressed in the Ginzburg–Landau form (Yue et al., 2006): where λ (N) is the magnitude of the mixing energy and ε (m) is a capillary thickness which determines the thickness of the diffusive interface. The free energy density (in Eq. 3) is composed of the surface energy (the first term) and the bulk energy (the second term) (Basirat et al., 2017). The chemical potential G (J/m3) is calculated by the variation of the free energy with respect to �: where Ω is the spatial domain occupied by the two fluid phases. The interfacial tension, defined as the total free energy at the interface per unit area of the interface, is in relation to the capillary thickness and the magnitude of the mixing energy through (Jacqmin, 1999; Yue et al., 2004): Assuming that the diffusive fluid flux is proportional to the gradient of the chemical poten- tial, fluid mass conservation can be described by the Cahn–Hilliard equation (Cahn and Hill- iard, 1959): where u is the fluid velocity vector and γ (m3s/kg) is the mobility. The mobility is neces- sary to be large enough to keep a constant interfacial thickness, but small enough not to suppress the flow (Jacqmin, 1999). The mobility can be formulated as a function of the capillary thickness and a tuning factor: � = ��2 (Akhlaghi Amiri and Hamouda, 2013). 2.2 Two‑Phase Flow Dynamics Momentum conservation is governed by the incompressible form of the Navier–Stokes equation, where the surface tension force is considered as a body force according to Yue (3)fmix(�,∇�) = 1 2 �|∇�|2 + � 4�2 ( �2 − 1 )2 (4)G = δ∫ Ω fmixdx δ� = −�∇2� + � �2 �(�2 − 1) (5)� = 2 √ 2� 3� (6) �� �t + � ⋅ ∇� = ∇ ⋅ (�∇G) 646 H. Gao et al. 1 3 et al. (2006) ( − δ ∫ Ω fmixdx δx = G∇� ) , and the gravity effects are neglected due to the horizon- tal setup in Sect. 3.2: where u (m/s) is the fluid velocity vector, p (Pa) is pressure, I is the second-order identity tensor, μ (Pa·s) is the dynamic viscosity, and ρ (kg/m3) is the density. It is worth noting that the incompressible formalism of the Navier–Stokes equation can be applied for the fluid mixture at the interface as long as the interface thickness is adequately small and the results are unaffected (Gao et al. 2021b). Besides, it is assumed that the fluid properties are not affected by the solute transport, meaning ρ and μ are independent of the solute concen- tration. This is reasonable because the solute has small mass compared to the bulk fluid phases. μ and ρ are calculated as: where Vf,w and Vf,nw indicate the volume fraction of the wetting fluid and the non-wetting fluid, respectively. A no-slip boundary is applied at the grain surfaces, which implies that u = 0 in the momentum equation for wall boundaries. The contact angle on the solid wall is defined as: where n is the (outward) normal vector to the wall and �w is the contact angle. With the dif- fusive fluxes between the bulk fluids (Eq. 6), the contact line on the wall results effectively in slip (Gao et al., 2021a). Thus, additional settings about the contact line moving, which are necessary for VOF-CST models reported by Graveleau et al. (2017), are not required here. 2.3 Solute Transport The continuous species transfer model, based on the phase-field method, is derived and explained in Gao et al. (2021a, 2021b). The species interphase mass transfer includes sev- eral major mechanisms, i.e., advection, molecular diffusion, and partitioning. The chemical species considered here is considered as solute in both of the fluid phases, which means the solubility limit is not exceeded and no precipitation occurs. In this section, we review in brief the equations implemented in the solute transport model. Firstly, for species k dis- solving in fluid phase α with the mass concentration of ca,k , the solute mass conservation equation writes (Gao et al., 2021a): where the dominant transport process due to the bulk fluid motion (i.e., advection) is described in Eq. (11). The local thermodynamic equilibrium of the solute at the interface is (7)� ( �� �t + � ⋅ ∇� ) = ∇ ⋅ [ −p� + � ( ∇� + (∇�)T )] + G∇� (8)∇ ⋅ � = 0 (9)� = �wVf ,w + �nwVf ,nw and � = �wVf ,w + �nwVf ,nw (10)� ⋅ ∇� = cos�w|∇�| (11) { �(Vf ,wcw,k) �t + ∇ ⋅ (uVf ,wcw,k + 0.5�∇Gcw,k) = 0 �(Vf ,nwcnw,k) �t + ∇ ⋅ (uVf ,nwcnw,k − 0.5�∇Gcnw,k) = 0 647Effect of Pore Space Stagnant Zones on Interphase Mass Transfer… 1 3 described by a partitioning coefficient (Pow,k), which is defined as the concentration ratio of the species in each phase, at equilibrium. where cnw,k and cw,k are the concentration of the solute in the non-wetting phase and the wetting phase, respectively. The single-field formulation of the species concentration for both phases in the entire domain is obtained by defining the global concentration (ck) as: Therefore, the sum of the mass conservation Eq. (11) for species k in both phases writes: The mechanisms of molecular diffusion are considered by a diffusion term, where the molecular diffusion flux for the mixture is defined as in (Haroun et al., 2010a): As Vf,w + Vf,nw = 1, here only Vf,w is used in the expression, and the molecular diffusion coef- ficient is expressed as: It is reported that the harmonic interpolation of the diffusion coefficient for the mixture is more robust than a simple mixing rule (Deising et al., 2016; Maes and Soulaine, 2020). The diffusion flux (Eq. 16) consists of two parts: one flux driven by the concentration gradient, and the other flux resulting from the solubility law, which is in the direction normal to the inter- face. This flux resulting from the solubility law governs the partitioning of the solute between the two phases (Gao et al., 2021a). With the definition of Eq. (13) and Eq. (14), one obtains: The final governing equation for the global concentration of species k is obtained by incor- porating Eq. (18) and Eq. (16) into Eq. (15): (12)Pow,k = cnw,k cw,k (13)ck = Vf ,wcw,k + Vf ,nwcnw,k (14) �c k �t + ∇ ⋅ [ �c k + 0.5�∇G ( c w,k − c nw,k )] = 0 (15)�k = Dk(Vf ,w∇cw,k + Vf ,nw∇cnw,k) = Dk[∇ck − ( cw,k−cnw,k ) ∇Vf ,w] (16)Dk = 1 Vf ,w Dw,k + Vf ,nw Dnw,k (17)cw,k−cnw,k = ck(1 − Pow,k) Vf ,w + Pow,k(1 − Vf ,w) (18) �ck �t + ∇ ⋅ [ �ck + 0.5�∇G ck ( 1 − Pow,k ) Vf ,w + Pow.k ( 1 − Vf ,w ) ] −∇ ⋅ [ Dk ( ∇ck − ∇Vf ,w ck ( 1 − Pow,k ) Vf ,w + Pow,k ( 1 − Vf ,w ) )] = 0 648 H. Gao et al. 1 3 3 Numerical Details We explain the model implementation in Sect. 3.1, the setup of the model in Sect. 3.2, and the data treatment process in Sect. 3.3. 3.1 Numerical Implementation The model is implemented in the software COMSOL Multiphysics™ (V5.5), using the basic module and the CFD module. COMSOL employs the finite element method for spa- tial discretization (Gao et al., 2021b). The linear solver PARDISO is used to solve the dis- cretized partial differential equations, and the backward Euler method is used for calcu- lating the time stepping. The initial time step and maximum time steps are automatically adjusted so that singularities are avoided (Akhlaghi Amiri and Hamouda, 2013). The mesh generator in COMSOL uses the advancing front algorithm to calculate the tessellation, and the mesh elements are in triangles (Gao et  al., 2021b). The computer has a single CPU with 8 cores, operating at 4.3 GHz, and 128 GB RAM. The numerical implementation of the model in COMSOL has been verified in the previous study by comparing it with one- dimensional analytical solutions for several processes, including solute advective diffusion and transient interphase mass transfer processes (Gao et al., 2021a). 3.2 Model Setup The effect of the stagnant zones on the interphase mass transfer is studied based on two scenarios: (1) a single cavity and (2) a porous medium. The domains for simulation in both scenarios are horizontal setups, and thus, the gravity is irrelevant. The proper- ties of the fluids are considered to be the same in both scenarios: ρwater = 1000  kg/m3, μwater = 1 × 10–3 Pa·s, ρNAPL = 840 kg/m3, and μNAPL = 6.5 × 10–3 Pa·s. The water properties correspond to pure water at 20 °C, and the NAPL properties correspond to a North Sea crude oil (Mahani et  al., 2014). The interfacial tension σ = 0.01  N/m, and the NAPL is assumed as the wetting phase with the contact angle of θw = 135° according to the simula- tion by Maes and Soulaine, (2018). The solute for simulation is a contaminant k, with the partitioning coefficient Pow = 0.5, the initial concentration of 1 kg/m3, and the molecular diffusion coefficient Dk = 1 × 10–9 m2/s (Maes and Soulaine, 2018). The contaminant k is initially dissolved in the NAPL, and it is being removed from the domain by mass transfer into the flowing water. The characteristic mobility, χ, is chosen between 10 and 100, so that the thickness and configuration of the interface is ensured to be invariant during the mass transfer process, according to Appendix 1. The capillary thickness (ε) and the mesh size (h) are chosen small enough to ascertain mesh convergence and model convergence (Gao et al., 2021a). On the other hand, the mesh size (h) cannot be too small to make the compu- tational cost unaffordable. 3.2.1 Study of Interphase Mass Transfer in a Cavity The interphase mass transfer process is first studied in a single cavity, with only one inter- face being involved. The intention of the numerical setup is to fix the interfacial area, initial solute concentration in NAPL and NAPL volume in all simulations, so that the changes in 649Effect of Pore Space Stagnant Zones on Interphase Mass Transfer… 1 3 the mass transfer coefficient are mainly studied in relation to different sizes of the stagnant zones (space occupied by the immobile water) and different flow velocities. The setup of the cavity is illustrated in Fig. 2, where the top rectangular represents a flow channel and the bottom rectangular represents a throat with stagnant fluids. Water is injected at the left- hand boundary at inflow velocities of 1 × 10–2 m/s, 1 × 10–3 m/s, and 1 × 10–4 m/s, with the right-hand boundary as the outlet with a constant pressure p = 0. The initial solute concen- tration in water and the concentration at inlet are both zero (left boundary with constant concentration ck = 0). The right boundary (outlet) is defined as ∇ ⋅ (Dk∇ck)=0. NAPL is located at the bottom of the throat (measuring 0.1 mm × 0.1 mm), with initial solute con- centration of 1 kg/m3. The region between the NAPL and the flowing channel is associated with the stagnant zone, with a length of L, as shown in Fig.  2. Simulations are carried out with the stagnant zones bearing three different lengths: L1 = 0.03  mm, L2 = 0.1  mm, and L3 = 0.18 mm, at each inflow velocity. Thus, in total nine simulations are done for the study in the cavity. The characteristic length, in this case, is equal to the width of the throat Lc = 0.1 mm. The two-phase flow system is first simulated for a short period of 0.01 s, and the system reaches a steady state with the interface configuration matching the given con- tact angle. Then, the simulation of the solute transport and interphase mass transfer begins, and the results are recorded. The mesh size for the model is set at h = 2 × 10–6 m, and the capillary thickness is ε = 2 × 10–6 m. The characteristic mobility is chosen χ = 10. The sim- ulation of interphase mass transfer is terminated after depletion of at least 80% of the initial mass of the species k, and a constant mass transfer coefficient for the system is ensured to be obtained. The domain in Fig. 2 is discretized by 16,088 elements. The time step is 0.1 s, and the total simulations time is 40 s. Due to the simple setup of the cavity, the runtime of each simulation lasts ca. 1 h. 3.2.2 Study of Interphase Mass Transfer in a Porous Medium The solute interphase mass transfer is studied next in a two-dimensional porous medium, for several saturation topologies with different volume fractions of the stagnant zones. It is worth noting that compared to a three-dimensional medium, the major difference is, e.g., pore coordination number, the permeability of the porous medium, and the connectiv- ity of the wetting phase, but physics in the processes remain comparable. The pore space is shaped by circular grains, which have been widely used in micro-model experiments and pore-scale numerical studies (Anna et al., 2014; Ferrari and Lunati, 2013; Godinez- Brizuela et  al., 2017; Liu et  al., 2014; Soulaine et  al., 2021). The size of each grain is chosen randomly from diameters of 225 μm, 250 μm, or 275 μm. The grains are first uni- formly distributed to fill the domain (Fig. 3b), followed by a random shift from its original place in both x- and y-directions. The distances of the random movements of the grains are Fig. 2 Plot of the geometry and boundary conditions of the cavity 650 H. Gao et al. 1 3 controlled to be shorter than 20 μm, so that the grains do not overlap. The 2D geometry is generated at a porosity of 47.56%, measuring 7.36 mm × 3.53 mm in size, as shown in Fig. 3a. The triangular mesh elements have a side length of h = 6 × 10–6 m, and two lay- ers of mesh on the solid wall are further refined to 1.5 × 10–6 m, as shown in Fig. 3c. The mesh consists of in total 1,249,369 elements. The capillary thickness for the model is set at ε = 7 × 10–6 m. The simulation is carried out in two steps: a preparation step to simulate the two- phase flow and obtain a steady state saturation topology, and a main step to simulate the interphase mass transfer. In the preparation step, the main domain with the porous material is initially filled with the NAPL, and an extensive rectangular zone serving as an inlet on the left-hand side is initially filled with water. The extensive rectangular zone (measuring 0.12 mm × 3.53 mm in size) has the same width as the study domain and a small length to save the computational cost. The initial location of the interface is provided in the red dashed line in Fig. 3a, important to ascertain numerical stability. Water is injected at a constant specific flux of 0.1 m/s, 0.05 m/s, 0.04 m/s, and 0.02 m/s from the left boundary, to displace the NAPL. The characteristic mobility is chosen χ = 100, χ = 50, χ = 40 and χ = 20, respectively. The right-hand side of the domain is the outlet with a constant pressure p = 0, and the upper and lower sides of the domain are symmetry (free-slip and non-penetrating) boundaries. The simulations of displacement Fig. 3 Plot of (a) the geometry and boundary conditions of the porous medium (pore space in dark gray color); (b) the initial uniform distribution of circle grains; and (c) finite element mesh generated with trian- gular elements 651Effect of Pore Space Stagnant Zones on Interphase Mass Transfer… 1 3 processes last until steady-state saturation topologies at the given flow rates are obtained. The time step is set at 0.002 s, and the total simulation time 1.2 s. The runt- ime ranges between 100 and 150 h for each simulation of two-phase flow (preparation step). Next, the main step for simulation of the mass transfer starts. At each of the four saturation topologies, three water inflow velocities of 1 × 10−2  m/s, 1 × 10−3  m/s, and 1 × 10−4  m/s are simulated. Thus, in total 12 simulations are done in this study. The inflow velocities in the main step are smaller than those in the preparation step so that the saturation topologies remain unchanged (Karadimitriou et al., 2017). The initial sol- ute concentration is 0 in water and 1  kg/m3 in NAPL. The left boundary is constant concentration ck = 0, and the right boundary is defined as ∇ ⋅ ( Dk∇ck ) =0. The character- istic length, in this case, is equal to the averaged width of the throat in the porous media Lc = 0.062  mm (estimated from distances between centers of the circular grains). The simulation of interphase mass transfer is terminated after depletion of at least 95% of the initial mass of the species k. The time step is set as 1 s, and the total simulation time ranges between 100 and 500 s, depending on different degrees of saturation. The runt- ime is ca. 5 h for each simulation of the interphase mass transfer (main step). 3.3 Data Treatment The macroscopic parameters of the system (average phase concentration, mass transfer rate, interfacial area) can be determined by post-processing the pore-scale simulation results. The average phase concentration of the dissolved solute k can be obtained as: where VNAPL is the region in the domain where Vf,w > 0.5 and Vwater is the region in the domain where Vf,w < 0.5. Since the NAPL in the domain does not flow, the interphase mass transfer rate can be directly obtained by monitoring the changing rate of total solute mass in the NAPL. The total solute mass in the NAPL can be calculated as: and then the interphase mass transfer rate at time step tn can be obtained as: where Δt is the time interval between tn+1 and tn. The diffusive interfacial area, as the region with Vf ,w ∈ [0.05, 0.95] according to Gao et al. (2021a), can be calculated with: where Vinterface is the region in the domain with 0.05 < Vf,w < 0.95. The diffusive inter- facial area ( Adif wn ) can be converted into a sharp interfacial area (Awn) and the correspond- ing specific interfacial area (awn) by: (19)cNAPL = ∫ VNAPL c k dV ∫ VNAPL dV and cwater = ∫ Vwater c k dV ∫ Vwater dV (20)MNAPL = ∫ VNAPL ckdV (21)dM∕dt|tn = MNAPL|tn −MNAPL|tn+1 Δt (22)Adif wn = ∫ Vinterface dV 652 H. Gao et al. 1 3 where bwn is the thickness of the diffusive interface, bwn = 4.1641ε calculated with Eq. (4) at G = 0 according to Gao et al. (2021a), and Vpm is volume of the porous medium. The flowing and stagnant zones are distinguished by the threshold local Péclet number, defined in Sect. 4. 4 Results and Discussion The results of the cavity study are shown in Sect. 4.1, and those investigating the porous medium set up in Sect. 4.2. 4.1 Interphase Mass Transfer in a Cavity First, the simulated velocity profiles in the cavity at inflow velocity uin = 1 × 10–2  m/s at three stagnant zones lengths (L1 = 0.03 mm, L3 = 0.1 mm and L5 = 0.18 mm) are shown as examples in Fig. 4a. The stagnant and flowing zones can be distinguished by a local thresh- old Péclet number. The local Péclet number is calculated as Pel = ‖u‖ ∗ Lc∕Dk , where u is the flow velocity obtained from the simulation results, Lc is a characteristic length equal to the throat width, and Dk is the molecular diffusion coefficient. The local threshold Péclet number can be chosen by analyzing its probability distribution (Aziz et al., 2018). The probabilities are calculated as the ratio of the pore space at the given interval of local Péclet number (e.g., log10(Pel)∈ [−0.1, 0] , log10(Pel)∈ [0, 0.1] , etc.) to the total pore space. The probability distributions of the local Péclet number for the examples in Fig. 4c show multimodal distributions (Hasan et al., 2020). The first large peak observed at the right- hand side indicates the regions of flowing water, and the first small peak from the left-hand side indicates regions of the immobile NAPL, where flow velocity is about four orders of magnitude lower than that in the flowing zones. The small peaks in the middle (for L2 and L3 in Fig. 4c) indicate the stagnant zones, where weak recirculating flow exists with the velocity being about two orders of magnitude smaller than that in the flowing zones. The local threshold Péclet number is chosen as Pel = 10–0.3 (indicated by red dashed line in Fig. 4c), to differentiating the flowing and stagnant water. The distinguished NAPL, stag- nant water, and flowing water are plotted in Fig. 4b. There is almost no stagnant zone when the interface is close to the flowing region (see L1 in Fig. 4b). The stagnant zone can be a recirculation zone of low velocity (see L2 in Fig. 4b), or a combination of a fully stagnant zone and a recirculation zone with a large length (see L3 in Fig. 4b). Then, the mass transfer processes at uin = 1 × 10–4  m/s at stagnant zones lengths of L1 = 0.03  mm and L5 = 0.18  mm are compared and shown as examples in Fig.  5. The average phase solute concentration in both NAPL and water is higher when the length of the stagnant zone (measured between the NAPL-water interface and the flow channel) is larger (Fig. 5a). This is because the larger length of the stagnant zone is, the longer the time of the overall mass transfer of the solute from NAPL into the moving water is, and thus the magnitude of solute mass being flushed by the moving water is lower. For con- venience, we name the depth of path of diffusive solute transport in the stagnant water the diffusion length. The interphase mass transfer rate is plotted in Fig. 5b. At the begin- ning stage for the first 3 s, the mass transfer rates are large, as the solute is transferred (23)Awn = Adif wn bwn and awn = Awn Vpm 653Effect of Pore Space Stagnant Zones on Interphase Mass Transfer… 1 3 from NAPL to the pure water. The mass transfer rates are close for the simulations of L1 and L3, and the mass transfer rates decrease fast. From 3 to 35 s, the mass transfer rate is significantly lower for larger stagnant zone. Because in stagnant water diffusion- dominated transport is much slower compared to advection-controlled transport in the flowing water, a higher solute concentrations develop in the stagnant water region (see Fig. 5c) compared to those in the regions of mobile water. Thus, the concentration gra- dient next to the interface, performing as the driving force, is decreased because of the presence of the stagnant zone. This leads to a lower interphase mass transfer rate. After 35 s, the interphase mass transfer rate of the simulations with stagnant zones (L3) can be Fig. 4 Simulation results of the velocity profiles in the cavity at inflow velocity uin = 1 × 10–2 m/s, with stag- nant zone lengths of L1, L2, and L3: (a) velocity distribution; (b) distinguished mobile and immobile zones with the NAPL indicated by grey color (the white arrows show the direction of the velocity vector, and the length of the arrows is proportional to the natural logarithm of the length of velocity vector); (c) probability distribution of the local Péclet number, with the interval of log10(Pel) = 0.1 654 H. Gao et al. 1 3 slightly larger compared to that without stagnant zones (L1), due to lower solute mass remaining in L1. The solute concentration distributions at t = 1 s, t = 3 s, and t = 8 s are plotted in Fig. 5c. The concentration is plotted in a logarithm so that the solute flux dur- ing mass transfer can be better resolved. It is observed that at t = 1 s, the simulation of L3 is still at the stage of initial diffusive transport in the stagnant water, and for L1 the solute flux is already approaching the region of flowing water because there is almost no stagnant zones. It is also observed that the recirculating flow can promote the solute transport in the stagnant zones (see Fig. 5c at t = 3 s) and affect the established solute concentration profile (see Fig. 5c at t = 8 s). Fig. 5 Simulation results of the mass transfer processes at uin = 1 × 10–4 m/s at stagnant zones lengths of L1 and L5: (a) plot of the phase averaged concentration; (b) plot of the interphase mass transfer rate versus the time; and (c) plot of the solute concentration distribution at t = 1 s, t = 3 s, and t = 8 s 655Effect of Pore Space Stagnant Zones on Interphase Mass Transfer… 1 3 Furthermore, the macroscopic parameters, i.e., average phase concentration, mass trans- fer rates, and interfacial area, are obtained from simulation results as explained in Sect. 3.3. In Fig. 6a–c, mass transfer rates per (sharp) interfacial area are plotted against the driving force (concentration difference as a proxy) for all simulations. The slope of curves indi- cates the corresponding mass transfer coefficients, according to Eq. (1). The crosses in the figure show the times of the beginning of the simulations, the arrival of the constant slope of the data curves, and the end of the simulations. It can be observed that in all cases, the mass transfer coefficients (slopes of the corresponding curves) decrease for the transient mass transfer process at the beginning, and this is due to the rapid drop in the magnitude of the driving force at the beginning. A linear trend can be observed after a short period of time: about 2.4 s for L1, 3.7 s for L2, and 5 s for L3 at each inflow velocity. This means that Fig. 6 Simulation results of interphase mass transfer in a cavity: Plot of the mass transfer rates per interfa- cial area versus the driving force for different flow rates at stagnant zone lengths of (a) L1, (b) L2, and (c) L3, where data points are collected for every 0.1 s. The crosses label the time of the beginning, the arrival of the constant slope, and the end. The dashed lines are the regression lines for the data points in the time intervals given in the legends; and (d) Sherwood number versus the dimensionless diffusion length for all simulations 656 H. Gao et al. 1 3 it takes more time for the system to reach the equilibrium of mass transfer when the length of the stagnant zone is larger. The presence of a stagnant zone induces diffusive transport processes that are slow and therefore equilibration takes longer. When the equilibration is approached, with a linear trend observed, a constant mass transfer coefficient can be obtained for the system, which is in agreement with Graveleau et al., (2017) and Maes and Soulaine (2018). The mass transfer coefficients obtained for the three cases and expressed as a dimensionless Sherwood number (according to Eq. 2) are plotted against the diffusion length (dimensionless by the characteristic length) in Fig. 6d. At each inflow velocity, it is found that the mass transfer coefficient is decreased for larger diffusion lengths, i.e., the position of the interface is at larger distances from the interface to flowing water. Compar- ing the Sh obtained from different inflow velocities at each diffusion length, it can be found that the mass transfer coefficient is larger for a higher water flow rate. The difference in the mass transfer coefficients at different flow rates is significantly decreased with larger diffusion lengths. The flow rate can hardly affect the mass transfer coefficient, because of the solute accumulated in the stagnant zone and the concentration profile established in the stagnant zone. This means the mass transfer process is mainly limited by the diffusive transport in the stagnant zone when the diffusion length is large (e.g., L/Lc > 1.5). Besides, for the simulations with almost no stagnant zones (at L1), the correlation between the Sher- wood number and the Reynolds number for the cavity can be formulated as: with R2 equals 0.99. 4.2 Interphase Mass Transfer in a Porous Medium The velocity and phase distribution patterns for pore-scale two-phase flow simulations are shown in Fig.  7. As expected, it is apparent that the domain with the larger volume of residual NAPL is obtained for lower preparation step injection rates. The phase distribution patterns obtained from the four simulations display NAPL saturations ranging from 0.236 to 0.334, as shown in Fig. 7b. At the same inflow velocity, the velocity in the channels is generally higher for higher residual NAPL saturation (see Fig.  7a, at uin = 1 × 10–2  m/s), because the overall space of the flowing channel in the domain is reduced. Furthermore, the process of determining the mobile and stagnant zones is identical to those described above in Sect. 4.1. The probability distributions of the local Péclet number are bimodal, as shown in Fig. 7c. The first peak on the left-hand side indicating the stagnant zone is larger at higher NAPL saturation, and the second peak on the left-hand side indicates the mobile zones. The volume of the stagnant water increases at higher NAPL saturation, and the flow pathways are more tortuous for higher NAPL saturations ranging between 0.236 and 0.334. It is obtained from the four saturation topologies: Ssw = 0.059 at SNAPL = 0.236, Ssw = 0.103 at SNAPL = 0.284, Ssw = 0.151 at SNAPL = 0.298, and Ssw = 0.205 at SNAPL = 0.334, as shown in Fig. 7. The stagnant water saturation increases with increasing NAPL saturation from SNAPL = 0.236 to SNAPL = 0.334. This is in agreement with Karadimitriou et al. (2017). The specific interfacial areas obtained from all simulations are close, ranging from 1240 (1/m) to 1330 (1/m). Thus, the effects of the interphase mass transfer due to the difference in the interfacial area are minimal in this study. The comparison of phase-averaged concentrations with time for the two examples (at uin = 1 × 10–4  m/s) is shown in Fig.  8a. Averaged concentrations decrease with time for both phases because solute mass transferred to the moving water is being flushed. It is (24)Sh = 2.62Re0.128 657Effect of Pore Space Stagnant Zones on Interphase Mass Transfer… 1 3 also shown that phase averaged concentrations are larger for SNAPL = 0.334 for both phases. This is mainly due to the larger initial solute mass for larger NAPL saturation. Addition- ally, temporal changes in the mass transfer rates are shown in Fig. 8b for the two examples. Mass transfer rates drop rapidly for the first 5 s because of the fast decrease in the driving force. The mass transfer rate is larger for SNAPL = 0.334 after 5 s, and this is mainly due to the upkeep of the NAPL concentration maintaining a stronger driving force for higher NAPL saturation. The concentration distributions at SNAPL = 0.236 and SNAPL = 0.334 are shown in Fig.  8c. At a given time, the concentration is higher in NAPL clusters of larger volume, and lower interfacial area. At t = 40  s, solute concentration is depleted to a quite low value in most NAPL clusters for SNAPL = 0.236 but still relatively high for SNAPL = 0.334. This means that the dissolved solute remains for a relatively long time in NAPL (visible as larger retention times) for larger numbers of multi-pore, NAPL clusters. Besides, the concentration of the solute in water slightly increases in the flow direction. Fig. 7 Simulation results of two-phase flow at steady state: (a) velocity; (b) distinguished mobile and immobile zones, with the NAPL in grey color; (c) probability distribution of the local Péclet numbers, with the interval of log10(Pel) = 0.1 658 H. Gao et al. 1 3 This is in agreement with (Hu et al., 2021). The gradient of the concentration of the solute in the water along the flow direction reduces with time. From Figs. 6 and 7, it is apparent that the stagnant water phase bears higher concentra- tions than the mobile water phase. This means that the decrease in the driving force due to the presence of stagnant zones (discussed in Sect. 4.1) can also be observed in the porous medium. Comparing the phase distribution for SNAPL = 0.236 and SNAPL = 0.334 (Fig. 7), the volume of most of the stagnant water clusters is generally smaller and only one or two menisci are in contact with stagnant water clusters for SNAPL = 0.236, while the volumes Fig. 8 Simulation results of interphase mass transfer processes (at uin = 1 × 10–4 m/s) at SNAPL = 0.236 and SNAPL = 0.334: (a) phase-averaged solute concentration versus time; (b) mass transfer rates versus time (c) concentration distribution patterns at the time of t = 5 s and t = 40 s 659Effect of Pore Space Stagnant Zones on Interphase Mass Transfer… 1 3 of the stagnant water clusters are larger and a larger number of menisci are in contact with stagnant water clusters for SNAPL = 0.334. Besides, the diffusion lengths for an individ- ual meniscus are also different according to the phase distribution for SNAPL = 0.236 and SNAPL = 0.334 (Fig. 7). Because each meniscus has a different length (interfacial area), and each position on the menisci can have different diffusion lengths, it is difficult to quantify the overall diffusion length of the menisci. In this case, we estimate the diffusion length for a meniscus as the shortest straight distance between the center of the meniscus and the connected flowing zone. The estimated diffusion length can be smaller than the real solute diffusion pathway length because the tortuosity of the pathway is not being considered. In Fig. 9a, an example for estimating diffusion lengths in a zoom-in region is shown for the simulation with SNAPL = 0.334. The number of the water/NAPL interfaces exposed to stagnant zones and the average diffusion length of those interfaces in all simulations are plotted against the stagnant water saturation in Fig. 9b. It is found that both the number of interfaces exposed to stagnant zones and the average diffusion lengths of those menisci are generally larger for higher stagnant water saturation. No strictly linear increasing trend is observed because the average diffusion length also depends on the topology of the stagnant water clusters. The stagnant water saturation can generally be considered as a proxy for the overall diffusion length of the solute in the pore space starting at the menisci. Furthermore, the mass transfer rates per interfacial area are plotted versus the driving force for all simulations in Fig. 10. It is observed that the changing of the mass transfer coefficient with respect to time can be mainly separated into two stages. At Stage 1 (for the first ca. 40 s), the mass transfer rates (as the slope of the data curves) decrease continu- ously with time, meaning a kinetic mass transfer process. This is mainly due to the rela- tively faster depletion of the solute in the small NAPL clusters which are mostly in direct contact with the flowing water, as shown in Fig.  8c. Thus, the contribution of the mass transfer happened on the interface associated with small NAPL clusters (with interfaces at Fig. 9 Plot of (a) the estimation procedure for the diffusion length of the interfaces exposed to stagnant water clusters, in an enlarged section of the porous medium for a simulation with SNAPL = 0.334. The inter- faces, labelled 3 to 7, are exposed to the stagnant water zones, and the lengths of the dotted lines repre- sent the estimated diffusion lengths for the individual interfaces (the white arrows show the direction of the velocity vector, and the length of the arrows is proportional to the natural logarithm of the length of velocity vector); and (b) the estimated number of interfaces exposed to the stagnant zones and the averaged diffusion length of these interfaces versus the stagnant zone saturation 660 H. Gao et al. 1 3 small diffusion length) to the overall mass transfer in the domain become less. The mass transfer happening on the interfaces at large diffusion length (mostly in large NAPL clus- ters) becomes dominant, and thus the overall mass transfer coefficient decreases. Though Stage 1 only lasts for about 40 s, it is shown in Fig. 8a that about 70% of the initial sol- ute mass in the domain for SNAPL = 0.236 and 50% for SNAPL = 0.236 is removed from the porous domain at Stage 1. Furthermore, at Stage 2 (after the ca. 40 s), the linear trends are observed for all simulations. The linear trend is approached slightly earlier for lower NAPL saturations and lower volumes of the stagnant zones. It is observed in each subfigure of Fig. 10 that the differences in the flow rates can hardly affect the mass transfer coefficient Fig. 10 Plot of the mass transfer rate per interfacial area versus interphase concentration difference (“driv- ing force”), with the slopes of the curves representing the mass transfer coefficients, and data points are collected for every 1 s (the crosses label the time of injection of 1.4 PVs, 2.7 PVs and 3.9 PVs of water) from simulations of (a) SNAPL = 0.236, (b) SNAPL = 0.284, (c) SNAPL = 0.298, and (d) SNAPL = 0.334. The colored dashed lines represent the regression lines at the beginning of Stage 1 (considering the first five data points), and the black dashed lines represent the regression lines at Stage 2 (for the data points in the time intervals given in the legends) 661Effect of Pore Space Stagnant Zones on Interphase Mass Transfer… 1 3 at Stage 2. This is again because the constant mass transfer coefficients obtained at Stage 2 are mainly associated with the large NAPL clusters (with higher solute concentration remaining) with interfaces at larger diffusion lengths. This again agrees with the findings from the cavity study in Sect. 4.1. Stage 2 lasts for up to 500 s until most of the solute is removed from the porous domain. Finally, the slope of the curves at the beginning of the kinetic mass transfer at Stage 1 and the constant slope of the data curves at Stage 2 are obtained as the mass transfer coef- ficients, which are further expressed as Sherwood numbers (Eq. 2) and plotted against the Reynolds number and the stagnant water saturation in Fig. 11a and Fig. 11b, respectively. It is found that the Sherwood numbers obtained from the beginning of Stage1 (Fig. 11a) are strongly correlated to the Reynolds numbers (flow rates), and the effect of the stagnant water saturation on the Sherwood numbers is small. Because the overall mass transfer at the beginning is dominated by the mass transfer that happened on the interfaces in direct contact with the flowing water. The correlation function is expressed as: On the contrary, it is found that the Sherwood numbers obtained from the Stage 2 is mainly dependent on the stagnant water saturation. For the simulations at the three flow velocities for each saturation, the obtained mass transfer coefficient is the same, as shown in Fig. 10 with the black dashed lines. The Sherwood number (mass transfer coefficient) at Stage 2 significantly decreases with larger stagnant zone saturation, as shown in Fig. 11b. The correlation between the Sherwood number and the stagnant water saturation can be expressed as: The obtained mass transfer coefficients range from 4.26 × 10–6 m/s to 8.36 × 10–6 m/s, which is smaller than the results from (Maes and Soulaine, 2018) that the mass transfer coefficient equal to 1.7 × 10–5 m/s. The difference is mainly because their study focuses on (25)Sh1 = 7.2Re0.11 (26)Sh2 = 2.5e−4.4Ssw Fig. 11 Plot of (a) the Sherwood numbers at the beginning of Stage 1 versus the Reynolds numbers; and (b) the constant Sherwood numbers at Stage 2 versus the stagnant water saturation. The dashed lines represent the correlation functions 662 H. Gao et al. 1 3 the mass transfer process from the mobile phase into the immobile phase, where the effects of the stagnant zones are mitigated. In this section, the NAPL distribution depends on the capillary forces in the porous medium, generated randomly and macroscopically considered homogeneous. Thus, the heterogeneity of the NAPL distribution reported as the main factor for changing the mass transfer coefficient (Agaoglu et  al., 2016; Aminnaji et  al., 2019), is not affected in this study. Therefore, the changes in the Sherwood number are only induced by the differences in the inflow velocity and the differences in saturation topology. The results demonstrate the importance of considering the presence and the characteristics of stagnant zones in the assessment of the interphase mass transfer coefficient. 5 Summary and Conclusions We studied the interphase mass transfer of a dissolved species from a residual immobile NAPL phase into the flowing water phase in a domain with water-saturated flow employing the pore-scale PFM-CST model. The study was carried out in domains with geometries of a single simplified cavity and a porous medium, composed of circular grains. The flowing and stagnant zones were distinguished by applying a threshold local Péclet number. Thus, the regions of flowing water, stagnant water, and NAPL phases could be discriminated and their saturations obtained. Also, the mass transfer rates, the phase-averaged solute concen- trations, and the interfacial areas were quantified. The major results and findings of the study are summarized here: • After a period of time for the transient process at the beginning, the normalized mass transfer rate per interfacial area approaches a linearly increasing trend with increasing interphase concentration difference (“driving force”), and a constant mass transfer coef- ficient can be obtained for the system. This confirms that an equilibrium condition of mass transfer exists for the interphase mass transfer of a dissolved solute under given conditions. • With the same initial solute concentration in the NAPL phase, the individual averaged concentrations in both phases are larger during the mass transfer process when the dif- fusion length associated with the interface is larger, and the retention time of the solute in NAPL is longer. • For the study involving a cavity, a decrease in mass transfer coefficient is found for larger diffusion lengths (larger volumes of stagnant zones) and smaller flow velocity. When the diffusion length is small, the mass transfer coefficient can be well correlated to the inflow velocity. When the diffusion length is large (e.g., L/Lc > 1.5), the mass transfer process is mainly limited by the diffusive transport in the stagnant zone, and changes in the inflow velocity can hardly affect the mass transfer coefficient. • For the study involving a porous medium, it is found that the number of menisci in contact with stagnant zones and the diffusion lengths of these menisci both increase with increasing stagnant zone saturations. The stagnant zone saturation can be used as a proxy for the overall diffusion length of the system. 663Effect of Pore Space Stagnant Zones on Interphase Mass Transfer… 1 3 • In the porous medium, the changing of the mass transfer coefficient can be separated into two stages. For the first stage, the mass transfer coefficient continuously decreases, due to the depletion of the solute in the small NAPL clusters that are in direct contact with the flowing water. The mass transfer coefficient at the beginning of the first stage can be well correlated to the inflow velocity. For the second stage, the mass transfer mainly happens on the interfaces associated with large NAPL clusters with larger diffu- sion lengths. The constant mass transfer coefficients at the second stage are determined by the stagnant zone saturation. This study shows that the stagnant zones play an important role in the interphase mass transfer of a dissolved solute. The stagnant zone saturation is a key macro-scale parameter to be included for developing correlations with the mass transfer coefficient. Future work is required to further investigate how the mass transfer coefficient changes with respect to the factors that determine the stagnant zone saturation, such as wettability and geometry of the porous medium. Appendix 1 The mobility (γ) from Eq.  (6) is an important parameter that affects the accuracy of the results from the simulations using the PFM. The mobility is formulated as a function of the capillary thickness (ε) and a tuning factor (χ): � = ��2 , where χ is called charac- teristic mobility. In Appendix A, we investigate the choices of characteristic mobil- ity and its effects on the results, by simulating a two-phase displacement process at different values of χ range from 1 to 104. The simulation is implemented on a sub- region of the original geometry of porous media (Fig.  12a). The boundary conditions are:u|inlet = 0.01m∕s, p|outlet = 0, ut=0 = 0, pt=0 = 0 . Top and bottom boundaries are the no- flow boundaries, labelled in red lines as shown in Fig. 12a. The small section on left has �t=0 = 1 , the main domain has �t=0 = −1 , and �|inlet = 1, meaning initially NAPL filled domain is drained by water. The channel wall boundary condition is no-slip. The properties of the fluids, the mesh sizes, and the capillary thickness are all the same as that explained in Sect. 3.2. The simulation results are studied by monitoring the diffusive interfacial area, as the region with Vf ,w ∈ [0.05, 0.95] , which can be calculated using Eq. (22). The diffusive interfacial area during the displacement from simulations at χ = 1, χ = 10, χ = 102, χ = 103 and χ = 104 are demonstrated in Fig.  12b. The diffusive interfacial area increases at beginning and starts to stabilized after breakthrough at about 0.06 s. The simu- lated fluid distribution patterns and the monitored diffusive interfacial area are similar for χ ranges from 10 to 103, the effects of choices of χ on the results is minor within this range. It is also observed that when χ ≤ 1, the simulated diffusive interfacial area is much larger than the results with χ ∈[10,103]. This is mainly due to the expansion of the interface thickness, which leads to inaccurate results. On the other hand, when χ ≥ 104, the simulated diffusive interfacial area is much smaller than the results with χ ∈[10,103]. This is mainly because the flow is affected by the strong Cahn–Hilliard diffusion, and the fluid distribution pattern is changed. Thus, the suitable range of characteristic mobility is χ ∈[10,103] for the flow system given in Sect. 3.2. 664 H. Gao et al. 1 3 Acknowledgements We acknowledge the German Research Foundation (Deutsche Forschungsgemein- schaft, DFG), Project Number 428614366, for funding this research. HS and NKK acknowledge funding by the DFG within the Collaborative Research Center 1313 (Project Number 327154368–SFB 1313). HS acknowledges funding by the DFG under Germany’s Excellence Strategy—EXC 2075—390740016. Funding Open Access funding enabled and organized by Projekt DEAL. Data availability All data used to support this work are reported in the manuscript and the supporting infor- mation in the respective tables and figures. Declarations Conflict of interests The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Com- mons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/. Fig. 12 Plot of (a) the geometry of a sub-region of the original porous media and the boundary conditions for study of the choices of mobility; and (b) the region of diffusive interface with V f ,w ∈ [0.05, 0.95] during the two-phase displacement, at different χ http://creativecommons.org/licenses/by/4.0/ 665Effect of Pore Space Stagnant Zones on Interphase Mass Transfer… 1 3 References Agaoglu, B., Copty, N.K., Scheytt, T., Hinkelmann, R.: Interphase mass transfer between fluids in sub- surface formations: A review. Adv. Water Resour. 79, 162–194 (2015). https:// doi. org/ 10. 1016/j. advwa tres. 2015. 02. 009 Agaoglu, B., Scheytt, T., Copty, N.K.: Impact of NAPL architecture on interphase mass transfer: a pore network study. Adv. Water Resour., Pore Scale Model. Exp. 95, 138–151 (2016). https:// doi. org/ 10. 1016/j. advwa tres. 2015. 11. 012 Akhlaghi Amiri, H.A., Hamouda, A.A.: Evaluation of level set and phase field methods in modeling two phase flow with viscosity contrast through dual-permeability porous medium. Int. J. Multiph. Flow 52, 22–34 (2013). https:// doi. org/ 10. 1016/j. ijmul tipha seflow. 2012. 12. 006 Alpak, F.O., Riviere, B., Frank, F.: A phase-field method for the direct simulation of two-phase flows in pore-scale media using a non-equilibrium wetting boundary condition. Comput. Geosci. 20, 881– 908 (2016). https:// doi. org/ 10. 1007/ s10596- 015- 9551-2 Aminnaji, M., Rabbani, A., Niasar, V.J., Babaei, M.: Effects of pore-scale heterogeneity on macroscopic NAPL dissolution efficiency: a two-scale numerical simulation study. Water Resour. Res. 55(11), 8779–8799 (2019). https:// doi. org/ 10. 1029/ 2019W R0260 35 Aydin Sarikurt, D., Gokdemir, C., Copty, N.K.: Sherwood correlation for dissolution of pooled NAPL in porous media. J. Contam. Hydrol. 206, 67–74 (2017). https:// doi. org/ 10. 1016/j. jconh yd. 2017. 10. 001 Aziz, R., Joekar-Niasar, V., Martinez-Ferrer, P.: Pore-scale insights into transport and mixing in steady-state two-phase flow in porous media. Int. J. Multiph. Flow 109, 51–62 (2018). https:// doi. org/ 10. 1016/j. ijmul tipha seflow. 2018. 07. 006 Basirat, F., Yang, Z., Niemi, A.: Pore-scale modeling of wettability effects on CO2–brine displacement dur- ing geological storage. Adv. Water Resour. 109, 181–195 (2017). https:// doi. org/ 10. 1016/j. advwa tres. 2017. 09. 004 Bond, W.J., Wierenga, P.J.: Immobile water during solute transport in unsaturated sand columns. Water Resour. Res. 26, 2475–2481 (1990). https:// doi. org/ 10. 1029/ WR026 i010p 02475 Brusseau, M.L.: Rate-limited mass transfer and transport of organic solutes in porous media that contain immobile immiscible organic liquid. Water Resour. Res. 28, 33–45 (1992). https:// doi. org/ 10. 1029/ 91WR0 2498 Cahn, J.W., Hilliard, J.E.: Free energy of a nonuniform system. III. nucleation in a two-component incom- pressible fluid. J. Chem. Phys. 31, 688–699 (1959). https:// doi. org/ 10. 1063/1. 17304 47 Chomsurin, C., Werth, C.J.: Analysis of pore-scale nonaqueous phase liquid dissolution in etched silicon pore networks. Water Resour. Res. (2003). https:// doi. org/ 10. 1029/ 2002W R0016 43 Coutelieris, F.A., Kainourgiakis, M.E., Stubos, A.K., Kikkinides, E.S., Yortsos, Y.C.: Multiphase mass transport with partitioning and inter-phase transport in porous media. Chem. Eng. Sci. 61, 4650–4661 (2006). https:// doi. org/ 10. 1016/j. ces. 2006. 02. 037 de Anna, P., Jimenez-Martinez, J., Tabuteau, H., Turuban, R., Le Borgne, T., Derrien, M., Méheust, Y.: Mixing and reaction kinetics in porous media: an experimental pore scale quantification. Environ. Sci. Technol. 48, 508–516 (2014). https:// doi. org/ 10. 1021/ es403 105b Deising, D., Marschall, H., Bothe, D.: A unified single-field model framework for Volume-Of-Fluid simu- lations of interfacial species transfer applied to bubbly flows. Chem. Eng. Sci. 139, 173–195 (2016). https:// doi. org/ 10. 1016/j. ces. 2015. 06. 021 Dillard, L.A., Blunt, M.J.: Development of a pore network simulation model to study nonaqueous phase liq- uid dissolution. Water Resour. Res. 36, 439–454 (2000). https:// doi. org/ 10. 1029/ 1999W R9003 01 Dillard, L.A., Essaid, H.I., Blunt, M.J.: A functional relation for field-scale nonaqueous phase liquid dis- solution developed using a pore network model. J. Contam. Hydrol. 48, 89–119 (2001). https:// doi. org/ 10. 1016/ S0169- 7722(00) 00171-6 Dou, Z., Zhang, X., Zhuang, C., Yang, Y., Wang, J., Zhou, Z.: Saturation dependence of mass transfer for solute transport through residual unsaturated porous media. Int. J. Heat Mass Transf. 188, 122595 (2022). https:// doi. org/ 10. 1016/j. ijhea tmass trans fer. 2022. 122595 Essaid, H.I., Bekins, B.A., Cozzarelli, I.M.: Organic contaminant transport and fate in the subsurface: Evo- lution of knowledge and understanding. Water Resour. Res. 51, 4861–4902 (2015). https:// doi. org/ 10. 1002/ 2015W R0171 21 Ferrari, A., Lunati, I.: Direct numerical simulations of interface dynamics to link capillary pressure and total surface energy. Adv. Water Resour. 57, 19–31 (2013). https:// doi. org/ 10. 1016/j. advwa tres. 2013. 03. 005 Gao, H., Tatomir, A.B., Karadimitriou, N.K., Steeb, H., Sauter, M.: A two-phase, pore-scale reactive trans- port model for the kinetic interface-sensitive tracer. Water Resour. Res. (2021a). https:// doi. org/ 10. 1029/ 2020W R0285 72 https://doi.org/10.1016/j.advwatres.2015.02.009 https://doi.org/10.1016/j.advwatres.2015.02.009 https://doi.org/10.1016/j.advwatres.2015.11.012 https://doi.org/10.1016/j.advwatres.2015.11.012 https://doi.org/10.1016/j.ijmultiphaseflow.2012.12.006 https://doi.org/10.1007/s10596-015-9551-2 https://doi.org/10.1029/2019WR026035 https://doi.org/10.1016/j.jconhyd.2017.10.001 https://doi.org/10.1016/j.ijmultiphaseflow.2018.07.006 https://doi.org/10.1016/j.ijmultiphaseflow.2018.07.006 https://doi.org/10.1016/j.advwatres.2017.09.004 https://doi.org/10.1016/j.advwatres.2017.09.004 https://doi.org/10.1029/WR026i010p02475 https://doi.org/10.1029/91WR02498 https://doi.org/10.1029/91WR02498 https://doi.org/10.1063/1.1730447 https://doi.org/10.1029/2002WR001643 https://doi.org/10.1016/j.ces.2006.02.037 https://doi.org/10.1021/es403105b https://doi.org/10.1016/j.ces.2015.06.021 https://doi.org/10.1029/1999WR900301 https://doi.org/10.1016/S0169-7722(00)00171-6 https://doi.org/10.1016/S0169-7722(00)00171-6 https://doi.org/10.1016/j.ijheatmasstransfer.2022.122595 https://doi.org/10.1002/2015WR017121 https://doi.org/10.1002/2015WR017121 https://doi.org/10.1016/j.advwatres.2013.03.005 https://doi.org/10.1029/2020WR028572 https://doi.org/10.1029/2020WR028572 666 H. Gao et al. 1 3 Gao, H., Tatomir, A.B., Karadimitriou, N.K., Steeb, H., Sauter, M.: Effects of surface roughness on the kinetic interface-sensitive tracer transport during drainage processes. Adv. Water Resour. (2021b). https:// doi. org/ 10. 1016/j. advwa tres. 2021. 104044 Godinez-Brizuela, O.E., Karadimitriou, N.K., Joekar-Niasar, V., Shore, C.A., Oostrom, M.: Role of corner interfacial area in uniqueness of capillary pressure-saturation- interfacial area relation under transient conditions. Adv. Water Resour. 107, 10–21 (2017). https:// doi. org/ 10. 1016/j. advwa tres. 2017. 06. 007 Graveleau, M., Soulaine, C., Tchelepi, H.A.: Pore-scale simulation of interphase multicomponent mass transfer for subsurface flow. Transp. Porous. Med. 120, 287–308 (2017). https:// doi. org/ 10. 1007/ s11242- 017- 0921-1 Haroun, Y., Legendre, D., Raynal, L.: Volume of fluid method for interfacial reactive mass transfer: applica- tion to stable liquid film. Chem. Eng. Sci. 65, 2896–2909 (2010a). https:// doi. org/ 10. 1016/j. ces. 2010. 01. 012 Haroun, Y., Legendre, D., Raynal, L.: Direct numerical simulation of reactive absorption in gas–liquid flow on structured packing using interface capturing method. Chem. Eng. Sci. 65(1), 351–356 (2010b). https:// doi. org/ 10. 1016/j. ces. 2009. 07. 018 Haroun, Y., Raynal, L., Legendre, D.: Mass transfer and liquid hold-up determination in structured packing by CFD. Chem. Eng. Sci. 75, 342–348 (2012). https:// doi. org/ 10. 1016/j. ces. 2012. 03. 011 Hasan, S., Joekar-Niasar, V., Karadimitriou, N.K., Sahimi, M.: Saturation dependence of non-fickian trans- port in porous media. Water Resour. Res. 55, 1153–1166 (2019). https:// doi. org/ 10. 1029/ 2018W R0235 54 Hasan, S., Niasar, V., Karadimitriou, N.K., Godinho, J.R.A., Vo, N.T., An, S., Rabbani, A., Steeb, H.: Direct characterization of solute transport in unsaturated porous media using fast X-ray synchrotron microto- mography. PNAS 117, 23443–23449 (2020). https:// doi. org/ 10. 1073/ pnas. 20117 16117 Held, R.J., Celia, M.A.: Pore-scale modeling and upscaling of nonaqueous phase liquid mass transfer. Water Resour. Res. 37, 539–549 (2001). https:// doi. org/ 10. 1029/ 2000W R9002 74 Hu, Y., Patmonoaji, A., Xu, H., Kaito, K., Matsushita, S., Suekane, T.: Pore-scale investigation on nonaque- ous phase liquid dissolution and mass transfer in 2D and 3D porous media. Int. J. Heat Mass Transf. 169, 120901 (2021). https:// doi. org/ 10. 1016/j. ijhea tmass trans fer. 2021. 120901 Imhoff, P.T., Jaffé, P.R., Pinder, G.F.: An experimental study of complete dissolution of a nonaqueous phase liquid in saturated porous media. Water Resour. Res. 30, 307–320 (1994). https:// doi. org/ 10. 1029/ 93WR0 2675 Jacqmin, D.: Calculation of two-phase navier-stokes flows using phase-field modeling. J. Comput. Phys. 155, 96–127 (1999). https:// doi. org/ 10. 1006/ jcph. 1999. 6332 Jia, C., Shing, K., Yortsos, Y.C.: Visualization and simulation of non-aqueous phase liquids solubilization in pore networks. J. Contam. Hydrol. 35, 363–387 (1999). https:// doi. org/ 10. 1016/ S0169- 7722(98) 00102-8 Karadimitriou, N.K., Joekar-Niasar, V., Babaei, M., Shore, C.A.: Critical role of the immobile zone in non- fickian two-phase transport: a new paradigm. Environ. Sci. Technol. 50, 4384–4392 (2016). https:// doi. org/ 10. 1021/ acs. est. 5b059 47 Karadimitriou, N.K., Joekar-Niasar, V., Brizuela, O.G.: Hydro-dynamic solute transport under two-phase flow conditions. Sci. Rep. 7, 1–7 (2017). https:// doi. org/ 10. 1038/ s41598- 017- 06748-1 Kechagia, P.E., Tsimpanogiannis, I.N., Yortsos, Y.C., Lichtner, P.C.: On the upscaling of reaction-transport processes in porous media with fast or finite kinetics. Chem. Eng. Sci. 57, 2565–2577 (2002). https:// doi. org/ 10. 1016/ S0009- 2509(02) 00124-0 Kennedy, C.A., Lennox, W.C.: A pore-scale investigation of mass transport from dissolving DNAPL drop- lets. J. Contam. Hydrol. 24, 221–246 (1997). https:// doi. org/ 10. 1016/ S0169- 7722(96) 00011-3 Knutson, C.E., Werth, C.J., Valocchi, A.J.: Pore-scale modeling of dissolution from variably distributed nonaqueous phase liquid blobs. Water Resour. Res. 37, 2951–2963 (2001). https:// doi. org/ 10. 1029/ 2001W R0005 87 Liu, H., Valocchi, A.J., Werth, C., Kang, Q., Oostrom, M.: Pore-scale simulation of liquid CO2 displace- ment of water using a two-phase lattice Boltzmann model. Adv. Water Resour. 73, 144–158 (2014). https:// doi. org/ 10. 1016/j. advwa tres. 2014. 07. 010 Maes, J., Soulaine, C.: A new compressive scheme to simulate species transfer across fluid interfaces using the Volume-Of-Fluid method. Chem. Eng. Sci. 190, 405–418 (2018). https:// doi. org/ 10. 1016/j. ces. 2018. 06. 026 Maes, J., Soulaine, C.: A unified single-field Volume-of-Fluid-based formulation for multi-component inter- facial transfer with local volume changes. J. Comput. Phys. 402, 109024 (2020). https:// doi. org/ 10. 1016/j. jcp. 2019. 109024 Mahani, H., Berg, S., Ilic, D., Bartels, W.-B.-B., Joekar-Niasar, V.: Kinetics of low-salinity-flooding effect. SPE J. 20, 8–20 (2014). https:// doi. org/ 10. 2118/ 165255- PA https://doi.org/10.1016/j.advwatres.2021.104044 https://doi.org/10.1016/j.advwatres.2017.06.007 https://doi.org/10.1007/s11242-017-0921-1 https://doi.org/10.1007/s11242-017-0921-1 https://doi.org/10.1016/j.ces.2010.01.012 https://doi.org/10.1016/j.ces.2010.01.012 https://doi.org/10.1016/j.ces.2009.07.018 https://doi.org/10.1016/j.ces.2012.03.011 https://doi.org/10.1029/2018WR023554 https://doi.org/10.1029/2018WR023554 https://doi.org/10.1073/pnas.2011716117 https://doi.org/10.1029/2000WR900274 https://doi.org/10.1016/j.ijheatmasstransfer.2021.120901 https://doi.org/10.1029/93WR02675 https://doi.org/10.1029/93WR02675 https://doi.org/10.1006/jcph.1999.6332 https://doi.org/10.1016/S0169-7722(98)00102-8 https://doi.org/10.1016/S0169-7722(98)00102-8 https://doi.org/10.1021/acs.est.5b05947 https://doi.org/10.1021/acs.est.5b05947 https://doi.org/10.1038/s41598-017-06748-1 https://doi.org/10.1016/S0009-2509(02)00124-0 https://doi.org/10.1016/S0009-2509(02)00124-0 https://doi.org/10.1016/S0169-7722(96)00011-3 https://doi.org/10.1029/2001WR000587 https://doi.org/10.1029/2001WR000587 https://doi.org/10.1016/j.advwatres.2014.07.010 https://doi.org/10.1016/j.ces.2018.06.026 https://doi.org/10.1016/j.ces.2018.06.026 https://doi.org/10.1016/j.jcp.2019.109024 https://doi.org/10.1016/j.jcp.2019.109024 https://doi.org/10.2118/165255-PA 667Effect of Pore Space Stagnant Zones on Interphase Mass Transfer… 1 3 Marschall, H., Hinterberger, K., Schüler, C., Habla, F., Hinrichsen, O.: Numerical simulation of species transfer across fluid interfaces in free-surface flows using OpenFOAM. Chem. Eng. Sci. 78, 111–127 (2012). https:// doi. org/ 10. 1016/j. ces. 2012. 02. 034 Miller, C.T., Poirier-McNeil, M.M., Mayer, A.S.: Dissolution of trapped nonaqueous phase liquids: mass transfer characteristics. Water Resour. Res. 26, 2783–2796 (1990). https:// doi. org/ 10. 1029/ WR026 i011p 02783 Powers, S.E., Abriola, L.M., Weber, W.J., Jr.: An experimental investigation of nonaqueous phase liquid dissolution in saturated subsurface systems: steady state mass transfer rates. Water Resour. Res. 28, 2691–2705 (1992). https:// doi. org/ 10. 1029/ 92WR0 0984 Powers, S.E., Abriola, L.M., Dunkin, J.S., Weber, W.J.: Phenomenological models for transient NAPL- water mass-transfer processes. J. Contam. Hydrol. 16, 1–33 (1994a). https:// doi. org/ 10. 1016/ 0169- 7722(94) 90070-1 Powers, S.E., Abriola, L.M., Weber, W.J., Jr.: An experimental investigation of nonaqueous phase liquid dis- solution in saturated subsurface systems: transient mass transfer rates. Water Resour. Res. 30, 321–332 (1994b). https:// doi. org/ 10. 1029/ 93WR0 2923 Quintard, M., Whitaker, S.: Convection, dispersion, and interfacial transport of contaminants: Homogene- ous porous media. Adv. Water Resour. 17, 221–239 (1994). https:// doi. org/ 10. 1016/ 0309- 1708(94) 90002-7 Schaffer, M., Maier, F., Licha, T., Sauter, M.: A new generation of tracers for the characterization of inter- facial areas during supercritical carbon dioxide injections into deep saline aquifers: Kinetic interface- sensitive tracers (KIS tracer). Int. J. Greenhouse Gas Control 14, 200–208 (2013). https:// doi. org/ 10. 1016/j. ijggc. 2013. 01. 020 Smedt, F.D., Wierenga, P.J.: Solute transfer through columns of glass beads. Water Resour. Res. 20, 225– 232 (1984). https:// doi. org/ 10. 1029/ WR020 i002p 00225 Soulaine, C., Debenest, G., Quintard, M.: Upscaling multi-component two-phase flow in porous media with partitioning coefficient. Chem. Eng. Sci. 66, 6180–6192 (2011). https:// doi. org/ 10. 1016/j. ces. 2011. 08. 053 Soulaine, C., Maes, J., Roman, S.: Computational Microfluidics for Geosciences. Front. Water 3 (2021) Tatomir, A.B., Jyoti, A., Sauter, M.: Monitoring of CO2 Plume Migration in Deep Saline Formations with Kinetic Interface Sensitive Tracers (A Numerical Modelling Study for the Laboratory). In: Vishal, V., Singh, T.N. (eds.) Geologic Carbon Sequestration: Understanding Reservoir Behavior, pp. 59–80. Springer International Publishing, Cham (2016) van Genuchten, M.T., Wierenga, P.J.: Mass transfer studies in sorbing porous media i. analytical solutions 1. Soil Sci. Soc. Am. J. 40, 473–480 (1976). https:// doi. org/ 10. 2136/ sssaj 1976. 03615 99500 40000 40011x Weber, W.J., McGinley, P.M., Katz, L.E.: Sorption phenomena in subsurface systems: concepts, models and effects on contaminant fate and transport. Water Res. 25, 499–528 (1991). https:// doi. org/ 10. 1016/ 0043- 1354(91) 90125-A Yang, L., Nieves-Remacha, M.J., Jensen, K.F.: Simulations and analysis of multiphase transport and reac- tion in segmented flow microreactors. Chemical Engineering Science: Special Issue on Multiphase Microfluidic Engineering; Edited by Krishna Nigam, Volker Hessel, Chun-Xia Zhao and Anton Mid- delberg; and: Special issue of selected papers from CFD in the minerals and process industries; Edited by David F. Fletcher, Niels G. Deen, Petar Liovic, M. Philip Schwarz, Peter J. Witt 169, 106–116 (2017). https:// doi. org/ 10. 1016/j. ces. 2016. 12. 003 Yue, P., Feng, J.J., Liu, C., Shen, J.: A diffuse-interface method for simulating two-phase flows of complex fluids. J. Fluid Mech. 515, 293–317 (2004). https:// doi. org/ 10. 1017/ S0022 11200 40003 70 Yue, P., Zhou, C., Feng, J.J., Ollivier-Gooch, C.F., Hu, H.H.: Phase-field simulations of interfacial dynamics in viscoelastic fluids using finite elements with adaptive meshing. J. Comput. Phys. 219, 47–67 (2006). https:// doi. org/ 10. 1016/j. jcp. 2006. 03. 016 Zhang, C., Werth, C.J., Webb, A.G.: A magnetic resonance imaging study of dense nonaqueous phase liquid dissolution from angular porous media. Environ. Sci. Technol. 36, 3310–3317 (2002). https:// doi. org/ 10. 1021/ es011 497v Zhao, W., Ioannidis, M.A.: Pore network simulation of the dissolution of a single-component wetting non- aqueous phase liquid. Water Resour. Res. (2003). https:// doi. org/ 10. 1029/ 2002W R0018 61 Publisher’s Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. https://doi.org/10.1016/j.ces.2012.02.034 https://doi.org/10.1029/WR026i011p02783 https://doi.org/10.1029/WR026i011p02783 https://doi.org/10.1029/92WR00984 https://doi.org/10.1016/0169-7722(94)90070-1 https://doi.org/10.1016/0169-7722(94)90070-1 https://doi.org/10.1029/93WR02923 https://doi.org/10.1016/0309-1708(94)90002-7 https://doi.org/10.1016/0309-1708(94)90002-7 https://doi.org/10.1016/j.ijggc.2013.01.020 https://doi.org/10.1016/j.ijggc.2013.01.020 https://doi.org/10.1029/WR020i002p00225 https://doi.org/10.1016/j.ces.2011.08.053 https://doi.org/10.1016/j.ces.2011.08.053 https://doi.org/10.2136/sssaj1976.03615995004000040011x https://doi.org/10.1016/0043-1354(91)90125-A https://doi.org/10.1016/0043-1354(91)90125-A https://doi.org/10.1016/j.ces.2016.12.003 https://doi.org/10.1017/S0022112004000370 https://doi.org/10.1016/j.jcp.2006.03.016 https://doi.org/10.1021/es011497v https://doi.org/10.1021/es011497v https://doi.org/10.1029/2002WR001861 Effect of Pore Space Stagnant Zones on Interphase Mass Transfer in Porous Media, for Two-Phase Flow Conditions Abstract Article highlights 1 Introduction 1.1 Interphase Mass Transfer 1.2 Effects of Stagnant Zones 1.3 Numerical Methods at Pore Scale 1.4 Objectives 2 Mathematical Model 2.1 Phase-Field Method 2.2 Two-Phase Flow Dynamics 2.3 Solute Transport 3 Numerical Details 3.1 Numerical Implementation 3.2 Model Setup 3.2.1 Study of Interphase Mass Transfer in a Cavity 3.2.2 Study of Interphase Mass Transfer in a Porous Medium 3.3 Data Treatment 4 Results and Discussion 4.1 Interphase Mass Transfer in a Cavity 4.2 Interphase Mass Transfer in a Porous Medium 5 Summary and Conclusions Acknowledgements References