A Two‐Phase, Pore‐Scale Reactive Transport Model for the Kinetic Interface‐Sensitive Tracer

Previous laboratory experiments with the kinetic interface sensitive (KIS) tracers have shown promising results with respect to the quantification of the fluid‐fluid interfacial area (IFA) under dynamic, two‐phase flow conditions. However, pore‐scale effects relevant to two‐phase flow (e.g., the formation of hydrodynamically stagnant/immobile zones) are not yet fully understood, and quantitative information about how far these effects influence the transport of the tracer reaction products is not yet available. Therefore, a pore‐scale numerical model that includes two‐phase, reactive flow, and transport of the KIS tracer at the fluid‐fluid interface is developed. We propose a new method to quantitatively analyze how the mass of the KIS‐tracer reaction product in the flowing water is affected by the presence of the immobile zones. The model employs the phase field method (PFM) and a new continuous mass transfer formulation, consistent with the PFM. We verify the model with the analytical solutions of transport involving advection, reaction and diffusion processes. The model is tested for two‐phase flow conditions in a conceptual 2D slit. The applicability of the model is demonstrated in NAPL/water drainage scenarios in a conceptual porous domain, comparing the results in terms of the spatial distribution of the phases and solute concentration. Furthermore, we distinguish the mobile and immobile zones based on the local Péclet number, and the corresponding IFA, and solute mass in these two zones is quantified. Finally, we show that the solute mass in flowing water can be employed to selectively determine the mobile part of the IFA.

dissolves in nonpolar liquids, such as supercritical CO 2 , or its analog liquid-phase n-octane (see Schaffer et al., 2013). Thus, a KIS tracer is dissolved in a nonpolar liquid usually acting as a nonwetting phase. With the KIS tracer dissolved, the nonpolar liquid is injected into the initially water saturated porous medium and the KIS tracers are adsorbed onto the fluid-fluid interface (Figure 1). The adsorption process is assumed to follow Langmuir's isotherm (Equation 1): where K L is the Langmuir adsorption coefficient of the (nonwetting phase) dissolved KIS tracer at the fluid-fluid interface, c a the concentration of the KIS tracer in the (bulk) nonwetting phase, c i the KIS concentration at the interface, usually assumed to be approximately equal to the maximum concentration of the saturated interface c i,max . The adsorbed 2-NSAPh molecules at the fluid-fluid interface, when getting in contact with the water molecules, undergo an irreversible hydrolysis reaction. (2) Due to the excessive supply of the tracer in the nonwetting phase (high c a ), its concentration at the interface c i can be assumed to remain constant at c i,max (Schaffer et al., 2013;Tatomir et al., 2018). Thus, the originally first order hydrolysis reaction can be simplified into a pseudo zero-order reaction. The measured FIFA can be expressed with Equation 3: where  c nw w R (mol·s −1 m −2 ) is the reaction rate of the KIS tracer per unit interfacial area, (mol) is the mass of the reaction product k, Δt is the reaction time, and wn A is the FIFA (m 2 ). In a series of batch experiments  c nw w R was determined at 1.04 × 10 −11 mol·s −1 m −2 (Tatomir et al., 2018). From Equation 3, it is obvious that the mass transfer rate is only controlled by the magnitude of the FIFA. The two reaction products are naphthalene-2-sulfonic acid (2-NSA) and phenol, both measurable in the water phase samples. Tatomir et al. (2018) and Schaffer et al. (2013) showed in static batch experiments with the FIFA kept constant, that the 2-NSA is the compound more easily measured in a tracer experiment because of its enhanced fluorescence (2-NSA concentration in the order of μgL −1 ). With 2-NSA characterized by high polarity, that is, highly hydrophilic, it will be distributed into the bulk water (here also the wetting phase) away from the interface. The octanol/water partition coefficient of 2-NSA is equal to logD ow = −2.87 at PH > 5, as stated in the work of (Schaffer et al., 2013). The high water solubility explains a negligible back-partitioning. Static batch experiments showed a linear increase of the 2-NSA concentration in water (Tatomir et al., 2018). Furthermore, the adsorption of 2-NSA on water-solid interfaces has not been observed to have any visible effect on mass transport in previous laboratory studies, and thus it is not considered further (Schaffer et al., 2013;Tatomir et al., 2018). With a zero-order reaction at the interface, numerical modeling of 2-NSAPh decay and transport in the nonwetting phase becomes redundant, and the most relevant component of mass transport modeling is the reaction product 2-NSA, present only in the wetting phase. On the pore-scale, modeling of the hydrolysis reaction and the 2-NSA transport imply several main steps: (1) zero-order reaction, which is equivalent to a constant production of 2-NSA at the interface, (2) distribution of produced 2-NSA in the water phase controlled by molecular diffusion and high water solubility, and (3)  The potential of applying KIS tracers in real porous media to determine the FIFA under dynamic conditions was first demonstrated by Tatomir et al. (2018). Tatomir et al. (2018) provided a proof-of-concept using controlled column experiments with well-characterized porous media composed of glass beads. The KIS tracers were dissolved in a NAPL (n-octane) to displace the water in, initially, fully saturated columns. At the column outlet measurements of fluid volume and 2-NSA concentration of the samples collected at defined time intervals provided breakthrough curves (BTCs) of fluid volumes and the 2-NSA concentration in the water phase. The analysis of the experimental data was conducted employing a macro-scale, reactive, two-phase flow and transport model (Tatomir et al., 2015(Tatomir et al., , 2018(Tatomir et al., , 2019(Tatomir et al., , 2020Tatomir, Jyoti, et al., 2016), with the interfacial reaction of the tracer being specified in the zones with both fluids co-existing, and the FIFA approximated explicitly as a function of saturation. The model predicted a nearly linear concentration BTC which could fit most of the experimental data points. The resulting specific FIFA (defined as the interfacial area per unit volume of porous medium) ranged between 500 and 750m −1 for glass beads with a mean diameter of 240 μm. Despite the successful application of the KIS tracers in laboratory column experiments, new questions came up, such as: (a) what do the FIFA tracers actually measure, that is, capillarity-associated interfacial area (terminal menisci), film-associated interfacial area or total interfacial area, (b) how much 2-NSA will end up in the water-film coating the grains and in the hydro-dynamically stagnant zones, and (c) how do these stagnant zones influence the resulting solute mass in the flowing water?
The stagnant zones, also referred to as the immobile zones, for the porous media flow are defined as the region where flow velocity is very weak and mass transport becomes diffusive (Aziz et al., 2019;Karadimitriou et al., 2016;van Genuchten & Wierenga, 1976). The experimental study at column scale found that the presence of immobile zones can affect the BTCs of the tracer in two-phase flow systems, that is, larger longitudinal dispersion, early breakthrough and long tailings (Bond & Wierenga, 1990;Bromly & Hinz, 2004;Khan & Jury, 1990;Smedt & Wierenga, 1984). The pore-scale micro-model experiments studying the development, role and importance of immobile zones was first conducted by Karadimitriou et al. (2016), who found a nonlinear contribution of the immobile zone to the dispersion coefficient and a nonmonotonic relation between immobile zone saturation and total saturation under transient transport conditions. Karadimitriou et al. (2017) further studied the impact of the Péclet number (ratio of the advective to the diffusive transport rates) under several given saturation topologies, and they found that the ratio of immobile zone saturation to total saturation is not influenced by the different flow rates, if the tortuosity remains identical. By pore-scale direct simulation with the volume of fluid method, Aziz et al. (2018) found the hydrodynamic transport can be non-Fickian in a homogeneous porous medium, and they proposed that relative permeability is possibly a proxy for the stagnant saturation. Hasan et al. (2019) confirmed these findings with pore-network modeling (PNM), and they found a link between the immobile zone saturation and the relative permeability. The flaws of the existing macroscopic models to describe the non-Fickian transport in two-phase flow process were proposed in these recent studies (Aziz et al., 2018;Hasan et al., 2019;Karadimitriou et al., 2016Karadimitriou et al., , 2017. In contrast to these studies with a conservative tracer injected simultaneously with the invading phase, the KIS tracer dissolved in the invading phase depends on the hydrolysis reaction at the interface and the detection of the by-product in the receding phase. This means, in our case, the effect of the immobile zones is manifested in the solute transport in the receding phase. The above questions, regarding the effect of the immobile zones on the distribution of the reacted solute in KIS-TT, can probably be answered by investigating KIS tracer transport mechanisms at the pore-scale. Experimentally, the pore-scale study of KIS tracer reactive transport in two-phase flow can be carried out in micro-models (Karadimitriou et al., 2016), considering its good fluorescent properties. The problem can also be studied by pore-scale modeling. Pore-scale modeling plays an important role in understanding the porescale phenomena and can provide fundamental insight to understand the macro-scale processes (Meakin & Tartakovsky, 2009). With pore-scale modeling, the flow and transport properties of the soil/rock matrix are not averaged but are directly resolved, and the pore-space geometry is explicitly represented either by using idealized geometries or by reconstructed geometries based on XMT images (e.g., Culligan et al., 2004;Peche et al., 2016;Tatomir, Halisch, et al., 2016). In addition, running pore-scale simulations is a more flexible approach than the experiments since there is freedom to tune accordingly the fluids' physical and chemical properties, and the geometry of the porous medium for a variety of boundary conditions. Considering all the above, the main objective of this study is to develop a pore-scale model that can simulate the behavior GAO ET AL.

10.1029/2020WR028572
3 of 25 of the reactive KIS tracer transport under two-phase flow conditions. The following sections review the relevant literature on pore-scale models including reactive transport.

Pore-Scale Numerical Methods
Pore-scale modeling methods mainly include PNM and direct numerical simulation (DNS) approaches. The DNS approaches include particle-based methods, for example, the Lattice Boltzmann method (LBM), smooth particle hydrodynamics, and grid-based computational fluid dynamics (CFD), such as the level-set method (LSM), the volume-of-fluid method (VOF), and the phase-field method (PFM) (Alpak et al., 2016;Meakin & Tartakovsky, 2009). PNM is a well-developed method for pore-scale studies (Blunt & King, 1990;Hasan et al., 2019;Joekar-Niasar et al., 2008;Raoof et al., 2013). PNM simplifies porous media into networks of pores and throats, where flow is governed by Poiseuille's law (Joekar-Niasar et al., 2008). PNM is computationally more efficient than DNS, which provides the full coupling between capillary and viscous forces (Alpak et al., 2016). With PNM it is computationally cheaper to consider domains with a size large enough to be considered as representative elementary volumes (REV), and thus the model can be applied for evaluation of continuum-scale problems. Despite these advantages, PNM is limited by its basis on simplified physics and simplified representations of the rock (Alpak et al., 2016;Basirat et al., 2017;Meakin & Tartakovsky, 2009;Yin et al., 2019). Thus, the DNS methods, based on first principles, are better in capturing transport phenomena and fluid dynamics on the micro-scale, in (natural) pore space with complex geometries (Alpak et al., 2016;Meakin & Tartakovsky, 2009). The LBM is a popular method for pore-scale studies of multiphase reactive transport, because of its advantages in several aspects, such as being able to deal with complex boundaries, to incorporate microscopic fluid-fluid and fluid-solid interactions, and being able to implement a parallelization of the algorithm (Kang et al., 2006;Liu et al., 2014Liu et al., , 2015. The LBM model solves the LB equations for fluid flow and solute transport. In LB simulations, the chemical reaction of species is treated as a homogenous reaction in the bulk fluid and a heterogeneous reaction at the interface as a kinetic boundary condition. The LBM modeling approach has been used in applications such as geological storage of CO 2 , which involves precipitation and dissolution at the fluid-solid surface (Kang et al., 2010) and CO 2 dissolution trapping (Chen et al., 2018). However, one of the disadvantages of the LBM models is that the relation between interaction forces and fluid dynamics requires complex calibration procedures, with many adjustment parameters, such as the adequate approximation of a specific physical system (Ferrari & Lunati, 2013;Frank et al., 2018). With the CFD approach, the Navier-Stokes equations are directly solved in a discretized domain by finite volume or finite element techniques, and the interface between two fluids is represented by an indicator function, such as the volume fraction in VOF and the phase variable in PFM. One challenge for the simulation of multiphase reactive transport with the CFD method is the handling of the concentration jump at the interface when a solute solubility is different in the fluids on either side of the interface (Maes & Soulaine, 2018). Haroun et al. (2010) managed to tackle this challenge with a new Continuous Species Transfer (CST) formulation developed in the VOF framework. In the CST formulation, a constant partition coefficient (or Henry's constant) is introduced to solve the thermodynamic equilibrium of the solute at the interface between two fluids (Haroun et al., 2010). The approach of Haroun et al. (2010) allows for the modeling of the discontinuous solute concentration across the interface, while respecting the continuity of the local diffusive solute mass flux. Graveleau et al. (2017) applied the VOF-CST model to simulate subsurface flow problems with moving contact lines. Maes and Soulaine (2018) identified that the CST generates large numerical diffusion in the phase concentration, which leads to inaccurate simulation of solute mass partitioning between two phases. Therefore, they proposed a new approach, termed compressive CST (C-CST), by adding a compressive term. The C-CST formula managed to significantly reduce numerical errors. However, one major disadvantage of the VOF method is that, since the volume fraction is a step function, accurate curvature and smooth physical quantities near the interface are hard to be obtained (Alpak et al., 2016;Sun & Tao, 2010). Another limitation of the VOF method is that a solid wall boundary is implemented indirectly with an additional moving contact line model (Basirat et al., 2017;Meakin & Tartakovsky, 2009).
In contrast to the VOF method, LSM and PFM with a smooth indicator function are better to treat the curvature and the physical quantities at the interface. LSM has been reported to generate more numerical errors when the interface experiences severe stretching or tearing, and the mass is not conserved (Sun & Tao, 2010; GAO ET AL.
10.1029/2020WR028572 4 of 25 Sussman & Puckett, 2000). Major advantages of the PFM are both the mass conservation and the ability to compute accurately the curvature at the interface (Akhlaghi Amiri & Hamouda, 2013). Furthermore, PFM treats the interface thermodynamically as a diffusive thin layer formed by the mixture of the two fluid phases, and the dynamics of the diffusive interface is governed by the free energy it contains (Yue et al., 2006). Thus, the PFM is physically more consistent than the VOF or the LSM models (Alpak et al., 2016). Another advantage is that considering a diffusive interface rather than a sharp interface around the contact line results effectively in slip, through the diffusive fluxes between the bulk fluids (Ding & Spelt, 2007). Besides, with a diffusive interface in PFM, the zero-order reaction of 2-NSA at fluid-fluid interface can be directly implemented as a homogenous reaction in the region of the diffusive interface, where the two fluids exhibit some limited mixing. The PFM, solving the coupled Cahn-Hilliard and Navier-Stokes equations, is already well developed and numerous simulation studies were performed addressing the two-phase flow system (Akhlaghi Amiri & Hamouda, 2013;Alpak et al., 2016;Jacqmin, 1999;Yue et al., 2004Yue et al., , 2006. In addition, numerous recent studies have applied PFM to study various subsurface flow problems related to fluid viscosity, capillarity, temperature, wettability, heterogeneity, and fractures (Akhlaghi Amiri & Hamouda, 2013Basirat et al., 2017;Rokhforouz & Akhlaghi Amiri, 2017. To our knowledge, PFM has not yet been applied to study reactive transport for two phase flow in porous media.
Here, we derive a new CST formulation in a PFM framework, and therefore the model is termed PFM-CST. At first, we verify the developed PFM-CST model, and, we consequently employ the code to study the KIS tracer reaction and transport processes in immiscible two-phase flow in porous media. We focus on the understanding of how the fluid immobile zones affect the magnitude of the total transported solute mass. The paper is organized as follows: Section 2 introduces the mathematical and numerical models. In Section 3, we present the details of numerical setup for model verifications and its application to a realistic pore space geometry to study the KIS tracer reaction and transport of its reaction by-product 2-NSA. The result of the model verifications and application are shown in Section 4. Section 5 lists the main conclusions.
In this section, we briefly review the basic underlying theory implemented in the model. Instead of a sharp interface between two fluids, the PFM treats the interface as a thin diffusive layer formed by the mixture of two fluids. The PFM introduces a smoothly changing phase variable  to describe the composition of the fluid mixture at the interface, and  remains constant in the bulk fluid phases. The PFM is based on the free energy density as a function of the phase variable and its gradient (Jacqmin, 1999). The free energy density (f mix ) (J/m 3 ) for isothermal mixing of two fluids can be expressed in the Ginzburg-Landau form (Equation 4) (Yue et al., 2006): where λ (N) is the magnitude of the mixing energy and ε (m) is a capillary width representing the thickness of the diffusive interface. The free energy density is made up of two components (on the right side of Equation 4): the first term accounts for the surface energy and the second term is the bulk energy. The surface energy represents the interaction between two fluids and the preference of two fluids for complete mixing (Yue et al., 2004). The bulk energy describes the two-phase immiscibility and involves two minimal values, where  = −1 and  = 1, representing the two bulk phases (wetting and nonwetting fluid), respectively (Jacqmin, 1999). These two competing energies determine the  profile across the interface (Jacqmin, 1999;Yue et al., 2004). The chemical potential G (J/m 3 ) can be defined as the variation of the free energy with respect to the dimensionless phase variable (Equation 5): where Ω is the region of space occupied by the two fluid phases. Considering the interface in one-dimension with boundary conditions: The interfacial tension σ (N/m) is the total free energy at the interface per unit area of the interface, and the relationship between the interfacial tension, the capillary width, and the mixing energy density can be obtained at equilibrium from Equation 7 (Yue et al., 2004): Fluid mass conservation is governed by the Cahn-Hilliard equation (Equation 8), which assumes that the diffusive fluid flux is proportional to the gradient of the chemical potential (Cahn & Hilliard, 1959): where γ (m 3 s/kg) is the mobility expressed as a function of the interface thickness and a tuning factor: The tuning factor χ (m·s/kg) is called the characteristic mobility governing the relaxation time of the interface (Akhlaghi Amiri & Hamouda, 2013). χ needs to be large enough to maintain a constant thickness of the interface, but small enough not to dampen the flow (Jacqmin, 1999).

Two Phase Flow Dynamics
Momentum conservation for an incompressible fluid is governed by the Navier-Stokes equation (Equation 9): where I is the identical tensor, p is pressure, u is the fluid velocity vector, the surface tension is considered as a body force (Yue et al., 2006), and ρ and μ are respectively the density and viscosity of the mixture (nonwetting fluid and wetting fluid) according to volume fraction of the fluids: where subscripts w and nw indicate the wetting and nonwetting phase, respectively. ρ and μ are independent of the concentration of KIS tracer and its reacted product, because the tracer and its reacted product have negligible mass fraction in the fluids. The volume fraction of the fluids can be obtained as: At the grain surfaces, a no-slip boundary is applied, which implies that u = 0 in Equation 9 for solid wall boundaries. The wetting condition on the solid wall is expressed by Equation 14: GAO ET AL.
10.1029/2020WR028572 6 of 25 where n is the (outward) normal vector to the wall and  w is the contact angle.

Solute Advection
To obtain the solute mass conservation equation, we first write the fluid mass conservation equation with the volume fraction of phases. For the wetting phase, by substituting    Due to   0 u , by dividing the constant coefficient one obtains: ) into the Cahn-Hilliard equation (Equation 8), one obtains the volume fraction based mass conservation equation for the nonwetting phase. Thus, the total fluid volumetric flux for the wetting phase F w and nonwetting phase F nw can be expressed as: Hence, for the solute k dissolving in fluid phase α (α = w, nw) with the concentration of , a k c , the solute mass conservation equation writes: where only the solute advection is included here, and the other processes (reaction, diffusion, solubility) will be discussed and added in next sections. Solute advection is the transport process due to the bulk fluid motion. The bulk fluid motion consists of two parts: the advective flux and the diffusive flux at the interface due to the gradient of the chemical potential. In this case, the advection of the solute is ensured to be consistent with the fluid motion. The advection term applied here is similar to that in the C-CST model developed by Maes and Soulaine (2018), where the compressive term for VOF method is replaced by the interface diffusion term for PFM.

Reaction at the Interface
A zero-order interfacial reaction of the KIS tracer is assumed in the model as a homogeneous reaction taking place at the diffusive interface, where the two fluids mix. The diffusive interface is located at the region with  ∈ (−1,1), and due to the smooth change in  there is no clear boundary.

Solute Distribution Across the Interface
As mentioned in the introduction, the reaction product (2-NSA) has high water solubility and distributes homogeneously into the water phase. This selective distribution mechanism is not included in the classical advection-diffusion equation. In this case, we applied the method developed by (Haroun et al., 2010), where a partitioning coefficient is introduced in the formulation. The partitioning coefficient is expressed as the concentration ratio of the solute in the two-phase fluids system at equilibrium.
where c nw,k and c w,k are the concentration of the product (2-NSA) in the nonwetting phase and the wetting phase respectively, and P ow,k is the partition coefficient.

Continuous Transport Equation
To obtain the single-field formulation of the (2-NSA) concentration for both phases in the entire domain, we define the global concentration as: The sum of the mass conservation Equation 18 for solute k in both phases writes: By substituting Equations 17 and 21 into the above equation, and adding the molecular diffusion term and the interfacial reaction term, one obtains the conservation of global concentration which writes as: where the solute molecular diffusion flux for the mixture is defined as (Haroun et al., 2010): As V f,w + V f,nw = 1, here only V f,w is used in the expression, and the molecular diffusion coefficient is ex- The diffusion term Equation 24 indicates an additional flux term resulting from the solubility law. This solubility flux is in the direction normal to the interface, which governs the distribution of the species between the two phases. With the definition of Equations 20 and 21, one obtains: Thus, the final governing equation for global concentration of solute k is reorganized:

Numerical Implementation
The model is implemented into COMSOL Multiphysics™. COMSOL is an interactive environment for simulating different scientific and engineering problems employing the finite element method for spatial discretization (Akhlaghi Amiri & Hamouda, 2013;Tatomir et al., 2018). The governing system comprises four governing equations and four fundamental variables. The coupled governing equations of the Navier-Stokes equation (Equations 9 and 10) and the Cahn-Hilliard equation (Equation 8) are solved for fluid velocity (u), pressure (p), and the phase variable (). And the governing reactive transport equation (Equation 26) is solved for the solute concentration (c k ). The partial differential equations are solved by the COMSOL linear solver PARDISO. Time stepping is solved with the backward Euler method known for its stability. The initial time step and maximum time steps are controlled to be small enough in order to avoid a singularity. The mesh elements consist of regular triangles with side length h. The mesh size is refined according to the complexity of the geometry and thickness of the interface. The computations were performed on a single CPU with 12 cores at 4.3 GHz, and 32 GB RAM.

Numerical Method and Details
We first verify the reactive transport model by applying the newly implemented continuous transport algorithm (Equation 26) (Section 3.1). Then, the model is tested for a drainage process in a 2D slit (Section 3.2). Finally, the model is applied to study the KIS tracer reactive transport in a realistic 2D porous medium geometry proposed by Keller et al. (1997) (Section 3.3).

Model Verification
For the model verification, we provide first verifications of the two-phase flow simulation with the phasefield method. The two-phase flow model is verified in two different cases: a pressure-difference-driven, co-current two-phase flow through a thin channel (Text S1) and a surface tension driven, deformed bubble retraction process in a quiescent domain (Text S2). The important parameters for the phase-field method, such as characteristic mobility, capillary width and mesh size, are studied. Besides, since the model is to be employed to simulate two-phase flow and reactive transport in porous materials, the phase distribution in the regions close to the grain surface needs to be resolved by the model in detail. Thus, we set up a benchmark simulation of the deposited water films on the grain surface (Text S3). Then, we provide a verification of the reactive transport simulation governed by the above Equation 26. The verification is accomplished by comparing the results of the numerical model with those of an analytical solution for a transient interfacial diffusion process (Section 3.1.1), a simultaneous reaction-diffusion process at the interface (Section 3.1.2), and an advection-diffusion process (Section 3.1.3). In these three one-dimensional model verifications, the mesh element size to capillary width ratio is fixed as h/ε = 1.

Diffusion Across the Interface
When a solute dissolves in a two-phase fluid system with nonequal solubilities in each fluid phase, a concentration jump will be established at the interface at equilibrium, assuming this interface is a sharp one. In this section, we study how this concentration jump is dealt with for a diffusive interface by the model. We consider a simple two-phase system without any reaction, and only the wetting phase contains the chemical component 2-NSA initially. The model is setup in 1D as follows: in a 0.2 mm long domain, the interface is located at x = 0.1 mm. Water is located at x > 0.1 mm and the nonwetting phase is located at x < 0.1 mm, as shown in Figure 2a  concentration of c w,0 = 1 mol/m 3 , and the nonaqueous, or in this case the nonwetting, phase has an initial solute (2-NSA) concentration c nw,0 = 0. The molecular diffusion coefficient is D w = D nw = 1 × 10 −9 m 2 /s for both phases. With this modeling setup, solute diffusion across the interface is allowed to take place. This transient diffusion process can be described by the analytical solution from Equation 27 for water and Equation 28 for the nonwetting phase (Bird, 2002): where erf( ) is the error function. The system with very small partition coefficient is the focus of our study, as the 2-NSA has P ow < 0.01, which means it has very small partition into the nonwetting phase.
Then, we provide a verification of how solute distribution is treated at the diffusive interface. With the very low partitioning coefficient (P ow < 0.01), the solute concentration in water is much larger than that in the nonwetting phase c w,k >> c nw,k . Additionally, the solute concentration in water almost remains at c w,k = 1 mol/m 3 because of the weak diffusion across the interface and the absence of the chemical reaction. In this case, Equation 21 can be simplified to  , , k f w w k c V c . This means that the solute concentration of the mixture (c k ) at the interface is approximately linearly proportional to the volume fraction of the water phase (V f,w ), with a ratio of c w,k . By comparing the simulated concentration curve to the corresponding volume fraction of the fluid, we are able to check if the solute concentration distribution at the diffusive interface is accurately solved by the model.

Interfacial Reaction and Molecular Diffusion
The model is verified for the condition that a zero-order reaction is active. We consider a one-dimensional quiescent channel with two phase fluids. The solute is assumed to react at the interface and diffuse into the water. The channel is 0.3 mm in length, with the interface located at x = 0.1 mm. Water is located at x > 0.1 mm and the nonwetting phase (n-octane) is located at x < 0.1 mm, as shown in Figure 2b. The initial 2-NSA concentration in the domain is c 0 = 0 mol/m 3 . The molecular diffusion coefficient is equal to D = 1 × 10 −9 m 2 /s for both phases. The zero-order reaction rate is dif c R = 1 mol/m -3 s. The capillary width for the model is set at ε = 1 × 10 −7 m, thus the reaction region at the interface has a thickness of 4.16ε according to Section 2.3.2.
The above process can be described by the reaction-diffusion equation: which can be solved analytically. For the analytical solution, the interface between two phases is defined as a no flow boundary for the solute, and the zero-order reaction is active in a narrow region with 0.1 mm < x < L. The analytical solution for Equation 29 is shown below (Carslaw & Jaeger, 1959): GAO ET AL.

10.1029/2020WR028572
10 of 25 , erfc is the complementary error function, and L is the thickness of the reaction region for the analytical solution. To be consistent with the model, the thickness of reaction region is set as L = 4.16ε.

Advection Diffusion
Then, the model is verified under the condition of two-phase flow, where the solute advection is active. The inaccurate modeling of solute advection can cause the artificial (unphysical) mass transfer near the interface (Maes & Soulaine, 2018. The artificial mass transfer is induced when the solute advection is not consistent with fluids motion. The artificial mass transfer is reported to be greater in the systems with larger Péclet numbers (Maes & Soulaine, 2018;Yang et al., 2017). To verify that there is no artificial mass transfer, we consider an advection-diffusion process with the solute concentration in both phases at equilibrium. The model is setup in 1D with a 0.2 mm long domain, where the nonaqueous phase locates at 0 < x < 0.02 mm and water locates at 0.02 < x < 0.2 mm initially, as shown in Figure 2c. Water has initial solute concentration c w,0 = 1 mol/m 3 . The nonaqueous phase has initial solute concentration c nw,0 = 0.5 mol/m 3 . The nonaqueous phase (c nw,0 = 0.5 mol/m 3 ) is injected continuously from the left boundary at in U = 0.05 m/s beginning with time t = 0. There is no reaction ( dif c R =0) and the partitioning coefficient is P ow = 0.5, meaning that the solute concentration in both phases is at equilibrium. In this case, the solute concentration along x can be analytically expressed as: The global Péclet number for the system is calculated as: where L c is the characteristic length, in this case equal to the domain length L c = 0.2 mm. The molecular diffusion coefficient D k is swapped 1 × 10 −6 , 1 × 10 −8 , and 1 × 10 −10 m 2 /s, to apply the verifications of the system at different Péclet numbers. The capillary width for the model is set at ε = 5 × 10 −7 m.

Model Test in a 2D Slit
The model is tested for the relevant processes (reaction, two-phase flow, advection, and diffusion) in a single, two-dimensional slit setup. We focus on the effect of fluid flow (advection) and molecular diffusion on the transport of solute for a spectrum of Péclet numbers. The slit is 0.8 mm long with a slit opening of 0.1 mm, and in 2D the top and bottom boundaries are the solid walls accounting for the pore space surface. The slit (0.02 mm < x < 0.8 mm) is initially filled with water and the inlet part (0 < x < 0.02 mm) is initially filled with n-octane, as the initial interface is necessary for the simulation. n-octane is injected from the left side at a constant velocity of U in = 0.01 m/s at t = 0. The contact angle between the two fluids and the boundary wall is set to 45°. The interfacial tension for n-octane/water is equal to σ = 0.0504 N/m (Tatomir et al., 2018), and dif c R = 1 mol/m −2 s. The initial concentration and boundary conditions are: The characteristic length equals to the half of the slit opening L c = 0.05 mm. As the L c is fixed, the Péclet number of the systems is determined by the ratio between the inflow rate and diffusivity. To obtain the systems with different Péclet numbers, we maintain the inflow rate and swap the molecular diffusion coefficients with D = 1 × 10 −8 , D = 1 × 10 −7 , and D = 1 × 10 −6 m 2 /s in three simulations. The mesh element size and capillary width for the model are h = ε = 2 × 10 −6 m.

Model Application for Flow in Porous Media
We applied the PFM-CST model to study KIS tracer reactive transport in a hydrophilic 2D porous medium. The focus of the study is to understand the effect of the induced immobile zones on transport of the reaction product (2-NSA). We assumed a 2D porous medium from a thin slice of Berea sandstone, first presented by Keller et al. (1997) for a micro-model experiment. The geometry has already been used as the computational domain in Basirat et al. (2017) and Rokhforouz and Akhlaghi Amiri (2019). The domain GAO ET AL.

10.1029/2020WR028572
11 of 25 measures 660 × 320 μm in size. The domain is divided into two parts: the main domain with the porous material and an extensive rectangular zone serving as an inlet for the nonwetting phase on the left side (Figure 3a). The domain is discretized with triangular elements with slide lengths h (Figure 3b). The mesh convergence study is implemented on a sub-section of the geometry, which shows the results convergences at h = ε = 1 × 10 −6 m (Text S4). The main domain is initially saturated with water and the inlet part is initially saturated with n-octane, as the initial interface is necessary for numerical stability (Basirat et al., 2017). At t = 0, n-octane with dissolved KIS tracer is injected from the left boundary with a constant velocity U in . The right-hand side of the domain is the outlet, and the upper and lower sides of the domain are no flow boundaries. The initial concentration of the reacted solute in the domain and the concentration at the inlet are both zero, and the solute is produced from the interfacial reaction during the drainage process, according to the KIS tracer theory. The relevant parameters of the nonwetting fluid are ρ octane = 703 kg/ m 3 and μ octane = 0.54 × 10 3 Pa·s, the contact angle θ w = 45°, and interfacial tension for n-octane/water σ = 0.0504 N/m. With the solution of the reaction-diffusion equation Equations 30 and 31, we know that the zero-order reaction rates can only influence the magnitude of the resulting solute concentration curves, but not its shape. This means that the reaction rate does not affect solute transport, and transport of the solute depends on its diffusivity and the bulk fluid motion. As the study focuses on providing an understanding of the transport and the spatial distribution of the tracer reacted product, which is not affected by the absolute value of the zero-order reaction rate, here we use a unit reaction rate of dif c R = 1 mol·m −3 s −1 . Time is expressed in a dimensionless manner using a characteristic time. The characteristic time is defined as  in / c c t L U , where L c is the characteristic length. The length of the domain is chosen as the characteristic length L c = 0.66 mm. Changes in Péclet number of the flow system can influence the transport process of the KIS tracer and its reaction products. The global Péclet number Pe g is determined as the ratio of the inflow velocity and solute diffusivity, as the characteristic length is known. The inflow velocities are chosen to be large enough so that the nonwetting phase can enter the domain, meanwhile the flow is ensured to be laminar. Thus, a total of six simulations are conducted at three different inflow velocities of U in = 0.01 m/s, U in = 0.05 m/s, U in = 0.1 m/s and the molecular diffusion coefficients of the solute D k = 1 × 10 −7 m 2 /s, D k = 1 × 10 −8 m 2 /s, resulting to global Péclet number from Pe g = 66 to Pe g = 6,600. The Pe g for the previous laboratory application of KIS tracer is covered (Tatomir et al., 2018). The inflow velocity chosen here conforms to the velocity spectrum for simulations done by Basirat et al. (2017) for drainage process with a less viscous nonwetting phase, on the same domain geometry.

Results and Discussion
We first show the result of the model verification in Section 4.1. Then, the results of simulating the 2D slit are shown in Section 4.2. Finally, the results of model application to study the KIS tracer reactive transport in a realistic 2D porous medium geometry are shown in Section 4.3.

Model Verification
Simulation results are shown and compared with analytical solutions, for a transient interfacial diffusion process in Section 4.1.1, for a simultaneous reaction-diffusion process at the interface in Section 4.1.2 and for an advection-diffusion process in Section 4.1.3.

Diffusion Across the Interface
The transient diffusion process depends on the partition coefficient of the solute. Figure 4a shows the comparison between breakthrough calculated by the analytical solution at t = 0.5 s for P ow = 1, P ow = 0.1, P ow = 0.01 and P ow = 0.001 and that of the numerical model. It is found that the model results match those GAO ET AL.  of the analytical solution well. With a lower partition coefficient, solute partitioning into the nonwetting phase is reduced. When P ow < 0.01, the diffusion across the interface has already become very small. Thus, P ow = 0.01 is applied for 2-NSA in this study. When the partition coefficient is P ow < 4 × 10 −4 , the model comes to some unphysical results ( Figure S10). The conditions with extremely small partition coefficient (P ow < 4 × 10 −4 ) cannot be solved with this model, and this is one limitation of the model. Figure 4b shows the resulting concentration profiles near the interface at t = 0.5 s. The solid curve shows the analytical solution considering a sharp interface. The dashed curves show the modeling results at ε = 1 × 10 −6 m and ε = 1 × 10 −7 m. It shows that with a thicker interface, the deviation from the sharp interface solution becomes larger, as expected. Besides, the modeled concentration curves (in dashed lines) are smoothly changing across the diffusive interface, which is different from the analytical solution and its discontinuous concentration profile for the sharp interface. We plot the corresponding volume fraction of water in the second y-axis (with scale on the right side) in Figure 4, and it is found that the concentration curve fits the corresponding volume fraction curves. This proves that the 2-NSA concentration is simulated accurately by the model for a mixture condition at the diffuse interface.
GAO ET AL.   Figure 5 shows the results of the reaction-diffusion process at t = 0.2, 1, and 2 s. The integral of the concentration curves over x (area under the curves) indicates the produced 2-NSA mass in the domain, which is increasing with time. Meanwhile, the produced 2-NSA diffuses from the interface toward the bulk volume of water following the concentration gradient. The modeling results fit very well the analytical solution, demonstrating that the reaction-diffusion process can be successfully simulated by our modeling approach. Figure 5 also shows that the mass transfer across the interface is negligible, which indicates that a partitioning coefficient of 0.01 is small enough to describe that the highly water-soluble product 2-NSA remains in the water phase (wetting phase). Figure 6 shows the results of the transport process at t = 1, 2, and 3 ms. The modeling results fit very well the analytical solution, demonstrating that the advection dominated transport process, with the system Péclet number up to Pe g = 10 5 , is simulated successfully by the model. As the advection of solute is consistent with the fluid motion, no artificial mass transfer is observed near the interfacial region.

Model Test in a 2D Slit
The wettability of the solid wall of the pore space is characterized with the contact angle, which is expressed by Equation 14 in the model. The contact angle on the solid wall is accurately simulated by the two-phase flow model, as shown in Figure 7a for the setting of  w = 45°. Additionally, the characteristic mobility χ, which determines the time scale of Cahn-Hilliard diffusion (Equation 8), needs to be chosen pragmatically (Yue et al., 2006). Straining flows can thin or thicken an interface, which must be resisted by large enough Cahn-Hilliard diffusion (Jacqmin, 1999). The capillary number for the system is calculated as: As the inflow velocity remains constant and the channel flow here is strongly affected by capillary forces (Ca = 5e−4), there should be no changes in interface configuration during the flow process (Prokopev et al., 2019;Soares & Thompson, 2009). It is found when χ is increased by two orders of magnitudes compared to the value of inflow velocity, the thickness and configuration of the interface is ensured to be invariant during the drainage process.
GAO ET AL.   . Comparison between concentration profiles of analytical and numerical model simulations of the transport process for Pe g = 10, Pe g = 10 3 , and Pe g = 10 5 at t = 1 ms, t = 2 ms, and t = 3 ms. Figure 7b shows the concentration profiles for Pe g = 1, Pe g = 10, and Pe g = 100 at t = 0.04 s. The difference in the concentration profiles result from the competition between advection and diffusion. At Pe g = 1, the transfer process is dominated by molecular diffusion, which is verified by the model's results that the concentration contour lines are almost straight and vertical to the inflow direction. For the Péclet number being Pe g = 10, advection starts to have an obvious effect and concentration contour lines begin to blend. For Pe g = 100, the transfer process is dominated by advection, and the solute mass produced at the interface shows a trend to concentrate at the center of the slit. For a stable laminar flow condition in the slit, flow at the center of the slit is faster than the regions near the solid walls, where a no-slip boundary is considered, and all streamlines are parallel to each other. When a two-phase system is considered, near the fluid-fluid interface, the receding phase tends to flowing gathering at the center of the slit to restore the stable flow condition, as shown by streamline plot in Figure 7b. Thus, as the effect of the advection, the water-based solute also shows a trend to concentrate at the center of the slit. A similar observation is found in the simulation results from the study for segmented channel two-phase flow by Yang et al. (2017). The same advection effect exists for all three conditions investigated, with different Péclet numbers. The effect becomes stronger when the diffusion is weak. This simulation provides a verification of the new two-phase reactive transport model.  Model results of solute transport in octane/water displacement process in a capillary tube: (a) volume fraction of the water phase for with contact angles 45° at t = 0.04 s (b) Plot of the concentration profiles at the Péclet numbers of Pe g = 1, Pe g = 10, and Pe g = 100, as well as the corresponding streamline profile.

Application of KIS Tracer Transport in Porous Media
We first show the resulting displacement patterns and solute concentration distributions for the three cases in Section 4.3.1. Then, we determine the threshold local Péclet number to distinguish the mobile and immobile zones of the pore space in Section 4.3.2. Finally, we analyze the quantified interfacial area and the reacted solute mass separately in the mobile and immobile zones, in Section 4.3.3. Figure 8 shows the distribution of the fluid phases and solute concentration before breakthrough for all cases. Since the 2-NSA concentration in the immobile water regions near the inlet where the reaction starts early is much larger than that at the moving front, we plot the concentration on a logarithmic scale to better illustrate the distribution of 2-NSA concentration. The fluid phase distribution is shown in Figure 8(left). At breakthrough, the nonwetting phase is more pore-filling when the inflow rate is larger. This means that a larger inflow rate is favorable for the drainage process. This result is in good agreement with Basirat et al. (2017). Besides, when the nonwetting phase reaches the outlet, the volume of the nonwetting phase in the domain shrinks, especially for Case 3 with the lowest inflow rate ( Figure S6). This is caused by the connected flow-paths of the nonwetting phase between inlet and outlet, formed until breakthrough. The GAO ET AL.  nonwetting phase pressure is suddenly decreased leading to the expansion of the wetting phase. This refers to the capillary end effect, which was also observed in the simulations of Basirat et al. (2017). The spatial distribution of concentration is shown in Figure 8(right). The solute concentration is changing about two orders of magnitude for each case. Trapped water clusters with smaller volumes and larger surrounding interfaces have higher concentration, as expected. Concentration at the front is generally smaller than that in trapped water clusters, mostly because the solute mass at the front is connected to the bulk water, and the solute concentration is diluted by advection and molecular diffusion toward the bulk water. It is observed that the spreading of 2-NSA in these trapped water clusters is much slower during drainage at larger Péclet numbers, as expected. The concentration distributions in Case 3 and Case 6 show the highest values. This is because, with lower inflow rate, the two fluid phases remain longer in contact with each other until the front reaches the outlet (breakthrough) which leads to a higher reaction by-product 2-NSA in the water phase. The mass balance analysis for both fluids and solute is conducted for the verification of the results (see Text S5).

Identification of the Mobile and Immobile Zones
The local Péclet number can be calculated as , where u is the velocity, Lc is the characteristic length equal to the average throat width Lc = 8 μm (Rokhforouz & Akhlaghi Amiri, 2019), and D is diffusivity. The mobile zones are generally defined as zones with Pe >> 1 (Chhabra & Shankar, 2018;Smedt & Wierenga, 1984). The threshold local Péclet numbers to separate the flowing zones and stagnant zones can be chosen by analyzing its probability distribution (Aziz et al., 2018). The drainage process for Case 2 at three different times is shown as an example in Figure 9. During the drainage process, the probability GAO ET AL.

10.1029/2020WR028572
17 of 25 distribution ( Figure 9c) turns from single-mode distribution to bimodal distribution, which is in agreement with Aziz et al. (2018) and Hansan et al. (2020). The first peak from the left-side indicting the formed capillary trapped regions is growing, and the second peak indicting the flowing regions is decaying. The reduction of flowing regions can be observed in Figure 9b. With the constant inflow boundary, the formed capillary trapped water is no longer to be drained. However, the capillary trapped water can still have slight local motion due to the surrounding pressure change during drainage ( Figure S7). This leads to the fact that the magnitude of the local Péclet number for capillary trapped region is smaller than the flowing regions, but larger than the regions of the dead-end pores. The magnitude of the local Péclet number for capillary trapped zones is more obviously displayed near the end of drainage (Figure 9c at t c = 0.45). To separate the majority of the capillary trapped regions from the flowing regions, we chose Pe l = 2.4 as the threshold for Case 2 according to Figure 9c. Furthermore, we do a sensitivity analysis by applying a spectrum of local Péclet numbers near the threshold (Pe l = 0.8 to Pe l = 4.8) to quantify the specific interfacial area and solute mass in the mobile zones. The results converge at Pe l = 2.4 ( Figure S8), which proves that the chosen threshold Pel is appropriate. With the same method, we obtain a threshold Pe l = 4.8 for Case 1 and a threshold Pe l = 0.48 for Case 3 ( Figure S9). Since the flow parameters are the same, the identified mobile and immobile zones are the same for Case 1, Case 2, Case 3, and Case 4, Case 5, Case 6, respectively. However, even the lowest threshold Pe l (Pe l = 0.48 in Case 3) is larger than that reported in Aziz et al. (2018). The reasons could be the differences between the geometries used in the two studies, the parameters of the two systems (e.g., injection velocity, solute diffusivity), and in the study of Aziz et al. (2018) the calculation are done when the flow patterns reach steady state.
With the concentration map of the 2-NSA in the mobile zones (enclosed by isoline in magenta color) during the displacement (Figure 9a), it is clear that the 2-NSA carried by the flowing water is mainly produced and located at the menisci at the moving front in the main flow channels. And the flowing water at the front is not connected to the dispersed water clusters trapped in dead-end pores or throats, which means there is no direct 2-NSA mass exchange between them. This part of 2-NSA mass in the flowing regions is measured for the KIS tracer test. With the model results and threshold local Péclet number, we can quantify the interfacial area and the reacted solute mass in the flowing regions to study the relationship between them.

Interpretation of the Tracer Reacted Mass in Flowing Regions
The specific interfacial area (SIA) is calculated by Equation 32: where ,im m wn A is the interface region within V fw ∈ [0.05, 0.95] for flowing and stagnant regions, b wn = 4.16ε is the thickness of the reaction region obtained from Equation 6, with the diffusive interface having a capillary width of ε = 1e−6 m, and V p is the total volume of the porous media that is studied. Mobile and immobile SIA are plotted in Figure 10 for changes in saturations during drainage. The mobile SIA values range between 200 and 1,200 1/m (Figure 10) while the immobile SIA can reach an order of magnitude higher up to 5,300 1/m ( Figure 10). The mobile SIA has strong oscillations during drainage, which is mainly caused by switching of the mobile and immobile interfaces at the front. When one meniscus at the front meets narrow pore necks, the fluid will become stagnant for several time steps, during which the meniscus that originally counted as the mobile interface will becomes temporarily immobile. Meanwhile, the nearby menisci, which are originally immobile, can also become mobile due to the equilibrium of the local pressure. As the menisci have different sizes, this switching of the mobile and immobile interfaces causes the oscillations of the curve of the accounted mobile interfacial area during the displacement. The overall change in mobile SIA is relatively small during drainage, and its values are similar for all cases. This can be attributed to the fact that the mobile interfacial area results from the menisci at the moving front, which depend on the width and geometry of the flow channels. On the other hand, the immobile SIA increases approximately linearly with saturation during the displacement process, due to the formation of trapped water clusters. Furthermore, the reacted solute mass can be directly quantified by: is the volume of void space for flowing and stagnant regions. The temporal change in total 2-NSA mass in the mobile zones (M m ) and immobile zones (M im ) is shown in Figure 10, where the mobile 2-NSA mass is quantified as the sum of the 2-NSA mass in the mobile zone, and in the part exiting the domain through the outlet (Equation 33). It is found that the 2-NSA mass in the mobile zones is much smaller than that in the immobile zones. The curve oscillations for the quantified solute mass are also caused by the switching of the mobile and immobile zones during the Haines Jumps.
We investigate if the 2-NSA mass transported by the mobile water phase can be related to the mobile interfacial area. Knowing the mobile and immobile interfacial areas, we can calculate the 2-NSA mass (M IFA ) produced by them respectively by Equation 34 (Figure 10). Because the total 2-NSA mass matches that calculated from the total interfacial area ( Figure S5), the additional fraction of immobile 2-NSA mass is derived from the retention of part of the 2-NSA mass produced at the mobile interface. We introduced a mobile mass retention term ξ in Equation 34 to formulate the relationship between mobile interface and 2-NSA in flowing water, as well as immobile interface with 2-NSA mass in immobile zones as Equations 35 and 36: and where M RT is the 2-NSA mass calculated by adding the mobile mass retention term (ξ). The curves are better fitted by adding ξ as shown in Figure 10. The deviation of the curves from M and M RT is found near the end of the drainage process, because a large decrease in the mobile zones occurs shortly before breakthrough. The mobile mass retention terms determined for the six studied cases have values ranging from 0.22 to 0.59. At the same inflow rate, the system with larger molecular diffusion coefficients (smaller Pe g ) has larger ξ.
With the same molecular diffusion coefficient, the system at larger inflow rate (larger Pe g ) has smaller ξ.
A zoom-in region for the retention process is shown in Figure 11. In the zoom-in region at t c = 0.129, the solute mass produced and pushed by three moving menisci (mobile interfaces) at the front are approaching an intersection, where the fluid velocities are slowed down due to the pore necks. At t c = 0.144, one meniscus first passes through the pore neck and the nonwetting fluid fills the pore, and the other two menisci are pushed back during pore filling. The reacted solute is separated into two parts: one part that keeps flowing and the other part that is trapped in stagnant zones, as demonstrated in the sub-figure at t c = 0.159. The retention of the reacted solute in the flowing water (at the front) occurs when flowing water is being capillary trapped. The retention process is also along with bursting and splitting of the moving menisci at the front (see Figure 11 at t c = 0.144), and a part of the moving menisci is arrested and turns into immobile interfaces in the immobile zones, including dead-end pores and capillary trapped zones. The retention process happens frequently with the formation of the trapped water and the immobile interface during the displacement, leading to an approximate constant fraction of solute mass produced by mobile interface to be retained, in the relationships of Equations 35 and 36.
Though the flowing water and trapped water are isolated as the receding phase, the solute mass exchange between mobile and immobile zones happens by means of mobile mass being capillary trapped and turning into immobile mass, during the dynamic process. This part of exchanged mass can be described by the mobile mass retention term. The results imply that when the mobile mass retention term is known, the reacted 2-NSA mass in the flowing water measured in KIS-TT can be used to determine the averaged mobile interfacial area during the drainage process (Equation 3).

Summary and Conclusions
We proposed and demonstrated a novel pore-scale model for simulating the KIS tracer reactive transport process for two-phase flow conditions in a porous medium. In the mathematical model, a continuous species transport formula consistent with PFM was derived. One major advantage of the model is that the fluid-fluid interface hydrolysis reaction of the KIS tracer, which is approximated by a zero-order reaction, can be directly implemented in the diffusive interface approach. We first provide three types of verification of the model. By comparing the solute concentration and water volume fraction across the interface in a simple diffusion process, we showed that the solute concentration at the diffusive interface can be accurately calculated by the model. We found that for low partitioning coefficients (P ow = 0.01), the model adequately describes the selective distribution of 2-NSA into water. The simultaneous process of 2-NSA interfacial reaction and diffusion was verified by comparing the numerical results and the analytical solution of the reaction-diffusion equation. Additionally, the advection of fluid and solute are proven to be consistent with no artificial interface mass transfer in an advection-diffusion process, when solute concentration of the two phases are at equilibrium. Next, the model was tested by simulating the drainage process (i.e., wetting phase displacement by the nonwetting phase) in a 2D silt. After testing the solute transport model for different Péclet numbers, we found that the concentration profile changes are diffusion-dominated at Pe g = 1 and advection-dominated at Pe g = 100.
Finally, the model was applied to a realistic 2D porous medium. Six cases with different global Péclet numbers were simulated. The mobile and immobile zones are separated using a threshold local Péclet number, which is obtained by analyzing the probability distribution of the local Péclet number during the drainage process.
To understand how immobile zones affect the detectable solute mass in the outflow-based measurements, we quantified the interfacial area and the reacted 2-NSA mass separately in the mobile and immobile zones. It was observed that the fraction of 2-NSA mass transported by the flowing water is related to the mobile interface. It was found that the mobile SIA is much lower than the total SIA (ca. 10%) and it varies little during the drainage process. We propose an interpolation relationship between the 2-NSA mass in the mobile zones GAO ET AL.
10.1029/2020WR028572 21 of 25 and the mobile interfacial area. The relationship is formulated by adding a mobile mass retention term. This term relates to the part of 2-NSA produced by mobile interface that becomes residual in the immobile zones.
This study showed that the mobile interface, or the moving interface, plays an important role in the solute transport during the displacement process, especially when the solute is interfacial-reacted and dissolves in the receding phase. The KIS-TT, which is capable of measuring the mobile interface, can be a valuable technique in understanding the coupled flow and reactive transport processes in dynamic two-phase flow systems. Future work is required to study if the mobile mass retention term varies when applied to different pore geometries, or when using different fluid -porous media systems, such as changes in fluid viscosity ratio or medium wettability. Future work is also required to improve our macroscopic model of KIS tracer reactive transport (Tatomir et al., 2015(Tatomir et al., , 2018, for further including the different transport patterns in the mobile and immobile zones observed at the pore-scale.

Data Availability Statement
All data used to support this work are reported in the manuscript and the Appendix in the respective tables and figures.