Modelling of kinetic interface sensitive tracers reactive transport in 2D two-phase flow heterogeneous porous media

Fluid-fluid interfacial area plays an important role for mass- and energy-transfer processes across the interface which is relevant in several hydrogeological and engineering applications, e.g. enhanced oil-gas recovery, CO2 storage in geological formations, unconventional geothermal systems, contaminant removal, etc. Kinetic interface sensitive tracers were designed to determine the size of the interface between two fluids by undergoing hydrolysis at the fluid-fluid interface. This study investigates by means of numerical modelling the influence of heterogeneity on the KIS tracer breakthrough curves in six idealized scenarios (S1-S6). It is an extension of the previous work conducted in “one-dimensional” column experiments by Tatomir et al. (2018) [1]. The changes in interfacial area are created by inclusion of heterogeneities at the Darcy-scale. The results show that KIS tracers can be used in two-dimensional experimental setup and can provide information about the size and dynamic evolution of interfacial area. Therefore, this is a first step for the dimensioning of an experimental flume.


Introduction
Developing real-time monitoring techniques able to capture the movement of fluids in porous media is critically important [2], [3]. In particular, recent interest has focused in quantifying the size of the fluid-fluid interfacial area (FIFA) [1], [4]- [9]. Generally, the FIFA can be determined in controlled laboratory conditions by performing microtomographic experiments on millimeter-or centimeter-size samples [10]- [13]. Another laboratory technique is applying tracer methods [14]- [16]. Usually both the μCT-and tracer experiments are conducted at equilibrium. In contrast to these techniques, KIS tracers were developed by [17] to be applied in dynamic flow in porous media systems. [9] showed for the first time the application of these tracers in a controlled dynamic experiment and proposed a method for the determination of the specific IFA. Both studies use n-octane, a non-aqueous phase liquid (NAPL), which is an analogue for the supercritical CO2.
Numerical models capable to simulate the multiphase flow reactive multi-component transport accounting the FIFA have been developed in [2], [18], [19]. The validation of these numerical models was shown by successfully matching the experimental results of NAPL infiltration in a water saturated column. A column experiment can be regarded as a one-dimensional system therefore, the next step towards the application of KIS tracers at field-scale is the dimensioning and testing in laboratory two-dimensional flumes filled with porous media.
Literature observations about laboratory experiments report that even apparently homogeneous samples give flow patterns which are not consistent with standard models [20], [21]. In this sense, small scale heterogeneities cause the flow to develop channels and fingers Fig. 1.
The aim of the current study is to investigate by means of numerical modelling the non-wetting phase plume development together with the KIS tracer transport in a two-dimensional heterogenous porous media system. At the same time aiming to determine the FIFA from the resulting breakthrough curves.

Mathematical and numerical model
The immiscible displacement of water from the nonwetting phase can be mathematically formulated using the mass balance eq. (1) using the generalized Darcy's law to express the velocity from eq.(2). (2) where denotes the phase (with w, as the wetting phase and n as the non-wetting phase), is the phase saturation, α is the phase density, ϕ is the porosity of the matrix is the phase source or sink term, K is the intrinsic permeability, is the relative permeability, is the phase dynamic viscosity, g is the gravity term, is the phase pressure, denotes the capillary pressure and is the apparent velocity of the fluid as given by the extended multiphase Darcy's law.
The system of two-phase flow partial differential equations (1) is closed with equations (3) and (4): It can be observed the system does not include the interfacial area. [22] proposed introducing a simplified equation for the balance of specific IFA, i.e., equation (5) which is added to the continuity equations : where the interface velocity is expressed: represents the production/destruction rate of specific interfacial area. Even though this term is very important experimental data are still not available for its parametrization. If a porous medium is fully saturated by one phase the IFA is null. Once the displacement process begins the IFA is being created leading to changes in . After certain point, the creation of interface stops, and the destruction begins. At this point = 0. Following this observation, we can formulate: where is the strength of the change of . And further it leads to : Joekar-Niasar et al. (2008) used pore-scale network models to generate computationally − − surfaces. These can be approximated by using polynomial expression [23], [24]: ) as defined in [17] are injected dissolved in the non-wetting phase undergo a hydrolysis reaction at the FIFA to produce one acid ( ) and one alcohol ( ℎ ): The hydrolysis of esters reaction at the water/NAPL interface commonly is represented by a first-order reaction: where → is the mass transfer coefficient or reaction rate coefficient of component from phase to phase , , being the wetting or non-wetting phase, , , is the solubility limit of component in phase , and is the actual concentration of in phase . Following the assumptions of [1], [17] the first order reaction in eq. (11) can be simplified to a pseudo-zeroorder kinetic rate as given in eq. (12).
The KIS tracer and its by-products (acid and alcohol) transport in each immiscible phase of the flow phase flow porous medium system can be written: For the hydrodynamical dispersion tensor we use the formulation of [25]: The system of equations (1) -(14) is implemented in the frame of two numerical simulators. The first is the free open-source academic simulator DuMu x [26]. The model implementation and numerical results for column experiments are shown in [19]. The second model is a commercial software, COMSOL Multiphysics v5.2 [2]. A benchmarking of the two models is performed in [1]. The results presented in this work were computed only with the Comsol implemented model.

Modelling scenarios
We investigate six scenarios of an idealized heterogeneous porous medium consisting of a coarse sand with fine sand inclusions (Fig. 2). The ratio of heterogeneous area (0.4 m x 0.4 m) and surrounding domain is constant in all scenarios except for the first one which is the homogeneous case. The heterogeneities produce deviations of the straight front which induce production and destruction of the FIFA.

Fig. 2.
The six modelling scenarios (S1-S6). S1 is the homogenous case, where the domain is filled only with coarse material. The low permeable inclusions (with higher entry pressure) are highlighted with blue.
The two-dimensional domain size is 1m x 1m. The intrinsic permeability y of the background material is 10e-12 m 2 , while the permeability of the fine material (intrusions) is 100 times lower, i.e., 10e-14 m 2 . The inclusions are represented with a different capillary pressure -saturation -specific IFA relationship than the surrounding porous material.
The entry pressure of the coarse-sand, , = 1700 Pa, and respectively , = 3000 Pa for the fine-sand. The Brooks-Corey pore index coefficient is for both = 3. The coarse-sand entry pressure values are based on the mercury intrusion porosimetry laboratory experiments conducted on glass-beads with the mean dimeter d50=240 μm [9]. Evidently, a finer sand, from the same material will have a higher entry pressure value.   The non-wetting saturation at x0.3 (0.3, 0.95) has the highest variability from the four BTC in particular for S2, S3 and S4. This is because the position of the observation point is directly above the lower permeable intrusion on the direction of displacement. In all cases NAPL arrival time is faster than in the homogeneous case S1. Little influence is seen at observation point x0.5, where Sn remains at 70%.

Results
For the determination of the specific interfacial areas the acid concentration BTCs are being plotted in, Fig. 5. The BTCs show a linear increase with time except at locations x03 and x04, where the fluctuation of the nonwetting saturation leads to a fluctuation in specific FIFA and implicitly of the concentration.
For the interpretation of the BTC results we propose the following relationship: Note that the concentration is assumed to have a linear relationship with time and specific FIFA. Table 1 summarized the observed changes in awn in the six scenarios. The highest relative change is observed for scenarios S4, S3 and S6 which indicates a relationship of the size of heterogeneity relative to the direction of flow.

Conclusions
We presented a mathematical model that accounts for immiscible two-phase flow, specific fluid-fluid interfacial area (FIFA) and reactive transport of KIS tracer and reaction products. This research improves the understanding on the production /destruction of the FIFA, for which there is yet no experimental data available to quantify this term. These numerical results are aimed to contribute to the dimensioning of laboratory setups where KIS tracer experiments will be conducted under controlled conditions. Six scenarios (S1-S6) of heterogeneous porous media were investigated to understand the KIS-tracer breakthrough curves and the interfacial area production. Scenario S4 has the highest rate of change in the acid concentration BTC, with a 9.62% decrease, followed by S3 (-7.94%) and S6 (7.31%).
Results show that the location of the observation points leads to differences in the arrival of the non-wetting phase and tracer concentration. Therefore, for the interpretation of the future 2D laboratory BTCs a multiple point sampling system should be considered with the capability to analyse simultaneous samples. Additional modelling studies are required for understanding the importance of the key parameters at macro-scale with the highest influence on the amount of FIFA and its dynamics.
For the future, the design of the field-tracer experiments, will further require multi-disciplinary work of hydrogeologists, engineers, chemist and modellers, as field conditions involve various scales (e.g., presence of fracture, fissures and faults) and types (e.g., mineral heterogeneity, pH and temperature) of heterogeneities.