Research Article  Open Access
PoreScale Modeling of MixingInduced Reaction Transport through a Single SelfAffine Fracture
Abstract
This porescale modeling study in single selfaffine fractures showed that the heterogeneous flow field had a significant influence on the mixinginduced reaction transport. We generated the single selfaffine fracture by the successive random additions (SRA) technique. The porescale model was developed by coupling the NavierStoke equation (NSE) and advectiondiffusion equation with reaction (ADER). Eddies were captured in the selfaffine fracture due to the increasing Reynolds number and the sudden expansion of aperture. The fluxweighted breakthrough curves (BTCs) of reaction product showed the typical nonFickian characteristics (i.e., “early arrival” and “heavy tail”). It was found that the reactant was involved in eddies and then reacted inside the eddycontrolled domain. Consequently, eddies played a significant role in delaying the mass exchange process between the eddycontrolled 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 mixinginduced 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 macroscopicscale approach, the advectiondispersion reaction equation (ADRE), based on a perfect mixing assumption, commonly overestimates the reaction rate and tends to generate a less reliable prediction. Many studies [3–9] have shown that the ADRE is generally incapable of describing the local mixing phenomenon and the neglected porescale incomplete mixing is a primary reason for the overestimated reaction rate. Therefore, the mixinginduced 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 selfaffine structures for length scales ranging from microns to meters [14]. The selfaffine structures can be generated by the successive random additions (SRA) technique [15, 16]. Therefore, in this study, the flow field in the single selfaffine 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 fluxaveraged breakthrough curves (BTCs) of a passive tracer. Due to the heterogeneity of geological formation, the BTCs might be nonGaussianshaped and could exhibit the nonFickian characteristics (i.e., “early arrival” and “heavy tail”). Although applying macroscopic advectiondispersion equation (ADE) to fit the nonFickian BTCs could yield a nonnegligible error, the nonFickian 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 nonFickian 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 mixinginduced reaction transport. However, the influence of eddies on the mixing process in a selfaffine fracture is still an open topic.
The main objective of this work was to investigate mixinginduced reaction transport through single selfaffine fractures by using a porescale model. The selfaffine fracture was generated by the SRA technique. The flow field in the single selfaffine fracture was numerically solved by NSE. The nonreaction transport and the mixinginduced reaction transport were simulated by the ADRE. The flow field and eddies formation were first examined under four different Reynolds numbers (Re). Sequentially, the nonFickian BTCs for nonreaction transport and the mixinginduced 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 SelfAffine Fracture Generation
Following the seminal work of Mandelbrot and Pignoni [21], the height of the selfaffine 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 selfaffine 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 selfaffine fracture wall (e.g., the SRA and the Fourier transformation). In this study, we used the SRA technique [15, 16] to generate the selfaffine fracture wall. The application of the SRA can be found in our previous works [23, 24]. From the generated selfaffine fracture wall to the single selfaffine fracture, we used the shear displacement model [25] to construct the single selfaffine 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.
(a)
(b)
2.2. Flow Model
The flow field in a single selfaffine 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/m^{3} and =1.002×10^{−3} Pas. The selfaffine fracture walls were considered as the nonslip boundaries. The steadystate 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 steadystate 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 selfaffine 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 advectiondiffusion 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 fluxweight breakthrough curves (BTCs) can be obtained by a ratio of effluent solute mass to fluid mass:By normalizing fluxweight 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 onetoone 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 secondorder 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 wellmixing 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 selfaffine 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 mixingcontrolled.
3. Problem Setup
This study first considered the cases with four different pressure gradients to simulate the flow field in the generated single selfaffine fracture, hereafter, denoted as flow field 1–flow field 4 (FF#1FF#4). For each FF#1FF#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 selfaffine 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 flowreaction 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 mixinginduced 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 selfaffine 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 selfaffine 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 outofplane 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 eddycontrolled domain at the fulllength 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 [28–32] 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 eddycontrolled 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 fluxweighted BTCs of compound A for the nonreaction transport were calculated under the different flow fields in the single selfaffine 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 nonFickian 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 nonFickian 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 eddycontrolled domain, the mass exchange between the main flow channel and eddycontrolled 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 mixinginduced reaction leads to the decreasing concentration peak in BTCs for the same Re and Pe. However, this mixinginduced 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 nonFickian tail in the BTCs of the reaction product. First, the reaction product was transformed at the outside of the eddycontrolled domain and then captured by eddies. Second, the reaction product was transformed inside the eddycontrolled domain, which required that the reactant was involved in the eddycontrolled domain. Comparing Figures 4 and 7, it was found that the instantaneous transformation was inside the eddycontrolled domain due to the involved reactant. This is important because that the effective chemical kinetics of reaction transport inside the eddycontrolled 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 eddycontrolled 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 mixinginduced reactions has recently won much attention and been reported in many literatures [36–41]. 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 selfaffine 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 mixinginduced 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 selfaffine 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#1FF#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 mixinginduced nonreaction transport as shown in Figure 8. This suggests that the chemical kinetics of the mixinginduced reaction transport can be qualitatively analyzed by the temporal evolution of the reactants mixing process.
5. Summary and Conclusion
To model the mixinginduced reaction transport through a single selfaffine fracture, we generated the single selfaffine fracture by the SRA technique and a porescale model was developed by coupling the NSE and the ADRE. This allowed us to directly capture fundamental physical processes of flow, nonreaction, and mixinginduced reaction without consideration of parameterizations at larger macroscopic scales. The selfaffine 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 mixinginduced reaction showed the nonFickian 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 eddycontrolled 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 eddycontrolled domain. This may result in the different effective chemical kinetics of reaction transport between inside eddycontrolled 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 mixinginduced 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 porescale analysis performed in this study highlights the critical role of the heterogeneous flow field on the mixinginduced reaction transport through the selfaffine 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).
References
 Y. Jin, X. Li, M. Zhao, X. Liu, and H. Li, “A mathematical model of fluid flow in tight porous media based on fractal assumptions,” International Journal of Heat and Mass Transfer, vol. 108, pp. 1078–1088, 2017. View at: Publisher Site  Google Scholar
 J. Cai, W. Wei, X. Hu, and D. A. Wood, “Electrical conductivity models in saturated porous media: A review,” EarthScience Reviews, vol. 171, pp. 419–433, 2017. View at: Publisher Site  Google Scholar
 J. Cao and P. K. Kitanidis, “Porescale dilution of conservative solutes: An example,” Water Resources Research, vol. 34, no. 8, pp. 1941–1949, 1998. View at: Publisher Site  Google Scholar
 V. Kapoor, C. T. Jafvert, and D. A. Lyn, “Experimental study of a bimolecular reaction in Poiseuille flow,” Water Resources Research, vol. 34, no. 8, pp. 1997–2004, 1998. View at: Publisher Site  Google Scholar
 D. S. Raje and V. Kapoor, “Experimental study of bimolecular reaction kinetics in porous media,” Environmental Science & Technology, vol. 34, no. 7, pp. 1234–1239, 2000. View at: Publisher Site  Google Scholar
 O. A. Cirpka and A. J. Valocchi, “Twodimensional concentration distribution for mixingcontrolled bioreactive transport in steady state,” Advances in Water Resources, vol. 30, no. 67, pp. 1668–1679, 2007. View at: Publisher Site  Google Scholar
 G. Chiogna and A. Bellin, “Analytical solution for reactive solute transport considering incomplete mixing within a reference elementary volume,” Water Resources Research, vol. 49, no. 5, pp. 2589–2600, 2013. View at: Publisher Site  Google Scholar
 M. Dentz and P. J. De Barros, “Mixingscale dependent dispersion for transport in heterogeneous flows,” Journal of Fluid Mechanics, vol. 777, pp. 178–195, 2015. View at: Publisher Site  Google Scholar
 J. Qian, H. Zhan, Y. Zhang, P. Sun, and Y. Liu, “Numerical Simulation and Experimental Study of Bimolecular Reactive Transport in Porous Media,” Transport in Porous Media, vol. 109, no. 3, pp. 727–746, 2015. View at: Publisher Site  Google Scholar
 Y. W. Tsang, “The Effect of Tortuosity on Fluid Flow Through a Single Fracture,” Water Resources Research, vol. 20, no. 9, pp. 1209–1215, 1984. View at: Publisher Site  Google Scholar
 R. W. Zimmerman and G. S. Bodvarsson, “Hydraulic conductivity of rock fractures,” Transport in Porous Media, vol. 23, no. 1, pp. 1–30, 1996. View at: Google Scholar
 D. J. Brush and N. R. Thomson, “Fluid flow in synthetic roughwalled fractures: NavierStokes, Stokes, and local cubic law simulations,” Water Resources Research, vol. 39, no. 4, pp. SBH51–SBH515, 2003. View at: Google Scholar
 L. Wang and M. B. Cardenas, “An efficient quasi3D particle trackingbased approach for transport through fractures with application to dynamic dispersion calculation,” Journal of Contaminant Hydrology, vol. 179, pp. 47–54, 2015. View at: Publisher Site  Google Scholar
 B. Berkowitz, “Characterizing flow and transport in fractured geological media: a review,” Advances in Water Resources, vol. 25, no. 8–12, pp. 861–884, 2002. View at: Publisher Site  Google Scholar
 R. F. Voss, “Fractals in nature: from characterization to simulation,” in The Science of Fractal Images, Springer, 1988. View at: Google Scholar
 H.H. Liu, G. S. Bodvarsson, S. Lu, and F. J. Molz, “A corrected and generalized successive random additions algorithm for simulating fractional levy motions,” Mathematical Geology, vol. 36, no. 3, pp. 361–378, 2004. View at: Publisher Site  Google Scholar
 R. Haggerty, S. A. McKenna, and L. C. Meigs, “On the latetime behavior of tracer test breakthrough curves,” Water Resources Research, vol. 36, no. 12, pp. 3467–3479, 2000. View at: Publisher Site  Google Scholar
 Z. Dou, Z.F. Zhou, and J.G. Wang, “Threedimensional analysis of spreading and mixing of miscible compound in heterogeneous variableaperture fracture,” Water Science and Engineering, vol. 9, no. 4, pp. 293–299, 2016. View at: Publisher Site  Google Scholar
 J. Luo and O. A. Cirpka, “How well do mean breakthrough curves predict mixingcontrolled reactive transport?” Water Resources Research, vol. 47, no. 2, 2011. View at: Google Scholar
 Z. Dou, Z. Zhou, Q. Guo, and X. Lu, “Effects of flow fluctuation on dilution and spreading in a selfaffine fracture,” Water Environment Research, vol. 89, no. 8, pp. 752–762, 2017. View at: Publisher Site  Google Scholar
 B. B. Mandelbrot and R. Pignoni, The Fractal Geometry of Nature, 1983.
 J. M. Boffa, C. Allain, and J. P. Hulin, “Experimental analysis of fracture rugosity in granular and compact rocks,” EPJ Applied Physics , vol. 2, no. 3, pp. 281–289, 1998. View at: Publisher Site  Google Scholar
 Z. Dou, Z. Zhou, and B. E. Sleep, “Influence of wettability on interfacial area during immiscible liquid invasion into a 3D selfaffine rough fracture: lattice Boltzmann simulations,” Advances in Water Resources, vol. 61, pp. 1–11, 2013. View at: Publisher Site  Google Scholar
 Z. Dou and Z.F. Zhou, “Lattice Boltzmann simulation of solute transport in a single rough fracture,” Water Science and Engineering, vol. 7, no. 3, pp. 277–287, 2014. View at: Google Scholar
 J. S. Y. Wang, T. N. Narasimhan, and C. H. Scholz, “Aperture correlation of a fractal fracture,” Journal of Geophysical Research: Solid Earth, vol. 93, pp. 2216–2224, 1988. View at: Google Scholar
 H. Ralph, S. William, F. G. Herring, and D. Madura Jeffrey, General chemistry: Principles and modern applications, Prentice Hall, 2007.
 K. A. Connors, Chemical kinetics: the study of reaction rates in solution, John Wiley & Sons, 1990.
 J. Qian, M. Liang, Z. Chen, and H. Zhan, “Eddy correlations for water flow in a single fracture with abruptly changing aperture,” Hydrological Processes, vol. 26, no. 22, pp. 3369–3377, 2012. View at: Publisher Site  Google Scholar
 J. Bouquain, Y. Meheust, D. Bolster, and P. Davy, “The impact of inertial effects on solute dispersion in a channel with periodically varying aperture,” Physics of Fluids, vol. 24, 2012. View at: Google Scholar
 L. Zou, L. Jing, and V. Cvetkovic, “Roughness decomposition and nonlinear fluid flow in a single rock fracture,” International Journal of Rock Mechanics and Mining Sciences, vol. 75, pp. 102–118, 2015. View at: Publisher Site  Google Scholar
 S. Briggs, B. W. Karney, and B. E. Sleep, “Numerical modeling of the effects of roughness on flow and eddy formation in fractures,” Journal of Rock Mechanics and Geotechnical Engineering, 2016. View at: Google Scholar
 M. Wang, Y.F. Chen, G.W. Ma, J.Q. Zhou, and C.B. Zhou, “Influence of surface roughness on nonlinear flow behaviors in 3D selfaffine rough fractures: Lattice Boltzmann simulations,” Advances in Water Resources, vol. 96, pp. 373–388, 2016. View at: Publisher Site  Google Scholar
 B. F. Armaly, F. Durst, J. C. F. Pereira, and B. Schoenung, “Experimental and theoretical investigation of backwardfacing step flow,” Journal of Fluid Mechanics, vol. 127, pp. 473–496, 1983. View at: Publisher Site  Google Scholar
 B. B. Dykaar and P. K. Kitanidis, “Macrotransport of a biologically reacting solute through porous media,” Water Resources Research, vol. 32, no. 2, pp. 307–320, 1996. View at: Publisher Site  Google Scholar
 M. B. Cardenas, D. T. Slottke, R. A. Ketcham, and J. M. Sharp Jr., “NavierStokes flow and transport simulations using real fractures shows heavy tailing due to eddies,” Geophysical Research Letters, vol. 34, no. 14, 2007. View at: Google Scholar
 M. Rolle and P. K. Kitanidis, “Effects of compoundspecific dilution on transient transport and solute breakthrough: A porescale analysis,” Advances in Water Resources, vol. 71, pp. 186–199, 2014. View at: Publisher Site  Google Scholar
 O. A. Cirpka and P. K. Kitanidis, “Characterization of mixing and dilution in heterogeneous aquifers by means of local temporal moments,” Water Resources Research, vol. 36, no. 5, pp. 1221–1236, 2000. View at: Publisher Site  Google Scholar
 A. M. Tartakovsky, G. D. Tartakovsky, and T. D. Scheibe, “Effects of incomplete mixing on multicomponent reactive transport,” Advances in Water Resources, vol. 32, no. 11, pp. 1674–1679, 2009. View at: Publisher Site  Google Scholar
 M. Dentz, T. Le Borgne, A. Englert, and B. Bijeljic, “Mixing, spreading and reaction in heterogeneous media: A brief review,” Journal of Contaminant Hydrology, vol. 120121, no. C, pp. 1–17, 2011. View at: Publisher Site  Google Scholar
 E. Ballarini, S. Bauer, C. Eberhardt, and C. Beyer, “Evaluation of the role of heterogeneities on transverse mixing in benchscale tank experiments by numerical modeling,” Groundwater, vol. 52, no. 3, pp. 368–377, 2014. View at: Publisher Site  Google Scholar
 M. Rolle, C. Eberhardt, G. Chiogna, O. A. Cirpka, and P. Grathwohl, “Enhancement of dilution and transverse reactive mixing in porous media: Experiments and modelbased interpretation,” Journal of Contaminant Hydrology, vol. 110, no. 34, pp. 130–142, 2009. View at: Publisher Site  Google Scholar
 T. Le Borgne, M. Dentz, D. Bolster, J. Carrera, J.R. de Dreuzy, and P. Davy, “NonFickian mixing: Temporal evolution of the scalar dissipation rate in heterogeneous porous media,” Advances in Water Resources, vol. 33, no. 12, pp. 1468–1475, 2010. View at: Publisher Site  Google Scholar
 P. K. Kitanidis, “The concept of the Dilution Index,” Water Resources Research, vol. 30, no. 7, pp. 2011–2026, 1994. View at: Publisher Site  Google Scholar
Copyright
Copyright © 2018 Zhi Dou et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.