Effects of macro-scale heterogeneity on the kinetic interface-sensitive tracer test for measuring the fluid–fluid interfacial area in dynamic two-phase flow in porous media

A novel reactive smart tracer method, termed the kinetic interface-sensitive (KIS) tracer test, has been demonstrated in laboratory column experiments to enable measurement of the specific capillary-associated fluid–fluid interfacial area in dynamic two-phase flow displacement processes in porous media. Development of the tracer method towards effective application in real field conditions requires investigation of the influence of the porous media heterogeneity on the front size and the specific interfacial area, and, consequently, in how far a kinetic interface-sensitive tracer experiment, and the corresponding breakthrough curves, are affected. This study employs a two-dimensional Darcy-scale two-phase flow reactive transport model to investigate numerically the KIS tracer transport in heterogeneous porous media. Simulations were carried out for the primary drainage process in a domain formed of fine and coarse porous media. Various heterogeneity patterns, having different numbers of inclusions and different geometrical distributions, were studied. It is shown that the shape of the breakthrough curves can be used as an indicator for quantifying the displacement front roughness, the specific interfacial area in the domain, and the domain heterogeneity, e.g., the existence of preferential flow pathways inside the porous media. The results indicate that when the displacement front roughness is small, the concentration breakthrough curves exhibit a linear increase. The slopes of the breakthrough curves linearly depend on the fraction of the bulk volume occupied by the low-permeability sand inclusions. The volume-averaged specific interfacial area and the size of the transition zone can be determined from the slopes of the breakthrough curves.


Introduction
Understanding front displacement mechanisms and quantifying the fluid-fluid interfacial area in two-phase flow in porous media is important for many geoscientific disciplines, e.g., hydrogeology, environmental sciences, petroleum engineering, water resources, soil mechanics, etc.For example, in the context of CO 2 storage in geological formations, the CO 2 residual and solubility trapping mechanisms are highly influenced by the magnitude of the interface separating the two fluid phases (i.e., brine and supercritical CO 2 ).Therefore, knowledge of the effective interfacial area, its magnitude, and the shape of the injected CO 2 plume is important to provide information on trapping effectiveness within the formation (Niemi et al. 2017).An active area of research is the development of new efficient monitoring techniques able to track the movement of fluids.Since tracer techniques provide direct information about the hydraulic, transport and reactive processes and parameters of a reservoir, they have been widely used to study the movement of gas and liquids.

Front instability during displacement
A great deal of theoretical and experimental research has been dedicated to understanding the immiscible 1 3 displacement process at the pore scale (Lenormand et al. 1988), but these findings cannot be directly applied to the continuum (Darcy) scale.When the invading fluid is less viscous than the defending fluid, the displacement front becomes unstable and forms "fingers", while a stable front has a front roughness that becomes constant with time.Despite decades of research, the prediction of unstable front growth over time remains poorly understood (e.g., Méheust et al. 2002;Heiß et al. 2011).This phenomenon (of forming highly mobile phase channelling through a less mobile phase) known as viscous fingering, results in tree-like structures spreading throughout the network (Lenormand et al. 1988).Front instabilities can be caused by various factors, such as gravity and viscosity differences, and variations in permeability.Channelling, which refers to the formation of preferential pathways caused by variations in rock structure and pore geometry, i.e., variations in permeability, in capillary pressure and relative permeability, should not be confused with viscous fingering (Bouquet et al. 2020).

Experimental studies on heterogeneity and front morphology
Complex entrapment architectures with residual, ganglia, and pools can be produced by the interaction of geologic heterogeneity and the unstable behaviour during displacement.Although the process of mass transfer across the fluid-fluid interfaces is fairly well understood, the same process occurring in field systems under complex entrapment morphologies continues to be the subject of much research (Page et al. 2007).Microscopic scale investigations employing techniques such as high-resolution X-ray microtomography (XMT) (e.g., Porter et al. 2010) or thin micromodels (e.g., Karadimitriou et al. 2016) can be used to gain understanding on the fundamental fluid-fluid interfacial processes and to quantify to fluid-fluid interfacial area.However, at the pore-scale, where a fluid can either invade a space or not, the heterogeneity of the medium is fundamentally different from that at the Darcy (macro) scale, where all pore space can be invaded (Heiß et al. 2011).At the macro scale, the pore-scale heterogeneities are averaged out in the representative element volume and the exact pore geometries are no longer considered.Much of experimental work was conducted in soil column experiments, which can be regarded as quasi one-dimensional (1D) systems.When dealing with field-scale problems, however, the displacement process comes with a great deal of complexity and even though field conditions are never completely reproducible in the laboratory, especially with respect to porous media heterogeneity, there is still a lot to be learned from such endeavours.In this sense, Darcyscale displacement experiments in intermediate-scale flow cells, e.g., sandboxes, were conducted with the purpose of understanding the key processes (e.g., Brusseau et al. 2002;Fagerlund et al. 2007;Page et al. 2007;Heiß et al. 2011;Rasmusson et al. 2017;Van De Ven et al. 2020).In particular, the drainage of a wetting fluid by a nonwetting nonaqueous phase liquid (NAPL) was studied for different soils and structures, such as homogeneous, layered, or highly heterogeneous soils (e.g., Kueper et al. 1989;Dawe et al. 1992;Illangasekare et al. 1995;Oostrom et al. 1999;Zhang and Smith 2002;Grant et al. 2007;Fagerlund et al. 2007;Heiß et al. 2011).Translucent sandbox experiments have the advantage of real-time tracking of the two fluid phases while allowing the quantification of the fluid saturations.Although, these experimental studies are dealing with many aspects related to miscible and immiscible displacement, such as the investigation of macroscopic quasi-trapping, the finger development, morphology of the displacement front, evolution of the transition zone, the residual trapping, or dissolution of the trapped NAPL, etc., none of them have been aimed at measuring the displacement front by means of a (smart) reactive tracer.

Darcy-scale modelling of immiscible displacement and front morphology
When simulating immiscible displacement processes on large spatial scales, Darcy-scale models can capture the morphology of the transition zone at the displacement front.By averaging the impacts of subgrid Darcy-scale heterogeneity, the Darcy-scale model is upscaled.The growth rate of the transition zone between high and low saturation, which is closely correlated with the morphology of the displacement front, is a critical characteristic that should be replicated by the upscaled model.Instead of being a true interface on the pore scale, the displacement front on the Darcy scale should be viewed as the outermost isoline of saturation of the displacing fluid.Prior research on the Darcy-scale upscaling of immiscible displacement has produced results that can be categorized based on the front's stability (e.g., Fig. 1).Noetinger et al. (2004) observed that, when the mobility ratio is favorable, the transition zone becomes stable throughout time and can be quantified by fluid properties, the correlation length, and the variance of the permeability field.For a neutrally stable displacement, the growth of the transition zone correlates with a dispersive term (e.g., Langlo and Espedal 1994;Neuweiler et al. 2003), and for the unstable case without gravity or capillary forces, it is anticipated that the growth will be linear over time.

The kinetic interface sensitive (KIS) tracer
Kinetic interface sensitive tracers constitute a novel category of reactive tracers developed by Schaffer et al. (2013), aiming at the quantification of the nonwetting/wetting phase fluid interface and tracking the front during the application of geological storage of the supercritical CO 2 .The developed stable KIS tracer is a nonpolar hydrolysable phenolic ester (i.e., phenyl naphthalene-2-sulphonate, 2-NSAPh, CAS 62141-80-4), dissolved in the nonaqueous phase (e.g., NAPL, supercritical CO 2 ) and then injected into the aquifer.Upon contact with the formation fluid (brine), the KIS tracer undergoes an irreversible hydrolysis reaction at the interface.
The KIS tracer hydrolyses into two products (an acid and a phenol) that are highly water soluble and the back partitioning into the nonwetting phase can be considered negligible.Because the reaction kinetics of the tracer is well studied in laboratory batch experiments (Schaffer et al. 2013;Tatomir et al. 2018), the measured reaction product concentration can be used to back-calculate the fluid-fluid interfacial area.The hydrolysis reaction at the interface can be described as: The n-octane as the NAPL (with KIS tracer dissolved) was injected to drain an initially water-saturated column packed with the glass beads.Water samples were collected and measured at defined time intervals at the column outlet to provide the breakthrough curves (BTCs) of fluid volumes, and the 2-NSA concentration in the aqueous phase.The experimental data were interpreted by employing the Darcyscale, two-phase flow, reactive transport model for the KIS tracer (Tatomir et al. 2018).The resulting maximum specific capillary-associated interfacial area was measured at between 500 and 540 m -1 for glass beads with a mean diameter of 240 μm (Tatomir et al. 2018).Furthermore, to understand how the hydrodynamically stagnant zones can affect the KIS tracer transport, Gao et al. (2021a) developed a pore-scale, two-phase, reactive transport model for the KIS tracer.By analysing the size of the fluid-fluid interfacial area and the reacted solute mass in the flowing and stagnant zones, Gao et al. (2021a) proposed that the KIS tracer test determines the mobile capillary-associated interfacial area locating at the moving front.Gao et al. (2021b) further consider the condition of variable surface roughness of the mineral grains.They found that a mobile mass retention term needs only to be considered when the "root mean square roughness" is larger than 3 μm, concluding that the KIS tracer method is also feasible for determination of capillary-associated interfacial area when the grain surface is rough (e.g.natural soils).

Motivation
Despite the previous experimental and numerical studies that improved the soundness of the KIS tracer method, the answer to how the Darcy-scale porous media heterogeneity affects the KIS tracer test and the corresponding BTCs is not yet known.Commonly, all natural aquifers exhibit a certain degree of heterogeneity; therefore, for further development of the KIS tracer technique towards field-scale applications, it is necessary to consider the factor of porous media heterogeneity.The presence of heterogeneities can change the shape of a moving front as illustrated in Fig. 1.Initially, the profile of the displacement front is considered for "smoothed conditions", i.e., having a small front roughness, but due to the porous medium heterogeneities and the fluid-matrix interaction, several fingering scenarios can develop, having different front roughnesses.Applying a KIS tracer to such a front would allow the studying of its temporal and spatial evolution by interpreting tracer BTCs.The concept implies that a larger interfacial area leads to an increase in the reacted mass of tracer in the wetting phase.

Objective
The main objective of the study is to understand how the front morphology, i.e., deformation of the displacement front due to the presence of heterogeneity, affects the concentration BTCs obtained from KIS tracer tests.The study aims to determine if there are correlations between the concentration BTCs and the specific interfacial area, the morphology of the front (the size of the transition zone, front length) and the heterogeneity of the porous media.
A Darcy-scale two-phase reactive transport model is applied to simulate and study the KIS tracer transport during the primary drainage process in a two-dimensional (2D) heterogeneous porous media domain.The domain consists of lenses occupied by two kinds of porous media (i.e., fineand coarse-glass beads), creating in total 15 different heterogeneity patterns.The heterogeneities form regular patterns and are selected arbitrarily in order to produce different fingering and channelling distributions, which can also be later implemented in laboratory experiments.
The paper is organized as follows-section 'Mathematical model' explains the mathematical model of the two-phase flow in porous media coupled with reactive tracer transport while accounting for the reaction across the fluid-fluid interface; section 'Numerical model setup and data processing' introduces the problem and the procedure to analyse the tracer data; the results are discussed in section 'Results and discussion'; and section 'Summary and conclusions' lists the main conclusions.

Two-phase flow in porous media model
The drainage process, defined as the wetting phase being displaced by the nonwetting phase, can be mathematically represented by the governing macro-scale conservation equations for immiscible two-phase flow in porous media: where α denotes the phase (with w denoting the wetting phase and n the nonwetting phase), S α is the phase saturation, ρ α is the phase density, ϕ is the porosity, K is the intrinsic permeability tensor, k rα is the relative permeability of the phase, μ α is the phase viscosity, p α is the phase pressure, q α is the volumetric source or sink term, v α is the apparent velocity of the fluid as given by the extended multiphase Darcy's law, and g is the gravity term.The system of partial differential equations is closed with the following equations: where p c is the capillary pressure.
The functional correlation between saturation and capillary pressure has been derived by many researchers (Leverett 1941;Brooks and Corey 1964;van Genuchten 1980).For this research, the macroscopic capillary pressure-saturation relationship determined from laboratory experiments is expressed by the Brooks-Corey model: where p d is the entry pressure, λ is the pore-size distribution parameter, and S e the effective saturation is defined as: where the S wr is the residual wetting phase saturation.
The relative permeability-saturation relationship is written according to the Burdine theorem (Burdine 1953;Helmig 1997): Further assumptions considered for the two-phase flow in the porous media system discussed herein are: (1) the two fluids are immiscible, incompressible and move under isothermal conditions, (2) the solid matrix is incompressible and nondeformable.

KIS tracer reactive transport model
In an immiscible two-phase flow system such as n-octane and water in a column packed with glass beads, the KIS tracer is present only in the nonwetting phase, while the hydrolysis reaction by-products are present only in the wetting phase, i.e., no back-partitioning is allowed or assumed.Therefore, theoretically, one transport equation is required for solving each component (see Eq. 1), ester/KIS tracer, acid, and alcohol).Experimentally, only the acid concentration is measured in the water collected in the effluent (e.g., Tatomir et al. 2018).Furthermore, the components do not influence each other, nor undergo further secondary reactions, and no additional mass transfer process is active (e.g., n-octane dissolved in the water phase, or water in the n-octane phase).Subsequently, due to the zero-order reaction kinetics and the abundant KIS tracer initially present in the nonwetting phase, the depletion and transport of the tracer molecules are not necessary to be considered.Thus, for the formulation of the KIS tracer reactive transport model in immiscible two-phase flow in porous media, one transport equation is sufficient to solve the transport of the acid hydrolysis by-product.
The reactive transport model of the KIS tracer transport in the immiscible drainage process was demonstrated in Tatomir et al. (2018).The transport of the acid is described by an advection-dispersion-reaction equation (Eq.10), and it is solved decoupled from the two-phase flow equation system (Eqs.2-5).
Here, c κ w is the concentration of the component acid (κ) in the aqueous phase (w), D κ pm,w is the dispersion coefficient, and r κ n→w is the reaction source term.The dispersion coefficient can be calculated by D κ pm,w = v w α L , where α L is the longitudinal dispersivity of the acid.α L was determined using fluorescent tracer experiments to be α L = 10 -3 m for water-saturated flow in a glass beads packed column (Tatomir et al. 2018).Note that dispersivity measured in single-phase flow conditions is not necessarily valid when another fluid-phase is present (Karadimitriou et al. 2016).The hydrolysis reaction rate of the KIS tracer across the fluid-fluid interface is expressed as: where a wn is the specific fluid-fluid interfacial area.
A constant reaction rate implies that the concentrations of KIS tracer reaction products depend on the fluid-fluid interfacial area and the duration the two fluids are in contact with each other.The larger the interfacial area, the larger the mass of the reacted product, and similarly, the longer the time the two fluid phases are in contact, the larger the mass of reacted solute.

Explicit interfacial area model
The specific fluid-fluid interfacial area is incorporated into the two-phase reactive transport model by assuming a functional relationship a wn = a wn (S w , p c ), which can be reduced to a single valued relationship for drainage and imbibition processes (Niessner and Hassanizadeh 2009).Tatomir et al. (2013) have used a polynomial expression capable of avoiding nonphysical capillary and saturation values: where a 0 , a 1 , a 2 , a 3 and p max c are the five parameters determined by fitting the preceding polynomial to experimental data or to models.The parameters a 0 , a 1 , a 2 , a 3 and p max c for the two porous media in this study (given in Table 1) are obtained by fitting the polynomial to the thermodynamically driven, implicit interfacial area model from Grant and Gerhard (2007).The model assumes that the interfacial area is directly influenced by the work done on the system, and thus the interfacial area is directly proportional to the area under the pressure-saturation curve (Grant and Gerhard 2007).The detailed procedure can be found in Tatomir et al. (2018).

Model application to immiscible displacement in a flow cell
Intermediate-scale flow cells are the next step in the development of the KIS tracer technique towards its application in real-field conditions.They allow tracking of the fluid phases and moving front in real-time, as well as the quantification of the saturations.Therefore, the simulation domain is an intermediate-scale flow cell, measuring 0.22 m × 0.2 m.

Boundary and initial conditions
The domain is initially fully saturated with water, and the n-octane with dissolved KIS tracer is injected from the top boundary.Therefore, the top boundary which is the inlet is assigned with Neumann boundary conditions with a constant specific flux of the nonwetting fluid at q n = 9.26 × 10 -6 m/s, and the bottom boundary, which is the outlet, with Dirichlet boundary conditions (at constant wetting phase pressure of p w = 1 × 10 5 Pa and nonwetting phase saturation of S n = 0).The left and right sides of the domain are no-flow boundaries, as shown in Fig. 2a.The initial conditions in the domain are p w = 1 × 10 5 Pa and S n = 0, and there is no solute present, i.e., initial concentration c = 0.The solute is produced in the regions where wetting and nonwetting fluids are both present (Eq.1).

Domain geometry
Two different types of porous media-coarse glass-bead sand and fine glass-bead sand-previously used in the column experiments conducted by Tatomir et al. (2018) and ( 12) 2022), are used to create different heterogeneity patterns in the study domain.The heterogeneity patterns are inspired by the laboratory experiments reported in the literature (i.e., Heiß et al. 2011) and are chosen so that they could be implemented in future experimental work.
The study domain is divided into 77 square lenses (11 horizontal and 7 vertical) some filled with fine sand, while the others with coarse sand.The square lenses have side lengths of 0.02 m.By combining fine-and coarse-sand lenses different types of heterogeneous porous media are created.In total, 15 different heterogeneity patterns are considered, as shown in Fig. 2b.
These heterogeneity patterns cover a wide range of geometries by varying the number of inclusions and their spatial distributions, i.e., (1) low-conductive inclusions in a higher-conductive medium forming regular or random patterns (patterns 3-10), (2) high-conductive inclusions with a low degree of connectivity in a low-conductive domain (patterns 11 and 12), and (3) high-conductive inclusions forming high-conductive straight channels (patterns 13 and 14) or tortuous channels (pattern 15; Fig. 2b).
The entire domain consists of three parts: (1) a rectangular region next to the inlet (y > 0.18 m), (2) a rectangular region which is the main study domain in the centre of the domain (0.04 m < y < 0.18 m), and (3) a rectangular region next to the outlet (y < 0.04 m).The two rectangular regions next to the inlet and outlet boundaries are set to avoid any boundary effects on the resulting breakthrough curves.The properties of the fluids and of the two porous materials are given in Table 1.

Parameters of the heterogeneous media
One parameter to describe the macro-scale heterogeneity of the porous medium composed of the two types of sand is the ratio (R fine ) of the domain volume occupied by the fine sand (V fine ) to the volume of the total study domain (V domain ).
Furthermore, the spatial distribution of the high-conductive inclusions can also affect the hydraulic parameters of the porous media and further affect the two-phase flow and tracer transport processes (Hunt et al. 2014).For the patterns with an identical R fine value, when the highly conductive inclusions become better connected (e.g., forming strong preferential flow pathways), the equivalent permeability of the pattern becomes significantly larger (Colecchio et al. 2021).Therefore, the departure of the equivalent permeability from the geometric mean permeability indicates the magnitude of the connectivity of the highly conductive inclusions.The equivalent permeability (K eq ) of each pattern is obtained by simulating a steady-state single-phase flow (with a constant specific flux q = q n at inlet) using: where L is the longitudinal distance between inlet and outlet, q is the specific flux, and p w,in and p w,out are the averaged wetting phase pressure at inlet and outlet, respectively.The geometric mean permeability of a pattern can be calculated as: where k f and k c are the permeabilities of fine sand and coarse sand, respectively, n f and n c are the numbers of lenses occupied by the fine sand and the coarse sand, and n t is the total number of lenses.Thus, the ratio between the effective permeability and geometric mean permeability can be calculated as: The R K is closer to 1 when the high-conductive inclusions are evenly distributed, and R K becomes larger when the connectivity of the high-conductive inclusions along the longitudinal direction (y axis) is increased.

Problem definition
Two scenarios are investigated: case 1: the so-called standard case with the parameters defined in Table 1; and case 2: with the permeability of the heterogeneous fine-sand inclusions reduced to 1 × 10 -13 m 2 in order to enhance the displacement front deformation (i.e., for the geometrical patterns 3-15), while all the other parameters are kept identical to those in case 1.Thus, for the 15 heterogeneity patterns, a total of 28 simulations (patterns 1 and 2 are homogeneous and are only computed for case 1) are performed for understanding how the (macro-scale) porous media heterogeneity affects the concentration BTCs obtained from the KIS tracer experiments.The kinetic rate coefficient, R c KIS n→w , of the tracer is assumed to be 1.5 × 10 -11 kg m -3 s -1 .To monitor the displacement process and analyse the resulting tracer breakthrough curves, a monitoring line is defined at y = 0.05 (see Fig. 2a).

Data processing
The procedures used to interpret the data and to obtain the concentration BTCs and the averaged specific interfacial areas are further discussed.The concentration BTC is ( 14) obtained by calculating the average solute concentration of acid in the water passing through a measurement line: where M BTC and V w BTC are the solute mass flux and the water volume flux passing through the measurement line.These two fluxes can be obtained by integration: and where L indicates the length of the measurement line.Additionally, the total interfacial area in the domain can be calculated as: where V b is the bulk volume of the porous media in the study domain.Thus, the averaged specific interfacial area in the domain can be calculated by: In this study, the transition zone is defined as the region at the front where strong gradients in saturations occur (Heiß et al. 2011).This zone is located between the contour line of S n = 0.01 and the shock saturation (S shock ).The saturation gradient at the front decreases rapidly when the nonwetting saturation is larger than the shock saturation, which is obtained from the fractional flow function (Eq.22): The shock saturation is the solution of the following equation (Bakharev et al. 2020): Additionally, the length of the saturation isoline at the front at S n = S shock is measured, which is termed "front length" in section 'Results and discussion'.The COMSOL built-in data processing tool is used to measure the temporal variation of saturation isoline.The length of the saturation isoline is an indicator of the extent of the front deformation.A longer isoline means the front undergoes a greater deformation.

Results and discussion
The resulting concentration BTCs corresponding to the different heterogeneous patterns are analysed, and their relations to the fraction of bulk volume occupied by each sand are discussed in the following subsection.Then the effects of the front morphology, i.e., the (strong) deformation of the front on the concentration BTCs, are further analysed and discussed.

Influence of porous media heterogeneities on the concentration BTCs (case 1)
The front length is plotted in Fig. 3a.As expected, for the homogeneous patterns (pattern 1 and pattern 2), the front remains undisturbed, and the front length equals the width of the domain (0.22 m).However, it can be observed that for the heterogeneous patterns 3-12, the front lengths oscillate between values of 0.22 and 0.24 m.Even though the hydraulic parameters of the fine sand inclusions are very different compared to the inclusions filled with coarse sand, the shapes of the fronts are not significantly affected.Therefore, similar to the homogeneous patterns, the fronts advance fairly compactly, and thus the variations in the overall size of the front are small.However, the presence of highly conductive channels (patterns 13-15) has a much larger impact on front lengths.The front lengths keep increasing for patterns 13-15, with values of up to 0.28 m.Pattern 14 which has only two straight highly conductive channels exhibits the largest front deformation.
The next matter to consider is how the KIS tracer reaction product acid concentrations BTCs of all heterogeneous patterns are being impacted (Fig. 3b).First, it is observed that the concentration BTCs from the homogeneous patterns, i.e., pattern 1 and pattern 2, both increase linearly.The concentration BTC of pattern 1 (domain filled only with coarse sand) has the lowest slope, while that of pattern 2 (only fine sand) has the steepest slope.The porous medium with fine sand generates a larger specific interfacial area during the drainage process, and thus the production rate of acid (according to Eqs. 1 and 11) is higher, which was to be expected.The simulated linearly increasing trend of the BTC is consistent with the observations of previous experimental studies undertaken by the authors using columns filled with homogeneous porous media (Tatomir et al. 2018).Furthermore, it is observed that the heterogeneous patterns 3-12 also exhibit linearly increasing concentration BTCs.This is mainly because the shapes of the fronts are not much distorted (with the front length increased up to 9%).However, the shapes of the concentration BTCs are more affected for patterns 13, 14 and 15 where the presence of the highly conductive channels leads to a stronger front deformation (with the front length increased up to 27%).
It can be also observed that the tracer first arrival (breakthrough) time occurs earlier for patterns 13, 14 and 15 compared to patterns 1-12, which have similar breakthrough times.Column experiments using homogeneous porous media have shown that the tracer breakthrough and the nonwetting phase arrival at the outlet occur almost at the same time (Tatomir et al. 2018).In this sense, the discrepancies between the concentration and nonwetting phase arrival times are caused by the strong advection of the acid reaction product in the highly conductive pathways (see Appendix).Additionally, a nonlinear behaviour of the concentration BTCs for the heterogeneity patterns with highly conductive pathways (i.e., patterns 13, 14 and 15) is only observed at the early stage of the breakthrough (before 5,000 s).After the early stage (in this case after 5,000 s) the linearly increasing trend is resumed.Figure 3c shows the variation of the transition zone in the domain (calculated according to Eqs. 22 and 23).It is observed that the transition zone is linearly increasing at the beginning of the drainage process when the NAPL (with dissolved KIS tracer) is entering the study domain.When the entire transition zone is inside the study domain (the isoline S n = S shock has entered the study domain) the slope of the increasing curve starts decreasing.There is no more increment of the transient zone resulting from its entrance into the study domain.For the cases presented here, the nonwetting fluid (n-octane) is less viscous than the wetting fluid (water), and, therefore, the transition zone exhibits an unstable behaviour and continues to grow.At the late stage, as the front begins leaving the study domain, the size of the transition zone starts to decrease.By comparing the size of the transition zones obtained from simulating the different heterogeneity patterns, it is found that the homogeneous coarse sand pattern has the largest transition zone, while the homogeneous fine sand one has the smallest.The size of the transition zone is generally larger when the pattern includes more coarse sand lenses.
Two examples of the simulations in case 1 at the time of breakthrough are presented in Fig. 4. The figure compares the randomly distributed heterogeneity pattern (pattern 4) with pattern 15 which exhibits two tortuous channels.The saturation distribution in the first row, pattern 15, shows an obviously rougher and more deformed front than pattern 4. Additionally, the concentration distributions in the second row reveal that both patterns have strong advection of the reacted tracer in the highly conductive pathways, resulting in relatively higher concentrations in these pathways compared to the surrounding area.The specific interfacial area is shown in the third row.It is observed that the fine-sand inclusions have much higher a wn values than the ones of coarse sand, depending on the a wn -S w curve of each sand.The transition zones are plotted in the fourth row.It can be seen that the changing rate of the saturation is much higher for the region near the front (isoline of S n = 0.01), plotted in darker colour.
By plotting the slopes of the concentration BTCs obtained from all fifteen simulated heterogeneity patterns versus the fraction of bulk volume occupied by the fine sand (R fine ), as shown in Fig. 5a, a clear trend is observed by which the slopes of the concentration BTCs linearly increase with the fraction of the bulk volume occupied by the fine-sand inclusions.This trend can be explained by the fact that fine sand has a larger specific interfacial area (a wn ) than coarse sand at similar saturations (according to the a wn -S w curves), leading to a higher total interfacial area (A wn ) in domains with more fine-sand inclusions.
Additionally, Fig. 5b presents the correlation between the slope of concentration BTC and the specific interfacial area averaged over the whole domain.The results imply that the average specific interfacial area in any heterogeneous pattern can be determined only by measuring the slope of the BTCs and having prior knowledge of the two homogeneous cases.However, it should be noted that these findings are limited to the situations where the tracer BTC has a constant slope, i.e., where the roughness of the displacement front is small.
In contrast, there is no correlation between the maximum front length and the slope of the concentration BTCs as shown in Fig. 5c.The increment of isoline length is up to 27% of the undisturbed length in case 1.Similarly, the size of the transition zone is found to be inversely proportional to the slope of the BTCs, as illustrated in Fig. 5d.This is mainly due to the fact that a larger concentration BTC slope indicates a higher fraction of fine sand (R fine ) in the domain, while in the fine sand, smaller transition zones are formed than in the coarse sand (see Fig. 3c).

Analysing concentration BTCs for the heterogeneity patterns with strong front deformation (case 2)
By increasing the permeability contrast between the background porous media (i.e., coarse sand) and that of the inclusions (i.e., fine sand) the deformation of the front becomes larger.In case 2, the permeability of the finer sand is set to 1 × 10 -13 m 2 , which is 100 times smaller than that of the coarser sand.All other hydraulic and fluid parameters are kept the same (see section 'Model application to immiscible displacement in a flow cell').The ratios of the equivalent permeability to the geometric mean permeability, i.e., R K , for all patterns from both case 1 and case 2 are plotted in Fig. 6a.
As seen from the results, R K from case 2 (<7.59) is generally much larger than that from case 1 (<2.02).In each case, patterns 13, 14 and 15, which feature better connectivity of the inclusions of the coarser sand in the longitudinal direction, show significantly larger R K .The front lengths for all simulations in both cases are plotted in Fig. 6b, and it is evident that the front length is generally larger when R K is larger.This is due to the stretching of the front when highly conductive channels are formed in the simulations with high R K values.The resulting concentration BTCs for all patterns of heterogeneity in case 2 are given in Fig. 6c.Comparing these to the BTCs from case 1 (Fig. 3b), one can observe that the patterns with sparse low-permeability inclusions (5-10) still have a linear increase, and the slopes of the BTCs remain mostly the same.However, the first arrivals (i.e., concentration breakthrough times) are obviously occurring earlier.Patterns 3, 4, 11 and 12 now exhibit small concentration peaks shortly after breakthrough, while in case 1 the BTCs were linearly increasing.For the patterns with conductive channels (13, 14 and 15), where the concentration BTCs show early peaks in both cases, there is a decrease in the slope of the BTCs at later times.This implies that when large changes in the front morphology are involved, the tracer concentration BTCs show some or all the aforementioned features, e.g., an earlier breakthrough time, a local maximum value in concentration shortly after breakthrough, and a lower/flatter slope of the BTC at later times.These nonlinear features of the concentration BTCs for heterogeneous porous media align with the experimental results conducted in sand-filled columns reported in Tatomir et al. (2022), suggesting that these nonlinear BTCs may be caused by heterogeneity created during the packing process of the column, despite using uniform grain size.
The arrival time/ breakthrough time of tracer is plotted versus the isoline length in Fig. 6d.The maximum increment in front length is up to 195% comparing to the homogeneous cases, indicating a strong deformation of the front in case 2. It is shown that when the front length is larger the tracer arrival time occurs earlier.Due to the presence of highly conductive channels in the domain, a more significant deformation of the front occurs, corresponding to an increase of the front length and to an earlier tracer arrival.Thus, the tracer breakthrough time is also related to the morphology of the displacement front.The existence of highly conductive channels cannot be correlated to K eq and there is no correlation between the K eq and the front length and the transition size (Appendix).
One example to illustrate the changing morphology of the invading front when the permeability of the fine sand is reduced is given in Fig. 7a for the heterogeneity pattern 4. The spatial distribution of the nonwetting phase for case 2 heterogeneity pattern 4 (Fig. 7a), with R K = 2.54, is compared with case 1 (Fig. 4), with R K = 1.45.The front shows a more intense folding/deformation when the permeability of the fine sand is reduced, and the isoline length for case 2 (Fig. 6d) is much larger than that from case 1 (Fig. 5c).Due to the lower permeability, the fine sand is not drained as soon as the front passes by, and thus some water-saturated regions in the domain become entrapped, and "holes" appear at the Fig. 6 Plot of a the ratio of the equivalent permeability to the geometric mean permeability, b the maximum front length versus R K , c the concentration BTCs at the measurement line, and d the arrival time versus the isoline length, of case 2 for all simulated heterogeneity patterns with the inclusion lens permeability of 10 -13 m 2 , when the deformation of the front is enhanced front.It can be observed in Fig. 7a that the preferential flow in the coarse sand inclusions increases the acid advective transport in and towards these pathways, leading to a much higher acid concentration in the coarse sand lenses.Figure 7b shows the change of the concentration BTCs for pattern 4. It is found that a nonlinear (deformed) concentration BTC can be divided into three stages.Stage 1 is a nonlinear increase in concentration due to the advection of tracer by-product in the preferential flow pathways.It can be seen that the concentration breakthrough happens at approximately t = 2,200 s, when the front has not yet arrived at the monitoring/control line.At approximately t = 3,000 s, the concentration BTC reaches a peak, which happens at approximately the same time as the arrival of the tip of the front (fingers) and marks the beginning of stage 2.During stage 2, the acid concentration continues to increase but at a smaller rate.This is because the interfacial areas responsible for the acid production, i.e., where the displacement front is present, are mainly the areas located in the preferential flow pathways, where the specific interfacial areas are smaller than in the fine-sand regions where the nonwetting phase is entering much slower.
During stage 2, most of the cross-sectional area at outlet (i.e., monitoring/control line) is not drained and contributes insignificantly to tracer BTC.As more areas become drained the slope of the tracer BTCs increases until reaching stage 3. Stage 3 is the period for which the largest proportion of the displacement front has passed through the measurement line (at approximately t = 5,400 s), and the slope of concentration BTC reaches a stable value that is close to the conditions when the front does not exhibit large deformations.During stage 3, the entire cross-sectional area at the monitoring line contributes to the BTC, which indicates that the interfacial area contributing to the tracer production remains constant.The tracer concentration is averaged over the monitoring line and represents the overall specific interfacial area in the domain.The similar behaviour of the concentration BTCs can also be found in patterns 3 and 11-15, where R K is larger than 2.37 (Fig. 6c).
The results imply that when a strong deviation from a linear behaviour of the front is involved, an accurate determination of the averaged specific interfacial area can only be accomplished based on the slope of the concentration BTC measured for stage 3. Furthermore, with the occurrence of preferential flow pathways, stage 3 is reached at much later times-for example, for patterns 13 and 14, until t = 7,000 s, the concentration BTCs are still at stage 2 (Fig. 6c).

Summary and conclusions
Heterogeneities can have a profound effect on the fluid displacement patterns occurring during immiscible displacement.This study investigated how macro-scale porous media heterogeneity and front morphology affect the KIS tracer breakthrough and its spatial-temporal concentration distribution.It was shown that KIS tracers can provide important information about the front morphology, and it is capable of tracking its dynamic evolution.
Fifteen heterogeneity patterns were created using two types of sands previously used in experimental studies (Tatomir et al. 2018(Tatomir et al. , 2022)).The patterns cover a wide range of heterogeneous conditions from an even distribution to distributions with higher longitudinal connectivity of the inclusions.The simulations were carried out for the primary drainage process (i.e., NAPL with dissolved KIS tracer being injected to displace water initially present in the porous material).The concentration BTCs were obtained by calculating the mass flux and the water volume flux across a monitoring line.Two cases were studied: case 1 with the laboratory-determined parameters of two glass beads, and case 2 where the permeability of fine sand is decreased to enhance the deformation of the front.The major findings of the study are listed in the following: • The subsequent correlations were found between the tracer BTCs and the properties of the porous media properties: (1) the slope of the concentration BTCs linearly depends on the relative fraction of bulk volume occupied by each type of sand; (2) the slope of the concentration BTCs is linearly related to the spatially averaged specific interfacial area in the domain; and (3) the size of the transition zone is linearly related to the slope of the concentration BTCs.• The front deformation (roughness) significantly increases when the permeability ratio between the two sands is increased (i.e., case 2) and when the inclusions allow the formation of preferential flow paths, while the front is only slightly affected when the inclusions are more uniformly distributed.It was observed that when the front length is larger, it has an inverse linear correlation with the tracer breakthrough time; thus, the tracers arrival time is shorter.• When the front morphology is not much affected by the heterogeneous inclusions (e.g. the increment of the isoline length up to 9% of the undisturbed length), the concentration BTCs mainly appear as linearly increasing curves.• When the displacement front exhibits strong deformation, the concentration BTCs, which otherwise display a linearly increasing trend, exhibit deviations from the linear profile.The nonlinear concentration BTCs show three main features (Fig. 7): an earlier first arrival time, a peak of concentration shortly after first arrival, and a lower slope during the following stage stabilizing eventually at a constant slope.Based on the slopes of the concentration BTCs, the KIS tracer test in heterogeneous porous media can mainly be divided into three stages-stage 1: nonlinear increase of tracer concentration until the arrival of the first fingers of the nonwetting front, which coincides with a concentration peak; stage 2: (after the peak) gradual increase in concentration reaching eventually a stable value for the slope, which represents the beginning of stage 3; stage 3: the slope of the BTC continues to increase linearly at a constant slope close to the conditions when the front roughness is small (i.e., it is not folded).
The results imply that the KIS tracer test is capable of measuring the averaged specific interfacial area in heterogeneous porous media.The measurement time needed is much longer for the concentration BTC to reach the stabilized slope when short-cuts and preferential flow pathways are involved.It is shown that the KIS tracer test can provide information not only on the averaged specific capillaryassociated interfacial area, but also on the volume fraction in the bulk volume that is occupied by the inclusions.Besides, the results also imply that, based on the shape of the obtained concentration BTCs, one can predict the extent of front deformation, as well as the extent of heterogeneity present in the studied porous media.
By considering porous media heterogeneity, this work represents a milestone for further development of the KIS tracer technique towards being deployed in field-scale applications, e.g., the tracer could be applied in singlewell push-pull tests or in well doublet configurations.Further research is required to develop the experimental procedure of applying the KIS tracer test in intermediatescale flow cells.The flow cells allow the real-time tracking of the fluid phases, the quantification of the fluid saturations and the measurement of concentrations.This is important for improving the understanding of reactive transport in two-phase flow in porous media processes, and for providing a further validation of the numerical model in complex flow fields induced by the heterogeneous porous media.
reaction rate coefficient, where n and w denote the nonwetting and wetting phase.Static batch experiments showing the hydrolysis of 2-NSAPh in n-octane/pure water fluid systems resulted in a linear relationship between concentration curves of 2-NSA over time, indicating a zeroorder kinetic rate.This represents a key assumption towards the suitability of 2-NSAPh as KIS tracer (Schaffer et al. 2013).Tatomir et al. (2016) established a mathematical and numerical framework for the design and application of the KIS tracer in laboratory experiments and theoretical studies.The KIS tracer is designed to measure the fluid-fluid interfacial area in a transient immiscible drainage process.The concept of tracer application was first validated by Tatomir et al. (2018) in controlled column experiments with a wellcharacterized porous medium composed of glass beads.

Fig. 1
Fig. 1 Schematic illustration of the development of a fingering front in a linear displacement starting at initial conditions (left): smoothed, conserved, increased, branched front shapes (classification after Bakharev et al. 2020).On the right, the zoom-in pore-scale region showing a fluid-fluid interface where the hydrolysis reaction of the kinetic interface sensitive tracer happens (the KIS tracer hydrolysis reaction, as in Eq. (1), is represented in the figure by A: phenyl naphthalene-2-sulphonate, B: Naphthalene-2-sulphonate acid, and C: Phenol) The set of partial differential equations are written into the PDE interface from the basic module and are solved by the multifrontal massively parallel sparse direct solver (MUMPS).The governing equations can use the pressure-saturation formulation(Helmig 1997) or the global pressure formulation(Douarche et al. 2022), with the former being applied in this study because it facilitates the sequential solution of the equations.Tatomir et al. (2016) implemented the COMSOL model and conducted a sensitivity analysis to determine the KIS tracer behaviour with regard to changes in flow and transport parameters.Tatomir et al. (2018) provided a validation of the COMSOL implemented model by comparing it to the results obtained from a drainage experiment in a column packed with glass beads.The COMSOL model was further verified by the same group with a code inter-comparison to DuMu X(Flemisch et al. 2011), an open-source academic code for simulation of flow and transport processes in porous media(Tatomir et al. 2018).The spatial discretization is based on the finite element method (FEM), and the time discretization is backward Euler.The initial time step is 0.001 s and maximum time steps is 5 s.The mesh elements are regular triangles with side length of 2.6 mm.In total the mesh consists of 19,662 elements.All computations were performed on a single CPU with eight cores, operating at 4.3 GHz, and 128 GB RAM.

Fig. 2
Fig. 2 Plot of a the simulation domain (with the boundary conditions, initial conditions and monitoring line labelled) and b the heterogeneity patterns used in the modelling of the KIS tracer experiments

Fig. 3
Fig. 3 Plot of the a length of displacement front (isoline at S n = S shock ); b concentration breakthrough curves at the monitoring line; and c size of the transition zone in the study domain.There is variation in time during the drainage process in case 1 for all simulated heterogeneity patterns

Fig. 4
Fig. 4 Plot of the distribution of the nonwetting phase saturation, acid concentration in water, specific interfacial area, and the changing rate of nonwetting phase saturation (in the transition zone between S n = 0.01 and S shock ) at the time of concentration breakthrough of case 1, for pattern 4 (t = 3,300 s) and pattern 15 (t = 3,100 s)

Fig. 5
Fig. 5 Plot of a the slope of the concentration BTCs as a function of the fraction of total bulk volume occupied by the fine sand, b the averaged specific interfacial area at the end of the drainage process versus the slopes of the concentration BTCs, c the maximum isoline length during drainage versus the slopes of the concentration BTCs, and d the maximum size of the transition zone during drainage versus the slopes of the concentration BTCs, for all heterogeneous patterns (for case 1)

Fig. 7
Fig. 7Plot of a the nonwetting phase saturation and concentration distribution at three different time steps at t = 2,200, 3,000 and 5,400 s, from simulations of pattern 4 for case 2 (the grey arrows are the streamlines computed for the wetting-phase velocity, showing the advective transport towards and into the conductive lenses), and b comparison of concentration breakthrough at the monitoring line for simulation of pattern 4 before (case 1) and after (case 2) the modification of the permeability of the inclusions

Table 1
Parameters