## Abstract

The irregular morphology of single rock fracture significantly influences subsurface fluid flow and gives rise to a complex and unsteady flow state that typically cannot be appropriately described using simple laws. Yet the fluid flow in rough fractures of underground rock is poorly understood. Here we present a numerical method and experimental measurements to probe the effect of fracture roughness on the properties of fluid flow in fractured rock. We develop a series of fracture models with various degrees of roughness characterized by fractal dimensions that are based on the Weierstrass–Mandelbrot fractal function. The Lattice Boltzmann Method (LBM), a discrete numerical algorithm, is employed for characterizing the complex unsteady non-Darcy flow through the single rough fractures and validated by experimental observations under the same conditions. Comparison indicates that the LBM effectively characterizes the unsteady non-Darcy flow in single rough fractures. Our LBM model predicts experimental measurements of unsteady fluid flow through single rough fractures with great satisfactory, but significant deviation is obtained from the conventional cubic law, showing the superiority of LBM models of single rough fractures.

#### Yang Ju^{1,2,3} , Qingang Zhang^{2} , Jiangtao Zheng^{1,4} , Chun Chang^{1} & Heping Xie^{5}

^{1}State Key Laboratory of Coal Resources & Safe Mining, China University of Mining & Technology, Beijing 100083, China. ^{2}State Key Laboratory for Geomechanics & Deep Underground Engineering, China University of Mining & Technology, Xuzhou 221116, China. ^{3}School of Mechanics and Civil Engineering, China University of Mining and Technology, Beijing 100083, China. ^{4}Department of Engineering Mechanics, School of Aerospace, Tsinghua University, Beijing 100084, China. ^{5}Key Laboratory of Energy Engineering Safety and Mechanics on Disasters, The Ministry of Education, Sichuan University, Chengdu 610065, China.

*Received: 18 August 2016, Accepted: 19 December 2016, Published: 01 February 2017*

© The Author(s) 2017

## Introduction

Natural rocks are generally composed of complex and heterogeneous fractures, which provide storage capacity and migration paths for oil, gas and water resources^{1,2}. Irregular morphology of rock fracture significantly complicates the fluid flow, resulting in unpredictable engineering processes for enhancing geothermal-reservoir mining, geological sequestration of carbon dioxide and groundwater remediation, etc. The irregular morphology of single rock fracture significantly influences subsurface fluid flow and gives rise to a complex and unsteady state that typically cannot be appropriately characterized using simple laws^{3,4,5,6}. In addition, coal-mine water-bursting disasters, coal-gas outburst accidents, dam disasters and rock-slope failures have shown to be closely related to fluid seepage, the dynamic evolution of rock fractures and coupled stress–fluid flow processes^{7,8}. A clear and detailed knowledge on the fluid flow and its interaction with stress in fractured media is critical when addressing the above engineering issues.

Being the basic element of the complex fracture network, a single fracture with its morphology controls the fluid flow initiation and development in the network. Some mechanical model approaches have been proposed to investigate the properties of fluid flows through single fracture and fracture networks, such as representative elementary volume (REV), discrete fracture network (DFN), hydrological-mechanical-chemical (HMC), thermos-hydro- mechanical (THM) approaches, and parallel plate and channel models^{1,9,10,11,12,13,14,15,16,17,18,19,20,21}. Liu *et al*.^{1} proposed the hydrological-mechanical-chemical (HMC) method to explain the enigmatically spontaneous changes in permeability that develop within single fracture in limestone under *in-situ* conditions. The parallel-plate model, which considers contact areas between matrix and fluid and artificial fractures, has been proposed to evaluate the effects of contact area and surface roughness on fluid flow in rock fractures^{22}.

To adapt and simulate the fluid flow in the dominant passageway, Tsang *et al*.^{19,20} presented a channel model for fluid flow through a tight fracture subjected to high stresses. However, the morphology of the contacted surfaces that stresses applied appeared to be so irregular that accurate definition on the structures of channel walls and the properties of fluid flow using mathematical or physics tools became extremely complicated. The fractal-dimension method (FDM) provides an effective way to accurately describe the fracture morphology, comparing traditional methods including the bump-height^{23} and the joint roughness coefficient (JRC)^{24,25}.

Barton *et al*.^{25} provided a revised method from a coupled joint behaviour model using the joint roughness coefficient (JRC) and verified that the properties of seepage flow were dominated by the morphology and connectivity of the passageway formed by the untouched fracture surface. However, the bump-height method of measuring every point’s bump-height in rock fracture is extremely difficult to apply in engineering practices, and the JRC method only qualitatively characters the ten known fracture types. By contrast, the fractal-dimension method can be used to quantitatively describe the fracture morphology with a much wider application. To the best of our knowledge, very few numerical studies have applied the fractal dimension method to estimate the effect of irregular fracture surface on the permeability and fluid velocity field.

The common cubic law (CCL), based on the smooth parallel-plate assumption that the aperture changes can result in a change of conductivity as much as three orders of magnitude at moderate compressive stress levels, has been widely applied to the analysis of seepage-flow behaviours in rock masses. The CCL has also been used to investigate the properties of fluid flow and the mechanisms of hydraulic-mechanical coupling in fractured rock^{1,9,10,11,12,13,16}. In CCL modelling, Darcy’s-law is implemented, which forms the basis of hydrogeology and is one of the most famous law that describes fluid flow through a porous medium^{26}. However, substantially differing from Darcy’s law, fluid flow in natural rough fractures can be influenced by a range of factors^{27} including surface roughness, fluid-matrix interface area, aperture, connectivity^{28,29,30}, unit width flux and hydraulic gradient. A comprehensive and detailed study on the fluid flow through natural rough fractures is needed.

Comparing to CCL of computational fluid dynamics (CFD), the Lattice Boltzmann Method (LBM) provides a powerful technique for modelling single/multiple phase flow in porous and fractured media with complex geometries^{27,31,32,33,34,35}. It is a kinetic-based mesoscopic approach that bridges the micro- and macro-scale, offering distinctive advantages in simulation fidelity and computational efficiency^{36}. The LBM has been widely applied to study fluid flow in porous media. Ju *et al*.^{37} presented the dynamic methane flow and distribution at microscale in porous sandstones subjected to force-induced deformation through LBM, and the method effectiveness in complex porous structure was validated by experimental observations.

Pazdniakou and Adler^{38} and Gao *et al*.^{39} used the LBM to investigate the dynamic permeabilities of porous media and the multicomponent fluid-flow in complex porous media, respectively. Fan and Zheng^{40} studied the seepage flow in a complex and rough fracture network using LBM. However, limited to an effective tool to describe the fracture morphology, little work has been published on the effects of fracture roughness on flow properties in single fractured rock.

The main objective of this study is to investigate the fluid flow in single rough fractures and the effect of irregular morphology of fractured rock by combining fractal-dimension method and LBM. The fractal governing function was embed for generating single rough fracture models with fractal dimension (*D*) varying from 1.0 to 1.5 and LBM for fluid (water) flow. Modelling results were validated by experimental measurements under the same conditions. The accuracy and efficiency of this numerical method with considering the non-Darcy flow field were analyzed and discussed with respect to the velocity-field distributions and equivalent permeabilities in the fractured models.

## Results

### Velocity distribution field

Figure 1 illustrates the velocity distribution of water over the entire fracture space, with detailed structural information from the fractal model of *D* = 1.5. To further investigate the influence of surface roughness and quantify the modelling results by experimental measurement^{7}, as shown in Fig. 1, we evenly selected five cross sections (A–E) along the flow pathway of a single rough fracture model with a fractal dimension varying from 1.0 to 1.5. Each cross section includes 4,000 lattice points (100 × 40), from which 14 representative points (marked by the black dots in Fig. 1) were symmetrically selected to display the velocity distribution across the section.

**Figure 1: Schematic representation of water velocity distribution and the five cross sections selected in the single rough fracture model of D = 1.5.** Water flows along x-axis direction (marked by the yellow arrows) perpendicular to y-z plane from the inlet (x = 0) to the outlet at x = 20 cm. The Inset image shows an example of the cross-section, in which the blue colour indicates matrix, the yellow indicates fracture, the red line indicates the model boundary and the 14 black dots indicate the symmetrically selected lattice points for further analysis on velocity distribution. The legend depicts the velocity magnitudes, in which blue indicates the minimum value and red represents the maximum number.

Figure 2 shows the distribution of water velocity over the five cross sections in fracture models with varying fractal dimensions (*D* = 1.0 to 1.5). Higher velocities in the centre of each cross-section were observed, with a decreasing trend from the centre to the both ends. These results show that in a single fracture with a constant fractal dimension, the variation of water velocity at different locations is small (see Fig. 2), even if the local roughness of the fracture is different. In the smooth flat fracture (*D* = 1.0), the flow velocity remained unchanged in the five different cross-sections.

**Figure 2: Velocity distribution of water flow along the 14 selected lattice points in the five selected cross sections of the single rough fracture models with varied fractal dimensions.** From (**a**–**f**), the fractal dimension D is equal to 1.0, 1.1, 1.2, 1.3, 1.4 and 1.5, respectively.

## Permeability of single rough fracture

Table 1 shows the permeabilities in fractures with varying fractal dimension as determined from the cubic law (Equation (3)) (*k*_{f}), the experiment measurement (*k*_{e}) and the LBM simulation (*k*_{0}).

**Table 1: Permeabilities as determined from the cubic law (k _{ f }), the LBM simulations (k_{ 0 }) and the experimental measurements (k_{ e })^{7}.**

## Discussion

Through the integral calculation of all the points over the entire cross section (y-z plane) of fracture with various fractal dimensions, the average velocity of water in different five cross-sections is obtained. Figure 3 illustrates the average velocity deviation between the LBM simulation and the experimental measurements over five cross-sections in various fracture models. The deviation between the numerical and experimental measurements is less than 10% for fractures with *D* = 1.0 to 1.4 (see Fig. 3). However, the deviation increases up to 30% for fracture with *D* = 1.5. The possible reason might be the discrepancy between the LBM numerical model of rock based on self-compiled programs and the physical cells of fractured rock.

**Figure 3: Experimental and simulated flow velocities over cross-sections of the stated fractal dimensions (D).** The D values are provided in the figure. The lines represent modelling results and the symbols represent experimental data. ∆: average velocity deviation between the modelling and experiments for each fractal dimension.

Nevertheless, from an engineering point of view, the case of deviation less than 10% would be acceptable, which will not significantly impact the general trend that the velocity evolves. The very rough fracture surface significantly influences the fluid velocity-field distribution. However, for a single fracture with constant fractal dimension, although the rough structure of the selected path segment was different, the average flow velocity did not change significantly, implying that the average flow velocity was independent of its local structural morphology.

Figure 4 shows the linear correlation between the average water velocity over the entire fracture space in single rough fracture models, as simulated using the LBM. Furthermore, a comparison between the rough fractures with various fractal dimensions suggests that as the fractal dimension increases — that is, as the fracture roughness increases — the average velocity of the flow in any segment decreases, as does the mean velocity of the water flow through the entire path of the fracture.

**Figure 4: The linear regression between velocity-field distribution of water and the various fractal dimensions.**

The deviation between *k*_{0} and *k*_{f} increases significantly with fractal dimension and exceeds 30% when *D* = 1.5 (see Table 1). Meanwhile, it is noted that the deviation between *k*_{0} and *k*_{f} is less than 5% when *D* = 1.0. Same deviation trend applied to *k*_{e} and *k*_{f}, but with much more significant increase in *D* = 1.5. The value keeps less than 15% as the fractal dimension is smaller than 1.4. The possible reason might result from the discrepancy between the LBM numerical model of rock based on self-compiled programs and the physical cells of fractured rock. In summary, it can be concluded that the permeability (*k*_{0}) decreases with increasing in the fractal dimension. Meanwhile, the results for the non-Darcy flow obtained using the LBM approach deviated significantly from the results obtained from CCL, indicating its inconsistency and incapability for describing and representing the complex flow behaviours in the fractal models.

The equivalent permeability coefficients (*k*_{f}, *k*_{e} and *k*_{0}) of water flow varying with the various fractal dimensions in the fracture models are plotted in Fig. 5.

**Figure 5: Variation in the equivalent-permeability coefficients (k _{f}, k_{e}, and k_{0}) of water flows in fracture models with varying fractal dimensions.**

One can easily formulate the following linear relationship between the equivalent-permeability coefficients (*k*_{0}) of a single fracture and the fractal dimension *D* of its rough fracture. These results suggest that the fractal equivalent permeability (*k*) decreases linearly as the fractal dimension of the rough structure (that is, roughness) increases, except for the case of the cubic law, where *k*_{f}, is constant. We found that the LBM simulation results have a good consistency with the experimental measurement. Therefore, it seems to be an effective way to quantitatively characterize the spatial distribution of flow velocity, permeability, and the influence of the roughness on the fluid flow behaviour in the single rough fracture models with various fractal dimensions.

## Methods

### Fractures build up by FDM

A series of single rough fracture models with different fractal dimensions were constructed using the Weierstrass–Mandelbrot function, as implemented in self-programming functions. The Weierstrass–Mandelbrot function is formulated as refs 41, 42

where *b* is a real number greater than 1, ∅_{n} is any angle and D ∈ (1, 2) is the fractal dimension. The fractal governing function, *C(t*) is then the real part of *W(t*)^{42}:

Considering the flow surface (see Fig. 6) along the fracture depth, we implemented the Weierstrass–Mandelbrot function and the physical cells established in previous study^{21} to build up the single rough fracture models with various fractal dimensions.

**Figure 6: Single rough fracture models with varying fractal dimensions (1.0–1.5) representing different surface roughness, and the dimensions of flow surface that water flows through.** D = 1.0 quantifies a smooth flat fracture, and the fracture roughness increases with D values. The fractal depth and width are 5 mm and 2 mm separately, and the total area of the flow surface (A-A) is 10 mm^{2}.

Figure 6 shows the single rough fracture models with various fractal dimensions, including the magnified inserts showing the detailed structure. The scale of the fractal model is 200 mm long, 100 mm wide and 5 mm thick in x, y, z axis, respectively. The scale and morphology of the fracture were kept identical to the cell used in the experiments. Further methodological details with experimental measurements and results can be found in Ju *et al*.^{7}.

### Velocity field and permeability in single rough fractures

Permeability refers to the ability of a fluid flows through the fractured or porous rock. Permeability is typically given as a function of the fracture aperture, *b*_{0}, under the conditions of isothermal and laminar flow between two parallel glass plates^{31}:

In this section, water velocity field and fracture permeability are investigated through LBM. Understanding the correlation between fracture morphology and permeability is thus important for accurately evaluating reservoir recovery and production rates. For that purpose, we adopted the distribution functions of flow velocity (equations (6, 7, 8, 9, 10, 11)) and the equation for permeability (equation (15)) to determine the permeability *k* of a fluid flowing through the single rough fracture. To simulate and analyse the effects of the rough surface on the fluid-flow behaviour in our models, the physical units including fluid pressure field *p*, macroscopic fracture aperture *L* and kinematic viscosity of fluid *v* were first transformed to the lattice units before determining the velocity-field distributions (Equations (12, 13, 14)). To enhance the accuracy of the simulation in the current context of fractured rocks and to reduce the computational time, we adopted the D3Q19 model^{43} to discretize the velocity at each lattice, and the single-relaxation-time Bhatnagar–Gross–Krook (BGK) approximation (equation (4)) was used to determine the movement and collision of fluid particles, which can be expressed as

where *f(r, ξ, t*) refers to the velocity distribution, which is a function of the spatial position vector *r*, velocity vector ξ, and time *t*. By discretizing the left-hand terms of equation (4) in time and space and replacing the right-hand term of equation (4) by a first-order rectangle approximation, we can convert the equation to

where *τ* = *τ*_{0}/*δ*_{t} is the dimensionless relaxation time and *δ*_{t}*F*_{a}(*r, t*) refers to the external force term. In the D3Q19 model, the distribution function at the equilibrium state is defined as

where *ω*_{a} is the weight coefficient, *c*_{s} is the sound velocity of a lattice and *c* is the lattice speed. *ρ* is the fluid density and *u* is the fluid velocity, which can be determined by equations (10) and (11):

Before determining the velocity-field distributions in our LBM simulation, the macroscopic parameters of the physical units including fluid pressure field *p*, fracture aperture *L* and kinematic viscosity of fluid *v* were first transformed to the lattice units, which are determined by equations (12), (13) and (14):

where L is the length of a lattice, N is the lattice number and τ is the relaxation time.

After obtaining *u*, the water permeability in different fracture models can be calculated by Darcy equation. Because the Reynolds number of the flow with the experimental measurement is lower than the critical value 2000, the viscous force prevails, indicating that the water flow in the rough fracture represents laminar flow^{7}. Thus, assuming the fluid flow meets the condition of laminar flow in the representative micro-scale segment, the rock permeability can be determined as

where *k* denotes the permeability of the medium, *x* refers to the direction of flow, −*d*_{p}/*d*_{x} is the pressure gradient along the flow direction, *μ* is the water viscosity and *U* is the average velocity per unit area.

## Boundary and initial conditions

The single rough fracture models were generated with a gridding size of 2000 × 1000 × 50 pixels in the LBM modelling, representing 0.2 × 0.1 × 0.005 m in physical size. The relevant parameters and boundary conditions used in the numerical simulation were identical to the experimental set-up and as follows:

- The density and viscosity of water is referred as 998.2
*kg/m*^{3}and 0.001003*Pa*·*s*, respectively^{7}. In order to make the simulation straightforward, we postulate the fluid flow within the fracture models as single phase flow. - The left boundary of the model was set as inlet, which was defined as a constant pressure boundary at 490 Pa. The right boundary of the model was set as outlet under atmospheric pressure (see Fig. 2). The initial velocity of the flow field was 0 m/s. The other parts of the model, with the exception of the fractal fracture, were set as ‘bounce-back boundaries’, indicating that the evolution of the particles was considered as head-on collisions of two particles.
Convergence: The simulation convergence was controlled by mesh generation, two- particle collision patterns, fluid property and iteration steps. The mesh resolution was set as 1 pixel to ensure convergence in the relatively small fractures. The modelling was stopped and the convergence results exported at iteration steps exceeded 8000 and the standard deviation of the average energy was less than 10

^{−4}.

## Experimental setups

To verify the effective of the LBM simulation, a series of single rough fracture physical cells with varying roughness were produced using the Weierstrass–Mandelbrot function and transparent polymethyl methacrylate (PMMA) material. A high-speed video camera was employed to record the fluid flow through the entire single rough fracture with a constant hydraulic pressure. The properties of fluid flow varying with the fracture roughness and the influences of the rough surface were analyzed. More details on the seepage experiments can be referred to Ju *et al*.^{7}.

#### Acknowledgements

The authors gratefully acknowledge the financial support of the National Natural Science Foundation of China (Grants 51374213 and 51674251), National Natural Science Fund for Distinguished Young Scholars of China (Grant 51125017), State Key Research Development Program of China (Grant 2016YFC0600705), Fund for Innovative Research and Development Group Program of Jiangsu Province (Grant 2014–27), Science Fund for Creative Research Groups of the National Natural Science Foundation of China (Grant 51421003), China Postdoctoral Science Foundation Funded Project and the Priority Academic Program Development of Jiangsu Higher Education Institutions (PAPD 2014).

#### Author Contributions

Yang Ju and Qingang Zhang designed and conducted the experimental tests and numerical simulation, analyzed the data, and wrote the paper; Jiangtao Zheng and Chun Chang discussed the experimental and numerical data, and improved figures; Heping Xie reviewed the experimental data and modified the paper. All authors discussed the results and commented on the paper.

#### Additional Information

- Competing financial interests: The authors declare no competing financial interests.
- Correspondence and requests for materials should be addressed to Y.J. (email: juy@cumtb.edu.cn or yju@icloud.com)
- How to cite this article: Ju, Y. et al. Fractal model and Lattice Boltzmann Method for Characterization of Non-Darcy Flow in Rough Fractures. Sci. Rep. 7, 41380; doi: 10.1038/srep41380 (2017).
- Publisher's note: Allaboutshale.com remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

## References

*1. Liu, J. et al. A fully-coupled hydrological–mechanical–chemical model for fracture sealing and preferential opening. Int. J. Rock Mech. Min. Sci. 43, 23–36 (2006).*

*2. Tao, Q., Ghassemi, A. & Ehlig Economides, C. A. A fully coupled method to model fracture permeability change in naturally fractured reservoirs. Int. J. Rock Mech. Min. Sci. 48, 259–268 (2011).*

*3. Dou, Z. & Zhou, Z. Lattice Boltzmann simulation of solute transport in a single rough fracture. Water Sci. Eng. 7, 277–287 (2014).*

*4. Alexis, D. A., Karpyn, Z. T., Ertekin, T. & Crandall, D. Fracture permeability and relative permeability of coal and their dependence on stress conditions. J. Unconv. Oil Gas Resour. 10, 1–10 (2015).*

*5. Natarajan, N. & Suresh Kumar, G. Evolution of fracture permeability due to co-colloidal bacterial transport in a coupled fracture-skin-matrix system. Geosci. Front. 3, 503–514 (2012).*

*6. Babadagli, T., Raza, S., Ren, X. & Develi, K. Effect of surface roughness and lithology on the water–gas and water–oil relative permeability ratios of oil-wet single fractures. Int. J. Multiphase Flow. 75, 68–81 (2015).*

*7. Ju, Y. et al. An experimental investigation on the mechanism of fluid flow through single rough fracture of rock. Sci. China, Tech. Sci. (China) 56, 2070–2080 (2013).*

*8. Wang, W., Li, X., Lin, B. & Zhai, C. Pulsating hydraulic fracturing technology in low permeability coal seams. Int. J. Min. Sci. Technol. 25, 681–685 (2015).*

*9. Baghbanan, A. & Jing, L. Stress effects on permeability in a fractured rock mass with correlated fracture length and aperture. Int. J. Rock Mech. Min. Sci. 45, 1320–1334 (2008).*

*10. Li, B., Jiang, Y., Koyama, T., Jing, L. & Tanabashi, Y. Experimental study of the hydro-mechanical behavior of rock joints using a parallel-plate model containing contact areas and artificial fractures. Int. J. Rock Mech. Min. Sci. 45, 362–375 (2008).*

*11. Baghbanan, A. & Jing, L. Hydraulic properties of fractured rock masses with correlated fracture length and aperture. Int. J. Rock Mech. Min. Sci. 44, 704–719 (2007).*

*12. Hudson, J., Stephansson, O. & Andersson, J. Guidance on numerical modelling of thermo-hydro- mechanical coupled processes for performance assessment of radioactive waste repositories. Int. J. Rock Mech. Min. Sci. 42, 850–870 (2005).*

*13. Min, K., Rutqvist, J., Tsang, C. & Jing, L. Stress-dependent permeability of fractured rock masses: a numerical study. Int. J. Rock Mech. Min. Sci. 41, 1191–1210 (2004).*

*14. Polak, A. et al. In Elsevier Geo-Engineering Book Series Vol. Vol. 2 (ed Stephanson, Ove) 721–726 (Elsevier, 2004).*

*15. Rutqvist, J. & Stephansson, O. The role of hydromechanical coupling in fractured rock engineering. Hydrol. J. 11, 7–40 (2003).*

*16. Xiao, Y. X., Lee, C. F. & Wang, S. J. Assessment of an equivalent porous medium for coupled stress and fluid flow in fractured rock. Int. J. Rock Mech. Min. Sci. 36, 871–881 (1999).*

*17. Hestir, K. & Long, J. C. S. Analytical expressions for the permeability of random two-dimensional Poisson fracture networks based on regular lattice percolation and equivalent media theories. J. Geophys. Res. (USA) 95, 21565–21581 (1990).*

*18. Kranz, R. L., Frankel, A. D., Engelder, T. & Scholz, C. H. Permeability of whole and jointed barre granite. Int. J. Rock Mech. Min. Sci. 16, 225–234 (1979).*

*19. Tsang, Y. W. & Tsang, C. F. Channel model of flow through fractured media. Water Resour. Res. 23, 467–479 (1987).*

*20. Tsang, Y. W. & Witherspoon, P. A. The dependence of fracture mechanical and fluid flow properties on fracture roughness and sample size. J. Geophys. Res. (USA) 88, 2359–2366 (1983).*

*21. Zhang, Q., Ju, Y., Gong, W., Zhang, L. & Sun, H. Numerical simulations of seepage flow in rough single rock fractures. Petroleum 1, 200–205 (2015).*

*22. Zhang, C., Yu, G. & Zhang, C. Rock matrix-fractured media model for heterogeneous and fractured coal bed. Trans. Nonferrous Met. Soc. China. 21, 621–625 (2011).*

*23. Pinardi, K. et al. Effect of bump height on the strain variation during the thermal cycling test of ACA flip-chip joints. IEEE Trans. Compon. Packag. Technol. 23, 447–451 (2000).*

*24. Barton, N. & Choubey, V. The shear strength of rock joints in theory and practice. Rock Mech. 10, 1–54 (1977).*

*25. Barton, N., Bandis, S. & Bakhtar, K. Strength, deformation and conductivity coupling of rock joints. Int. J. Rock Mech. Min. Sci. 22, 121–140 (1985).*

*26. Neuman, S. P. Transient flow of groundwater to wells in multiple-aquifer systems (1968).*

*27. Dou, Z., Zhou, Z. & Sleep, B. E. Influence of wettability on interfacial area during immiscible liquid invasion into a 3D self-affine rough fracture: Lattice Boltzmann simulations. Adv. Water Resour. 61, 1–11 (2013).*

*28. Crandall, D., Bromhal, G. & Karpyn, Z. T. Numerical simulations examining the relationship between wall-roughness and fluid flow in rock fractures. Int. J. Rock Mech. Min. Sci. 47, 784–796 (2010).*

*29. Cai, J., Yu, B., Zou, M. & Mei, M. Fractal analysis of surface roughness of particles in porous media. Chin. Phys. Lett. (China) 27, 024705 (2010).*

*30. Qian, J., Chen, Z., Zhan, H. & Guan, H. Experimental study of the effect of roughness and Reynolds number on fluid flow in rough-walled single fractures: A check of local cubic law. Hydrol. Processes. 25, 614–622 (2011).*

*31. Eker, E. & Akin, S. Lattice Boltzmann simulation of fluid flow in synthetic fractures. Transp. Porous Media. 65, 363–384 (2006).*

*32. Yan, Y. & Koplik, J. Flow of power-law fluids in self-affine fracture channels. Phys. Rev. E, Stat. Nonlinear Soft Matter Phys. (USA) 77, 036315 (2008).*

*33. Yao, J., Wang, C., Yang, Y., Hu, R. & Wang, X. The construction of carbonate digital rock with hybrid superposition method. J. Pet. Sci. Eng. 110, 263–267 (2013).*

*34. Yin, P. & Zhao, G. Numerical study of two-phase fluid distributions in fractured porous media. Int. J. Numer. Anal. Methods Geomech. 39, 1188–1211 (2015).*

*35. Dong, B., Yan, Y. Y. & Li, W. Z. LBM simulation of viscous fingering phenomenon in lmmiscible displacement of two fluids in porous media. Transp. Porous Media. 88, 293–314 (2011).*

*36. Li, Q. et al. Lattice Boltzmann methods for multiphase flow and phase-change heat transfer. Prog. Energy Combust. Sci. 52, 62–105 (2016).*

*37. Ju, Y., Wang, J., Gao, F. & Xie, H. Lattice-Boltzmann simulation of microscale CH4 flow in porous rock subject to force-induced deformation. Chin. Sci. Bull. 59, 3292–3303 (2014).*

*38. Pazdniakou, A. & Adler, P. M. Dynamic permeability of porous media by the lattice Boltzmann method. Adv. Water Resour. 62, Part B, 292–302 (2013).*

*39. Gao, J., Xing, H., Tian, Z. & Muhlhaus, H. Lattice Boltzmann modeling and evaluation of fluid flow in heterogeneous porous media involving multiple matrix constituents. Comput Geosci. 62, 198–207 (2014).*

*40. Fan, H. & Zheng, H. MRT-LBM-based numerical simulation of seepage flow through fractal fracture networks. Sci. China, Tech. Sci. (China) 56, 3115–3122 (2013).*

*41. Mandelbrot, B. B. The fractal geometry of nature. (W H Freeman, 1982).*

*42. Xie, H. Fractals in Rock Mechanics. 36–70 (Science press, 1993).*

*43. He, Y., Wang, Y. & Li, Q. Lattice Boltzmann method: theory and applications. (Science Press, 2009).*

This work is licensed under a Creative Commons Attribution 4.0 International License. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in the credit line; if the material is not included under the Creative Commons license, users will need to obtain permission from the license holder to reproduce the material. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/

© The Author(s) 2017