Lattice Boltzmann simulation of solute transport in a single rough fracture

Zhi DOU*, Zhi-fang ZHOU

School of Earth Sciences and Engineering, Hohai University, Nanjing 210098, P. R. China

Lattice Boltzmann simulation of solute transport in a single rough fracture

Zhi DOU*, Zhi-fang ZHOU

School of Earth Sciences and Engineering, Hohai University, Nanjing 210098, P. R. China

In this study, the lattice Boltzmann method (LBM) was used to simulate the solute transport in a single rough fracture. The self-affine rough fracture wall was generated with the successive random addition method. The ability of the developed LBM to simulate the solute transport was validated by Taylor dispersion. The effect of fluid velocity on the solute transport in a single rough fracture was investigated using the LBM. The breakthrough curves (BTCs) for continuous injection sources in rough fractures were analyzed and discussed with different Reynolds numbers (Re). The results show that the rough fracture wall leads to a large fluid velocity gradient across the aperture. Consequently, there is a broad distribution of the immobile region along the rough fracture wall. This distribution of the immobile region is very sensitive to the Re and fracture geometry, and the immobile region is enlarged with the increase of Re and roughness. The concentration of the solute front in the mobile region increases with the Re. Furthermore, the Re and roughness have significant effects on BTCs, and the slow solute molecule exchange between the mobile and immobile regions results in a long breakthrough tail for the rough fracture. This study also demonstrates that the developed LBM can be effective in studying the solute transport in a rough fracture.

solute transport; single rough fracture; Lattice Boltzmann method; self-affinity; breakthrough curve

1 Introduction

The importance of solute transport in fractured rock has long been recognized in many fields, such as water resources assessment, tracer testing at the field scale, and groundwater pollutant remedy. In particular, the process of non-aqueous phase liquid (NAPL) bioremediation involves immiscible and miscible two-phase flow (Dou et al. 2012). The rate of mass transfer between NAPLs and pure water is strongly dependent on the rate of solute transport. The solute transport in a fracture is a very complicated process and is influenced by a number of factors (Dou et al. 2013). For most studies on solute transport in a fracture, it is usually assumed that the fracture surface is smooth. However, in practice, the fracture wall is extremely rough and the roughness of the fracture surface has an important impact on solutetransport. Meanwhile, it is well known that the solute transport relies not only upon the roughness of the fracture surface but also on the flow field conditions. The flow field conditions are important for modeling many subsurface transport problems (Crandall et al. 2010; Cai et al. 2010). The variation of the aperture in a rough fracture causes a heterogeneous fluid velocity field that affects solute transport (Qian et al. 2011).

Due to the varying fluid velocity at the cross-section of a fracture, some of the solutes remain static in the immobile fluid regions, and other mobile solutes go through the middle cross-section of the fracture. The solute in the immobile fluid regions might result in a long breakthrough tail. The phenomenon in which solute moves quickly through the middle of the fracture but slowly along the fracture wall can be described as the film transport. This film transport leads to a high degree of solute dispersion and a long tail of the breakthrough curve (BTC) (Yeo 2001). The roughness of the fracture wall has a significant effect on the thickness of film and film transport (Or and Tuller 2000). Many field and laboratory studies on solute transport have been conducted to investigate the influence of roughness on solute transport. The traditional parallel plate model undergoes a serious challenge in describing the properties of the solute transport in a natural fracture with rough surface walls. Thompson (1991) presented a two-dimensional numerical simulation of solute transport in a fracture with fractal surface walls and found that the effective solute velocity is significantly less than the calculated fluid velocity due to the occurrence of different contact areas and resulting transverse dispersion in the fracture. Zhang (2000) studied the solute transport in heterogeneous soils and found that the heterogeneity of soils has a significant effect on solute molecular travel time, which refers to the time that a solute molecule takes to move from the inlet to the outlet. Auradou et al. (2001) experimentally studied the solute fronts in rough self-affine fractures. Their results revealed that the solute front varied with the displacement mechanism. Cardenas et al. (2009) simulated the solute transport in a rough asymmetric fracture and confirmed that the roughness of the fracture wall could lead to a long tail of the breakthrough curve. In this paper, the travel time refers to the time that the solute molecule takes to break through the whole fracture.

In recent years, the lattice Boltzmann method (LBM) has received broad recognition. The LBM is a powerful technique for the computational modeling of a wide variety of complex fluid flow problems, including single- and multi-phase flows in complex geometries. This method naturally accommodates a variety of boundary conditions. Some studies have reported that the LBM is capable of simulating the solute transport and solving the advection-dispersion equation. Zhang et al. (2002) used the LBM to solve the two-dimensional (2D) advection and anisotropic dispersion equation (AADE) based on the uniform flow field. A similar study was reported by Zhou (2009). However, research on solute transport in a self-affine rough fracture with a non-uniform flow field is still limited.

The main objective of this study was to investigate the influence of a non-uniform flowfield on solute transport in a single rough fracture. The self-affine rough fracture wall was generated by the successive random addition method. Two distribution functions of lattice Boltzmann equations were used to represent the flow and concentration fields, respectively. The capability of this method for simulating the solute transport in a rough fracture with consideration of the non-uniform flow field was examined by Taylor dispersion. The characteristics of solute transport and the corresponding BTCs are analyzed and discussed.

2 Method

2.1 Coupling flow and concentration fields with lattice Boltzmann method

Two distribution functions are used to represent the flow and concentration fields, respectively. The evolution equation of each distribution function satisfies the following lattice Boltzmann equation with a single relaxation time:

where fi(X ,t) and ci(X ,t) are the distribution functions for fluid particles and solute particles in the i direction with a velocity eiat position X and time t, respectively, where in the two-dimensional nine-speed (D2Q9) model adopted in this study, the nine possible particle velocities (e0,e1,…,e8) are given by

where c=ΔxΔt, with Δx being the lattice space and Δt being the time step; τ and τ′ are two different dimensionless relaxation times; and feqi(X ,t) and ceqi(X ,t) are the equilibrium distribution functions for the fluid and solute, respectively, defined as

where csis the lattice sound speed (lu/ts); ueqis the local equilibrium velocity vector; and iω is the weight, which, for the D2Q9 model, is given by

ρ and C are the macroscopic fluid density and solute concentration, respectively, and can respectively be expressed as

The fluid pressure p can be expressed as

Depending on the different dimensionless relaxation times (τ and τ′) in Eqs. (1) and (2), the corresponding kinematic viscosity and diffusion coefficient are obtained by

Using the Chapman-Enskog expansion, the Navier-Stokes equation can be recovered with second-order accuracy from the LBM for the fluid field (Eq. (1)). Similarly, the solute particle distribution function (Eq. (2)) can be recovered with the macroscopic advection-diffusion equation. Thus, the governing equations in the model that couples the flow and concentration fields are given by

In our simulations, any lattice node in the computational domain represents either a solid node or a fluid node. For the solid node, before the streaming, the bounce-back algorithm instead of the collision step is implemented with a non-slip wall boundary condition. All simulations were realized with the C++ code.

2.2 Self-affine rough fracture

An ideal self-affine profile can be generated using the simple algorithm suggested by Voss (1988). For a two-dimensional self-affine rough fracture, a self-affine set is invariant under the affinity transformation. Thus, the height of a self-affine rough profile along the horizontal direction (x direction), h(x), satisfies the following formula:

where H is the roughness or Hurst exponent, and λ is a scaling factor. There is ample evidence that the Hurst exponent of natural fractures in rock is H? 0.8 ± 0.05. Schmittbuhl et al. (1993) recorded the height of a granitic fault surface as a function of position along one-dimensional profiles. Their results revealed that the profiles followed an anisotropic scaling invariance (self-affinity) and the Hurst exponent was H≈ 0.85. Boffa et al. (1998) used high-resolution roughness measurements to study the fractured granite and basalt. They found that the fracture profiles for the granite and basalt displayed a self-affine geometry witha characteristic Hurst exponent H≈0.8. There are three methods commonly used to construct a synthetic self-affine fracture wall surface: the successive random addition method, the randomization of the Weierstrass function based on Mandelbrot, and the Fourier transformation. In this study, the successive random addition method was applied to generate rough fracture walls. This method is a straightforward way of generating self-affine surfaces and consists of three major steps. First, hi(x) are initialized with hi= 0 at points x1,x2,… ,xn, where hi(x) is the heightalong the x direction after the ith iteration, and n is the number of points in the x direction. Second, a standard Gaussian variable δh0iwith variance σ02= 1 is added to the initialized height to obtain hi+δh0i. Then, the mid-height between hi(x) and hi+1(x ) is estimated by interpolation:

and all the mid-heights for all points (x1,x2,… , xn) are renumbered from 1 to 2n-1. This process is repeated a great number of times and, at the nth iteration, the variance of the centered Gaussian variable hi+δh0iis equal to. Thus, the height of a self-affine rough profile along the x direction hi(x) is equal to the value of hi(x) at the nth iteration. Fig. 1 shows two self-affine rough fracture walls with different Hurst exponents.

Fig. 1 Self-affine fracture walls with different Hurst exponents

3 Results and Discussion

3.1 Taylor dispersion simulation

A simple demonstration of the capability of LBM for simulation of solute transport in a rough fracture is the simulation of Taylor dispersion (Fig. 2). Stockman et al. (1997) used the lattice-gas model to simulate Taylor dispersion:

where D is the dispersion coefficient, b is the aperture within a single fracture, Dmis the molecule diffusion coefficient, and v is the mean velocity of the flow. Taylor dispersion occurs as a result of the velocity variations across the fracture aperture due to diffusion. Here we used two regular lattice Boltzmann equations to demonstrate Taylor dispersion. Fig. 2 shows Taylor dispersion in a single fracture, where the red contours represent the fracture wall surfaces and C0is the initial solute concentration.

Fig. 2 Taylor dispersion in single fracture

A series of case studies with different apertures were conducted to examine the ability of LBM to simulate Taylor dispersion in a single fracture. The channel widths were 5 to 50 lu and the lengths were 100 to 1000 lu. The average velocities for all simulations were set to 0.007 lu/ts. The molecule diffusion coefficient was set to Dm= 0.0013 lu2/ts. Fig. 3 shows the analytical solution (Eq. (17)) and the dispersion coefficient estimated through fitting the breakthrough curves computed with LBM to solutions to the advection-dispersion equation (ADE). It can be seen that the results of LBM are consistent with the analytical solution.

3.2 Characteristics of solute transport in single rough fracture

3.2.1 Effect of fluid velocity on solute transport in single rough fracture

A series of simulations of solute transport in a single rough fracture with different fluid velocities were conducted with the continuous injection method. Two self-affine fracture walls with H = 0.8 and H = 0.7 for the upper and lower fracture walls, respectively, were generated. A computational domain of 550 lu × 130 lu representing a two-dimensional variable-aperture fracture was used. The flow direction was from left to right (Fig. 4). For a concentration field, the constant concentration and full developed boundary conditions were used for the inlet and outlet, respectively. For the flow field, the periodic boundary condition was used for the inlet and outlet, respectively. Fig. 4 shows the steady-state velocity field in a two-dimensional variable-aperture fracture, where the red line represents the fracture wall, the black arrow with a tail represents the velocity field, and the dashed blue line represents the parallel single fracture. It can be seen from Fig. 4 that the velocity field is non-uniform and varies with the apertures of the fracture, and the rough fracture wall leads to a large fluid velocity gradient across the aperture. Due to the limitation of computational capacity, the mean aperture b for the variable-aperture fracture was 63 lu, as shown in Fig. 4. It is noted that the mean aperture has a significant effect on solute transport. Thus, the mean aperture was constant for the oncoming single fractures in this study.

Fig. 3 Dispersion coefficients obtained from LBM and Eq. (17)

Fig. 4 Flow field in a single rough fracture

In order to investigate the effect of fluid velocity on solute transport in single rough fractures, two cases with different Re values (Re = 13 and Re = 86) were examined.

Fig. 5 shows the process of solute transport in the rough fracture for Re = 13, where the red contour represents the fracture wall. It can be seen from the figure that the roughness of a fracture wall has a significant effect on the pathway of solute transport. As shown in Fig. 4, the flow field in the rough fracture is strongly dependent on the geometry of the aperture. Solutes go preferentially through the middle region of the aperture due to the relatively fast local fluid velocity. The solute fronts are first observed in the middle region of the aperture and the front spreading pathway varies with the local fluid velocity due to the roughness. Thus, the fracture geometry has a significant effect on the mobile region for solute transport. As more solutes are injected, the concentration of the solute front decreases due to the dominant molecular diffusion. This may result in a shorter breakthrough time. After solutes break through the whole fracture, solutes continue to diffuse toward the rough fracture walls (Figs. 5(c) and 5(d)) due to the transverse concentration gradient. This concentration gradient induces a transverse solute flux that enhances macro-dispersion. However, some regions with little solute (immobile regions) are observed (Fig. 5(d)) and remain low-concentration. The occurrence of these immobile regions is due to two factors: fluid velocities in these regions are close to zero, associated with the low solute flux, and the concentration increases slowly due to the solute molecular diffusion. The solute transport in the immobile region is controlled by the molecular diffusion mechanism. In theory, the solute molecules will be diffusing throughout the injection time and will be in equilibrium at infinite injection time. However, in this study, if the outletconcentration difference for a period of 8 000 ts was less than 0.1%, it was assumed that the final steady state for solute transport had arrived. In fact, the equilibrium state of solute molecular diffusion in these regions will last a very long time in a field test, even with a small Re. Thus, the mobile and immobile regions form as shown in Fig. 5(d).

Fig. 5 Process of solute transport in rough fracture for Re = 13

Fig. 6 shows the process of solute transport in the rough fracture for Re = 86. In this case, the geometry of the aperture in the rough fracture is completely consistent with that in the case with Re = 13. It can be seen from Fig. 6 that the increasing fluid velocity has a significant effect on solute transport in the rough fracture. The solute concentration gradient rises due to the dominant advection mechanism. The solute front with a high concentration gradient spreads in the middle of the aperture where the fluid velocity is relatively high. Although the advection is dominant for solute transport, the solutes can diffuse toward the rough fracture wall where the solute concentration is relatively low. This slow diffusion process results in the increment of marco-dispersion, which is similar to the case with Re = 13. The flow field has a significant effect on the distribution of solute concentration. Comparing Fig. 6(d) to Fig. 5(d), it can be seen that the boundaries between the immobile and mobile regions in Fig. 6(d) become distinct. This is because the greater Re value associates with the stronger advection effect on solute transport. The slow process of solute molecular diffusion takes limited time to reach the final steady state. It can also be found that the area of the immobile region in Fig. (6) is greater than that in Fig. (5). Furthermore, the distribution of the immobile region along the rough fracture wall is sensitive to the fracture geometry. The lesser roughness leads to a smaller immobile region for solute transport in the rough fracture.

Fig. 6 Process of solute transport in rough fracture for Re = 86

3.2.2 Breakthrough curve

For the smooth parallel fracture with continuous injection of a constant concentration, the classical analytical solution with the longitude dispersion coefficient DLis as follows (Bodin et al. 2003):

with the following initial and boundary conditions:

where Cf(x ,t) is the solute concentration at position x and time t.

Fig. 7 shows the BTCs for different Re values in the rough and smooth fractures, where the dimensionless time T=t/T0, with T0being the total time. For the solute transport in the smooth fracture, the Re value is set to 13, which is the same as one value in the rough fracture. The analytical solutions are calculated by Eq. (18). The results of LBM are in good agreement with the analytical solutions. The maximum error between analytical solutions and results of LBM is less than 0.15%. Fig. 7 shows that not only the Re but also the roughness has a significant effect on BTCs. For the BTC, the long breakthrough time indicates the long travel time. As expected, the smaller Re leads to the longer travel time in the rough fracture. Moreover, the travel times for the cases with Re = 13 and Re = 86 are less than that in the smooth fracture for Re = 13. The travel time for the solute molecule strongly relies on the fluid velocity in the middle of the aperture for the rough fracture. Even though Re is constant for both rough and smooth fractures, the travel time is different. This is because the fluid velocity in the middle of the aperture for the rough fracture is greater than that for the smooth fracture due to the variation of the aperture in the rough fracture. The travel time is sensitive to the roughness.

Fig. 7 Breakthrough curves for different Re values and types of fracture

In addition, long tails at the end of BTCs for the cases of the rough fracture are observed. For the type of rough fracture in this study, a reasonable explanation for the long breakthrough tail, suggested by Raven et al. (1988), is that the diffusive exchange of solutes between the mobile and immobile regions occurs.

4 Conclusions

In this study, the developed LBM was used to simulate the solute transport in a self-affinerough fracture. The ability of the developed LBM to simulate the solute transport was examined through Taylor dispersion. This study demonstrated that the LBM was very effective in simulating the solute transport in a rough fracture. The main conclusions are as follows:

(1) Since the roughness led to the non-uniform fluid velocity field in the fracture, the pathways of solute transport were separated into the mobile and immobile regions. The mobile region was distributed along the central aperture in the rough fracture and the advection mechanism was dominant for the solute transport in the mobile region. The immobile region was observed near the rough fracture walls and the molecular diffusion mechanism was dominant for the solute transport in the immobile region.

(2) The distribution of the immobile region was sensitive to both the Re and fracture geometry, and the immobile region was enlarged with the increase of the Re and roughness.

(3) The concentration of the solute front in the mobile region increased with the Re. The concentration of the solute front decreased due to the dominant molecular diffusion. This may result in a shorter breakthrough time.

(4) The Re and roughness have significant effects on BTCs, and the slow solute molecule exchange between the mobile and immobile regions results in a long breakthrough tail for the single rough fracture.

The current study was limited to a two-dimensional single rough fracture. Three-dimensional fractures will be examined in our future research. In this study, the solute transport was studied in a single rough fracture with different Re values. Future studies will further consider the different roughness values in a single fracture.

Auradou, H., Hulin, J. P., and Roux, S. 2001. Experimental study of miscible displacement fronts in rough self-affine fractures. Physical Review E, 63(6), 066306. [doi:10.1103/PhysRevE.63.066306]

Bodin, J., Delay, F., and de Marsily, G. 2003. Solute transport in a single fracture with negligible matrix permeability: 2. Mathematical formalism. Hydrogeology Journal, 11(4), 434-454. [doi:10.1007/ s10040-003-0269-1]

Boffa, J. M., Allain, C., and Hulin, J. P. 1998. Experimental analysis of fracture rugosity in granular and compact rocks. The European Physical Journal-Applied Physics, 2(3), 281-289. [doi:10.1051/ epjap:1998194]

Cai, J. C., Yu, B. M., Zuo, M. Q., and Mei, M. F. 2010. Fractal analysis of surface roughness of particles in porous media. Chinese Physics Letters, 27(2), 024705. [doi:10.1088/0256-307x/27/2/024705]

Cardenas, M. B., Slottke, D. T., Ketcham, R. A., and Sharp, J. M. 2009. Effects of inertia and directionality on flow and transport in a rough asymmetric fracture. Journal of Geophysical Research: Solid Earth, 114(B6), B6204. [doi:10.1029/2009jb006336]

Crandall, D., Bromhal, G., and Karpyn, Z. T. 2010. Numerical simulations examining the relationship between wall-roughness and fluid flow in rock fractures. International Journal of Rock Mechanics and Mining Sciences, 47(5), 784-796. [doi:10.1016/j.ijrmms.2010.03.015]

Dou, Z., Zhou, Z. F., Huang, Y., and Wu, W. 2012. Numerical study of non-aqueous phase liquid transport in a single filled fracture by lattice boltzmann method. Journal of Hydrodynamics, Ser. B, 24(1), 130-137. [doi:10.1016/s1001-6058(11)60227-8]

Dou, Z., Zhou, Z., and Sleep, B. E. 2013. Influence of wettability on interfacial area during immiscible liquidinvasion into a 3D self-affine rough fracture: Lattice Boltzmann simulations. Advances in Water Resources, 61, 1-11. [doi:10.1016/j.advwatres.2013.08.007]

Or, D., and Tuller, M. 2000. Flow in unsaturated fractured porous media: Hydraulic conductivity of rough surfaces. Water Resources Research, 36(5), 1165-1177. [doi:10.1029/2000WR900020]

Qian, J. Z., Chen, Z., Zhan, H. B., and Guan, H. C. 2011. Experimental study of the effect of roughness and Reynolds number on fluid flow in rough-walled single fractures: A check of local cubic law. Hydrological Processes, 25(4), 614-622. [doi:10.1002/hyp.7849]

Raven, K. G., Novakowski, K. S., and Lapcevic, P. A. 1988. Interpretation of field tracer tests of a single fracture using a transient solute storage model. Water Resources Research, 24(12), 2019-2032. [doi:10.1029/WR024i012p02019]

Schmittbuhl, J., Gentier, S., and Roux, S. 1993. Field measurements of the roughness of fault surfaces. Geophysical Research Letters, 20(8), 639-641. [doi:10.1029/93gl00170]

Stockman, H. W., Li, C. H., and Wilson, J. L. 1997. A lattice-gas and lattice Boltzmann study of mixing at continuous fracture junctions: Importance of boundary conditions. Geophysical Research Letters, 24(12), 1515-1518. [doi:10.1029/97gl51471]

Thompson, M. E. 1991. Numerical simulation of solute transport in rough fractures. Journal of Geophysical Research: Solid Earth, 96(B3), 4157-4166. [doi:10.1029/90jb02385].

Voss, R. F. 1988. Fractals in nature: From characterization to simulation. The Science of Fractal Images, 21-70. New York: Springer. [doi:10.1007/978-1-4612-3784-6_1]

Yeo, W. 2001. Effect of fracture roughness on solute transport. Geosciences Journal, 5(2), 145-151. [doi:10.1007/bf02910419].

Zhang, R. D. 2000. Generalized transfer function model for solute transport in heterogeneous soils. Soil Science Society of America Journal, 64(5), 1595-1602. [doi:10.2136/sssaj2000.6451595x].

Zhang, X. X., Crawford, J. W., Glyn Bengough, A., and Young, L. M. 2002. On boundary conditions in the lattice Boltzmann model for advection and anisotropic dispersion equation. Advances in Water Resources, 25(6), 601-609. [doi:10.1016/S0309-1708(02)00027-1].

Zhou, J. G. 2009. A lattice Boltzmann method for solute transport. International Journal for Numerical Methods in Fluids, 61(8), 848-863. [doi:10.1002/fld.1978].

(Edited by Yan LEI)

This work was supported by the National Natural Science Foundation of China (Grants No. 51079043, 41172204, and 51109139) and the Natural Science Foundation of Jiangsu Province (Grant No. BK2011110).

*Corresponding author (e-mail: douz@hhu.edu.cn)

Received May 28, 2013; accepted Dec. 6, 2013