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 nonDarcy flow through the single rough fractures and validated by experimental observations under the same conditions. Comparison indicates that the LBM effectively characterizes the unsteady nonDarcy 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 geothermalreservoir 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, coalmine waterbursting disasters, coalgas outburst accidents, dam disasters and rockslope 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), hydrologicalmechanicalchemical (HMC), thermoshydro 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 hydrologicalmechanicalchemical (HMC) method to explain the enigmatically spontaneous changes in permeability that develop within single fracture in limestone under insitu conditions. The parallelplate 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 fractaldimension method (FDM) provides an effective way to accurately describe the fracture morphology, comparing traditional methods including the bumpheight^{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 bumpheight method of measuring every point’s bumpheight 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 fractaldimension 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 parallelplate 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 seepageflow behaviours in rock masses. The CCL has also been used to investigate the properties of fluid flow and the mechanisms of hydraulicmechanical coupling in fractured rock^{1,9,10,11,12,13,16}. In CCL modelling, Darcy’slaw 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, fluidmatrix 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 kineticbased mesoscopic approach that bridges the micro and macroscale, 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 forceinduced 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 fluidflow 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 fractaldimension 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 nonDarcy flow field were analyzed and discussed with respect to the velocityfield 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 xaxis direction (marked by the yellow arrows) perpendicular to yz plane from the inlet (x = 0) to the outlet at x = 20 cm. The Inset image shows an example of the crosssection, 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 crosssection 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 crosssections.
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 (yz plane) of fracture with various fractal dimensions, the average velocity of water in different five crosssections is obtained. Figure 3 illustrates the average velocity deviation between the LBM simulation and the experimental measurements over five crosssections 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 selfcompiled programs and the physical cells of fractured rock.
Figure 3: Experimental and simulated flow velocities over crosssections 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 velocityfield 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 velocityfield 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 selfcompiled 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 nonDarcy 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 equivalentpermeability 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 equivalentpermeability 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 selfprogramming 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 (AA) 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 fluidflow 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 velocityfield 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 singlerelaxationtime 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 lefthand terms of equation (4) in time and space and replacing the righthand term of equation (4) by a firstorder 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 velocityfield 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 microscale 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 setup 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 ‘bounceback boundaries’, indicating that the evolution of the particles was considered as headon 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 highspeed 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 NonDarcy 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 fullycoupled 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 cocolloidal bacterial transport in a coupled fractureskinmatrix 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 oilwet 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 hydromechanical behavior of rock joints using a parallelplate 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 thermohydro 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. Stressdependent 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 GeoEngineering 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 twodimensional 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 matrixfractured 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 flipchip 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 multipleaquifer systems (1968).
27. Dou, Z., Zhou, Z. & Sleep, B. E. Influence of wettability on interfacial area during immiscible liquid invasion into a 3D selfaffine 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 wallroughness 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 roughwalled 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 powerlaw fluids in selfaffine 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 twophase 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 phasechange heat transfer. Prog. Energy Combust. Sci. 52, 62–105 (2016).
37. Ju, Y., Wang, J., Gao, F. & Xie, H. LatticeBoltzmann simulation of microscale CH4 flow in porous rock subject to forceinduced 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. MRTLBMbased 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