Abstract

This pore-scale modeling study in single self-affine fractures showed that the heterogeneous flow field had a significant influence on the mixing-induced reaction transport. We generated the single self-affine fracture by the successive random additions (SRA) technique. The pore-scale model was developed by coupling the Navier-Stoke equation (NSE) and advection-diffusion equation with reaction (ADER). Eddies were captured in the self-affine fracture due to the increasing Reynolds number and the sudden expansion of aperture. The flux-weighted breakthrough curves (BTCs) of reaction product showed the typical non-Fickian characteristics (i.e., “early arrival” and “heavy tail”). It was found that the reactant was involved in eddies and then reacted inside the eddy-controlled domain. Consequently, eddies played a significant role in delaying the mass exchange process between the eddy-controlled domain and the main flow channel, which resulted in the “heavy tail” in BTCs. As the Reynolds number increased, the breakthrough time increased while the concentration peaks of BTCs decreased. Furthermore, the dilution index presenting the exponential of the Shannon entropy of a concentration probability distribution was used to quantify the degree of reactant mixing. The results showed that the quantification of dilution for nonreaction transport was in good agreement with the outcomes of mixing-induced reaction transport. The high Reynolds number and Peclet number had a negative influence on the mixing process at the early time whereas they led to the enhanced mixing process at the late time.

1. Introduction

Reactive transport dominates the fate of groundwater contaminants in fractured rocks. Characterizing this fundamental transport process is of key importance to not only groundwater remediation but also many natural and engineered systems in the subsurface such as enhanced oil recovery, geothermal energy development, and electrical conductivity of geological formations [1, 2]. As a classical macroscopic-scale approach, the advection-dispersion reaction equation (ADRE), based on a perfect mixing assumption, commonly overestimates the reaction rate and tends to generate a less reliable prediction. Many studies [39] have shown that the ADRE is generally incapable of describing the local mixing phenomenon and the neglected pore-scale incomplete mixing is a primary reason for the overestimated reaction rate. Therefore, the mixing-induced reaction transport through the fracture still needs to be better understood at the pore scale.

Since transport behavior of solute is intrinsically coupled with the fluid flow, the literatures on the fluid flow in the rough fractures are briefly reviewed for the sake of completeness. The roughness of fracture wall effecting the fluid flow has attracted much attention in the last serval decades [10, 11]. Based on the ideal smooth parallel plate model, the classical cubic law, which is the analytical solution of Reynolds equation derived from the linearization of the NSE, is widely used to simplify the fluid flow in single rough fractures. However, many numerical and experimental studies reported [12, 13] that the classical cubic law was valid only when the flow was laminar and may result in nonnegligible errors when the fracture wall was rough. The assumption of the ideal smooth parallel plate model at the macroscopic scale is the primary factor which leads to the fact that the cubic law may be insufficient for describing fluid flow and reaction transport. In addition, the geological fracture walls can be described as self-affine structures for length scales ranging from microns to meters [14]. The self-affine structures can be generated by the successive random additions (SRA) technique [15, 16]. Therefore, in this study, the flow field in the single self-affine fractures was solved directly by NSE at the pore scale.

The traditional macroscopic reaction transport models rely on upscaling the transport equation and effective transport parameters. The common procedure is to determine effective transport parameters by fitting model results to flux-averaged breakthrough curves (BTCs) of a passive tracer. Due to the heterogeneity of geological formation, the BTCs might be non-Gaussian-shaped and could exhibit the non-Fickian characteristics (i.e., “early arrival” and “heavy tail”). Although applying macroscopic advection-dispersion equation (ADE) to fit the non-Fickian BTCs could yield a nonnegligible error, the non-Fickian BTCs may provide important information for understanding the mixing process. For example, “early arrival” in BTCs implies that the preferential path (channelling) may occur in fractures while “heavy tail” indicates that some slow mass exchange process controls transport at the late time [17, 18]. Thus, one can expect that the non-Fickian BTCs with “heavy tail” is capable of serving as an indicator of mixing process occurring in the heterogeneous geological formation. Luo and Cirpka [19] numerically examined the validity of macroscopic transport models, which are capable of describing all details of BTCs, for predicting a reaction transport in heterogeneous media. They found that directly applying the macroscopic transport model to reaction transport is inappropriate due to the presence of the concentration variations. Furthermore, the mixing process, which differs from the spreading process, represents the distribution of a solute over an increasingly larger volume. This is the only process that allows mass exchange between different streamlines and results in the decrease of its peak concentration [20]. In fact, the mixing can also be considered as the process that smears out concentration contrasts. The reactant mixing process is crucial for accurately predicting the mixing-induced reaction transport. However, the influence of eddies on the mixing process in a self-affine fracture is still an open topic.

The main objective of this work was to investigate mixing-induced reaction transport through single self-affine fractures by using a pore-scale model. The self-affine fracture was generated by the SRA technique. The flow field in the single self-affine fracture was numerically solved by NSE. The nonreaction transport and the mixing-induced reaction transport were simulated by the ADRE. The flow field and eddies formation were first examined under four different Reynolds numbers (Re). Sequentially, the non-Fickian BTCs for nonreaction transport and the mixing-induced reaction transport were analyzed. Finally, reaction products under different Re and Peclet number (Pe) were calculated. The corresponding reactant mixing process, which was quantified by the dilution index, was investigated.

2. Modeling Approach

2.1. Single Self-Affine Fracture Generation

Following the seminal work of Mandelbrot and Pignoni [21], the height of the self-affine rough fracture wall can be described as where is the constant, called scaling factor, and the exponent is a measure of the fracture roughness and defined as the Hurst exponent with the range from 0 to 1. For generating self-affine fractures, the specific Hurst exponent must be assumed. A large number of experimental studies reveal that the Hurst exponent of the natural rock fracture wall varies between 0.47 and 0.85. The Hurst exponent was proved to depend on the mineral component of rocks. The characteristic Hurst exponent of the fracture wall was around for granite and basalt, whereas it was for sandstone [22]. However, our primary objective of this study is not to perform an exhaustive analysis for the different Hurst exponent of the fracture walls. The Hurst exponent suggested to be universal by the previous works [23, 24] was assumed in the current study.

Once the Hurst exponent of the fracture wall is specific, several methods are available for generating the self-affine fracture wall (e.g., the SRA and the Fourier transformation). In this study, we used the SRA technique [15, 16] to generate the self-affine fracture wall. The application of the SRA can be found in our previous works [23, 24]. From the generated self-affine fracture wall to the single self-affine fracture, we used the shear displacement model [25] to construct the single self-affine fracture, associating with the local aperture as a function of longitudinal distance (see Figure 1). In order to improve the possibility of the occurrence of eddies, we used a relatively large standard deviation  mm and mean aperture  mm. This results in the coefficient of variation, , equal to 0.52. The total length of the generated fracture was equal to L=0.16 m.

2.2. Flow Model

The flow field in a single self-affine fracture is solved directly by using the continuity equation and NSE, which is assumed as the isothermal, incompressible, and homogenous single Newtonian steady flow:where is the density of fluid, is the velocity vector, is the total pressure, and is the dynamic viscosity of fluid. In this study, we used standard water properties at 20°C, e.g., =998.2 kg/m3 and =1.002×10−3 Pas. The self-affine fracture walls were considered as the nonslip boundaries. The steady-state flow was induced from left to right by a given pressure drop over the entire fracture. The flow model was solved based on the Finite Element Method (FEM) and implemented in COMSOL Multiphysics. There are about 148,000 triangular elements in the discretized fracture domain where the maximum and minimum element sizes were imposed as and , respectively. The results of the solved flow field served as the basis and were coupled with the solute transport model. The results of the mesh independence showed the steady-state flow rate under the pressure gradient =100 Pa/m varies from to when the amount of triangular elements was changed from 148,000 to 365,000. This indicates that the 148,000 triangular elements are sufficient to obtain the stable and accurate simulation results.

2.3. Solute Transport Model: Nonreaction and Reaction

The solute transport model based on the ADRE was used here to simulate the nonreaction and reaction transport in the single self-affine fracture. Based on the classical Fick’s law, the general equation of the ADRE was given aswhere is the solute concentration at the space position at time , is the flow velocity tensor, the is the molecular diffusion for the solute species , and is the reaction rate of the solute species . As a special case of ADRE, (4) could be recovered to the classical advection-diffusion equation when . Thus, to simulate the nonreactive solute transport, we assume that the reaction rate is equal to zero and the nonreactive tracer is instantaneously introduced in a tiny rectangular region (0.5 mm×2 mm). There are no absorption and reaction between the fracture wall and the tracer. The fracture matrix is assumed to be impermeable. The boundary conditions of the inlet and outlet for the nonreactive solute transport are set as zero concentration and zero longitudinal concentration gradients, respectively. The flux-weight breakthrough curves (BTCs) can be obtained by a ratio of effluent solute mass to fluid mass:By normalizing flux-weight BTCs and time, where is the dimensionless concentration, is the pore volume (PV), is the flow rate per length, and is the area of the fracture. Considering average flow velocity in the fracture, an important dimensionless number, Peclet number, is given by the ratio of the characteristic diffusion time scale to the characteristic advection time scale :

The reactive solute transport was also performed in this study with consideration of the simple irreversible reaction . Depending on the mass balance law, the reactants and reaction products can be described asWe assume that the compounds A and B cannot coexist and react in a one-to-one molar ratio at the any node, which means that the rate of production of C is equal to the rate of loss of each reactant. Thus the reaction is a second-order irreversible reaction. The general formula for such reaction rate is given by [26]where k is the reaction rate constant. The Damkohler number (Da) is defined as the ratio of advection (residence) time scale to reaction time scale . For the chemical kinetics in well-mixing conditions, the reaction time scale and Da are defined by [27]In this study we assume that the reaction between A and B is instantaneous, associating with a relatively large Damkohler number (). For the reactive solute transport, the reactant A, as the invading reactant, is injected instantaneously as the same rectangular pulse near the inlet boundary of the self-affine fracture, while the reactant B is introduced as the surrounding ambient water and is injected continuously with a constant concentration at the inlet boundary where the Dirichlet boundary condition was implemented. To avoid the complete consumption of the invading reactant A, the concentration of the invading reactant A is set to be 500 times higher than the reactant B, by which we can obtain the BTCs for reactant A at the outlet. Under such condition, the reaction between A and B is upon the magnitude of the dilution (mixing) of A. This implies that the reaction transport is mixing-controlled.

3. Problem Setup

This study first considered the cases with four different pressure gradients to simulate the flow field in the generated single self-affine fracture, hereafter, denoted as flow field 1–flow field 4 (FF#1-FF#4). For each FF#1-FF#4, both of the nonreaction and reaction transport simulations with the initial injection mode of the instantaneous rectangular pulse were performed in the same single self-affine fracture. Since the mean aperture and the length of fracture were specific, the advection time scale and diffusion time scale only depended on the average velocity and diffusion, respectively. Due to the sudden expansion of the aperture, the eddies usually form in the vicinity of the rough fracture wall [28]. Thus, to capture the influence of eddies on the spreading and mixing of the solute, the relatively large diffusion time scale was assumed for both the nonreaction and reaction transport, which resulted from a relatively large diffusion .

For the reaction transport, the reaction is assumed to be the irreversible instantaneous reaction . In the flow-reaction system, the properties of the instantaneous reaction are determined by the Da. For the large Damkohler number (), the advection time scale is much larger than the reaction time scale, which indicates that once the reactants encounter each other at any node, the reaction is instantaneous with respect to the reactants transport. In other words, only if reactants are locally mixed, is the relatively low concentration of the reactant completely depleted, and the reaction stops. Thus, the further reaction depended upon the further reactant mixing. Such reaction is that of mixing-induced reaction. For all the reaction transports in our study, the Da ranging from to is much larger than 1 and the reaction is considered to be instantaneously irreversible. Table 1 shows the dimensionless parameters for nonreaction and reaction transport in the single self-affine fracture.

4. Results and Discussion

4.1. Flow Field and Eddies Formation

Understanding flow field in fractures is fundamental and necessary to study the nonreaction and reaction transport. For the sake of completeness, we briefly introduced and discussed the flow fields. Since the occurrence of eddies was highly sensitive to the characteristics of the flow field, we used the four different pressure gradients to drive the flow field in the same single self-affine fracture. Consequently, the four different structures of the flow fields were obtained as shown in Figure 2. Here, we introduced the Reynolds number (Re) to qualify the characteristics of the flow field. The Reynolds number (Re) is defined as a the ratio of inertial forces to viscous forces: where is the flow rate and W is the width of fracture in the out-of-plane direction (equal to 1 m in the 2D problem).

Figure 2 shows a segment (22 mm in length) of a 160 mm fracture with a Hurst exponent of 0.8. At low Reynolds number (Re=30), there is no eddy in the segment. In fact, the flow field under Re=30 shows that the main flow channel occupies all of the cross sections and there is no noticeable eddy-controlled domain at the full-length scale of the fracture. As the Re increases, it can be seen that eddies tend to occupy significant portions of the fracture where the aperture shows the sudden expansion. The main flow channel becomes narrow due to the expansion of the eddy domain. Several previous works [2832] showed similar results for the growth of eddies due to the increasing Reynolds number. At high Reynolds number (Re=170), it can be seen clearly that the eddy-controlled domain occupies a large portion of the fracture.

Due to the occurrence of eddies, the streamlines can be detached from the bulk flow and then are involved in the downstream location again [33]. As a result, the regions of detached flow do not contribute to bulk flow. The previous study on porous media showed that the stagnant areas created by eddies gave rise to enhance the dispersion coefficient [34]. The velocity in the stagnant area is very low and the streamlines meander even change directions. This may cause the delayed solute transport since the local Pe is small and molecular diffusion becomes the dominating mechanism [35].

4.2. Breakthrough Curves

The flux-weighted BTCs of compound A for the nonreaction transport were calculated under the different flow fields in the single self-affine fracture, as shown in Figure 3. For the nonreaction transport, the salient tails can be observed in all BTCs, which indicates that the nonreaction transport is non-Fickian or anomalous. It can be seen in Figure 3 that as Pe and Re increase, the tails in BTCs become heavy. Since there is no visible eddy in the flow field with the small Pe and Re, the variable aperture distribution in the fracture is primarily responsible for the non-Fickian characteristic (i.e., tail). When Re is large, eddies play a significant role in enhancing the tails in BTCs. Figure 4 clearly shows the retention of compound A in the vicinity of the rough fracture walls where eddies form and grow. Once the compound is captured by the eddy-controlled domain, the mass exchange between the main flow channel and eddy-controlled domain determines the length of the tail in BTCs for a nonreaction transport.

Although eddies form and develop as the Re increases, there is no tail in the BTCs of compound A for the reaction transport, as shown in Figure 5. This is because the trapped compound A in eddies is eventually consumed by the reaction. This result further demonstrates that the transport of the nonreactive compound is significantly delayed by eddies. Comparing Figures 4 and 5, it can be seen that the mixing-induced reaction leads to the decreasing concentration peak in BTCs for the same Re and Pe. However, this mixing-induced reaction has little influence on the variation of concentration peaks in BTCs with the increasing Re and Pe.

As expected, the BTCs of reaction product C exhibit the heavy tails at the late time. It can be seen in Figure 6 that the length of tails in BTCs is longer in FF#4 than in FF#1, indicating that eddies have influence not only on the nonreaction compound transport but also on the reaction product transport. Since eddies could be considered as a segregated domain for the compound transport, there are two possibilities causing the non-Fickian tail in the BTCs of the reaction product. First, the reaction product was transformed at the outside of the eddy-controlled domain and then captured by eddies. Second, the reaction product was transformed inside the eddy-controlled domain, which required that the reactant was involved in the eddy-controlled domain. Comparing Figures 4 and 7, it was found that the instantaneous transformation was inside the eddy-controlled domain due to the involved reactant. This is important because that the effective chemical kinetics of reaction transport inside the eddy-controlled domain might be significantly different from that in the main flow channel.

Analogous to the nonreaction transport, the length of tails in BTCs is highly dependent on the mass exchange rate between the main flow channel and eddy-controlled domain. It is also found that the breakthrough time of the reaction product decreases as Pe and Re increase. This is mainly due to the narrow main flow channel caused by the growing eddies. Another finding is that the peak concentration in BTCs decreases as Pe and Re increase. This implies that the effective chemical kinetics of reaction transport is influenced by eddies. Since the reaction used in this study is instantaneous, irreversible, and dependent upon the reactant mixing, this influence of eddies on chemical kinetics can be directly reflected by the magnitude of reactant mixing. It should be noted that the analysis of BTCs is only sometimes informative of incomplete mixing process [36], especially in the heterogeneous geological formations.

4.3. Reactants Mixing and Reaction Product

The influence of incomplete mixing on mixing-induced reactions has recently won much attention and been reported in many literatures [3641]. The key point for this is that a fundamental difference exists between spreading and mixing processes. Especially in the heterogeneous flow field, the latter occurs at the small spatial scale and cannot be quantified completely in terms of spatial moments of concentration distribution which represents the spatial extent of compound (i.e., spreading). Since the reaction transport may be influenced by incomplete mixing, there are various methods for quantifying mixing such as scalar dissipation rate [42] and dilution index. Here, the dilution index, introduced by Kitanidis [43] based on entropy concepts, was used to quantify the global compound dilution (mixing) in the self-affine rough fracture and given bywhere is the concentration probability function of the compound A defined by The maximum value of the original dilution index, , is technically equal to the volume of the entire fracture, representing a nonuniform concentration distribution in fracture. Then the dimensionless dilution index can be obtained by .

In order to investigate the mixing process under the different flow fields, the direct simulations of nonreaction transport for compound A were conducted with the same initial concentration and boundary conditions. The dimensionless dilution index was calculated in each flow field before the compound arrived at the outlet boundary, as shown in Figure 8. It can be seen that the dilution index increases with pore volume and the flow field has a significant influence on the temporal evolution of the dilution index. At the early time (PV<0.22), the dilution index is higher in FF#1 than in FF#4. This implies that high Re and Pe have a negative influence on the mixing process at the early time. However, at the late time (PV>0.22), the dilution index is lower in FF#1 than in FF#4. This indicates that high Re and Pe lead to the enhanced mixing process at the late time. Increasing velocity causes the increasing stretching and deformation of the plume, which in turn enhances the diffusive mass exchange between the plume and ambient water. Thus, the mixing process becomes complete as Re and Pe increase.

Since the BTCs cannot accurately reflect the chemical kinetics for the mixing-induced reaction transport, a global measure of the chemical kinetics of the reaction transport is defined in terms of the temporal evolution of the reaction product mass within compound A through the single self-affine fracture. The mass of the reaction product C can be obtained asThe temporal evolution of the mass of the reaction product C is given in Figure 9 for FF#1-FF#4. In each case, the mass of the reaction product C is normalized by the maximum mass , where is the breakthrough time of the reaction product C. It can be seen in Figure 9 that the amount of the reaction product is lower in FF#4 than in FF#1 for PV<0.3, while it is higher in FF#4 than in FF#1 for PV>0.3. In general, these results are in good agreement with the temporal evolution of dilution index of compound A for the mixing-induced nonreaction transport as shown in Figure 8. This suggests that the chemical kinetics of the mixing-induced reaction transport can be qualitatively analyzed by the temporal evolution of the reactants mixing process.

5. Summary and Conclusion

To model the mixing-induced reaction transport through a single self-affine fracture, we generated the single self-affine fracture by the SRA technique and a pore-scale model was developed by coupling the NSE and the ADRE. This allowed us to directly capture fundamental physical processes of flow, nonreaction, and mixing-induced reaction without consideration of parameterizations at larger macroscopic scales. The self-affine fracture wall led to the tortuous streamlines in the flow field. Eddies were observed in the vicinity of fracture walls where the aperture showed the sudden expansion for the high Reynolds number. The BTCs of nonreaction and mixing-induced reaction showed the non-Fickian characteristics such as “early arrival” and “heavy tail”. Eddies were primarily responsible for the “heavy tail” in BTCs since the mass exchange process between the eddy-controlled domain and the main flow channel was significantly delayed by eddies. However, it should be noted that the reactant was first involved in eddies and then reacted inside the eddy-controlled domain. This may result in the different effective chemical kinetics of reaction transport between inside eddy-controlled domain and the main flow channel.

It was found that as Re and Pe increased, the breakthrough time increased while the concentration peaks of BTCs decreased. The degree of reactant mixing was quantified by the dilution index. The quantification of dilution for nonreaction transport was in good agreement with the outcomes of mixing-induced reaction transport. The temporal evolution of reactant mixing indicated that high Re and Pe had a negative influence on the mixing process at the early time whereas they led to the enhanced mixing process at the late time. The pore-scale analysis performed in this study highlights the critical role of the heterogeneous flow field on the mixing-induced reaction transport through the self-affine fracture. However, this role needed to be further derived quantitatively and extended in 3D single fractures and network fractures with consideration of the permeable fracture matrix.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

The study is financially supported by the National Natural Science Foundation of China (Grant nos. 41602239 and 41572209), the Natural Science Foundation of Jiangsu (Grant no. BK20160861), the International Postdoctoral Exchange Fellowship Program (Grant no. 20150048), and the Fundamental Research Funds of the Central Universities (Grant no. 2016B05514).