Fractal model and Lattice Boltzmann Method for Characterization of Non-Darcy Flow in Rough Fractures

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.

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 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.
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 ).

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. Nevertheless, from an engineering point of view, the case of deviation less than 10% would be acceptable, which will not significantly impact the general Scientific RepoRts | 7:41380 | DOI: 10.1038/srep41380 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.
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. 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.

Fractures build up by FDM.
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      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-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-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 eq 0 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: (1) 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. (2) 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. (3) 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