Numerical Simulation of the Nonlinear Flow Properties in Self-Affine Aperture-Based Fractures

Abstract

In order to study the effect of fracture geometry on the nonlinear flow properties in aperture-based fractures, a fractal model based on the self-affinity is proposed to characterize the three-dimensional geometry of rough-walled fractures. By solving the N–S (Navier–Stokes) equation directly, the relationships between the Forchheimer-flow characteristics, fractal dimension, and standard deviation of the aperture have been obtained. The Forchheimer equation is validated to describe the nonlinear relationship between flow rate and pressure gradient. For lower flow rate, the influence of the fractal dimension almost can be ignored, but the linear coefficient increases and the hydraulic aperture decreases with increasing standard deviation of the aperture, respectively. For larger flow rate, the nonlinear coefficient increases with the growth of the standard deviation of the aperture and fractal dimension. Thus, an empirical relationship between the nonlinear coefficient, fractal dimension, and standard deviation of aperture is proposed. In addition, the critical Reynolds number decreases with the increase of the standard deviation of the aperture and the fractal dimension, and the numerical results are generally consistent with the experimental data.

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

Fractures are prevalent among the natural rock masses under geological action. Compared with the intact rock, the permeability of fractures is much larger, and fractures become the dominant channels for fluid flow [1, 2]. Fracture-dominated flow plays an important role in many practical situations such as seepage control for fractured media [3, 4], hazardous waste isolation [5], slope engineering [6, 7], underground tunneling [8, 9], and many geological disasters.

Many efforts have been taken to investigate the hydraulic properties of fractures. Snow [10] has conceptualized the rough-walled fracture as the smooth parallel plate model, and the famous cubic law was derived. However, the surface roughness and aperture distribution of natural fractures are random and irregular [11, 12]. Consequently, there is no absolute smooth fracture in natural rock masses and it is difficult to satisfy the establishment conditions of cubic law, which may lead to a large deviation in the prediction of fracture permeability. Based on the laboratory tests and numerical analysis, many researchers [13–21] found that the relationship between flow rate and pressure gradient in rough-walled fractures is nonlinear deviating from the cubic law.

To investigate the relationship between the nonlinear flow properties and the roughness of the fracture surface, Zhang and Tian [22] carried out numerical simulation of a single fracture with different roughness and tortuousness, and the results showed a certain deviation from the modified cubic law resulted by the flow tortuosity. Chen et al. [16] studied the relationship between the hydraulic aperture and Forchheimer equation’s nonlinear coefficient by conducting flow tests under different confining pressures on twelve groups of granite cracks with different roughness. Yin et al. [19] analyzed the effects of the shearing process under a series of normal loads on nonlinear flow behavior in 3D rough-walled fractures with different roughness by the shear-flow test. The abovementioned investigations on the impact of roughness on the nonlinear flow properties all use the JRC (Joint Roughness Coefficient) to represent the fracture surface roughness. However, the value of the JRC is obtained based on observation experiences; thus, the roughness of fracture surfaces cannot be accurately and quantitatively represented.

Some other influence factors of rough fractures on the nonlinear flow properties were also considered. For example, Xia et al. [23] found that different contact conditions and the root-mean square height of fractures had an evident impact on apparent transmissivity based on laboratory observations. Tsang [24] explored the impact of tortuous of flow path on flow behavior through experiments and found that the smaller the aperture distribution, the greater the influence of tortuosity. Xiong et al. [17] designed a saturated seepage test of fracture under low flow rate and evaluated the role of roughness and aperture in affecting the nonlinear flow properties. It is not difficult to find that there are many factors affecting the nonlinear flow properties of rough fractures. However, the effect of the standard deviation of the fracture aperture on the nonlinear flow properties is rarely discussed.

The main objective of this study is to study the nonlinear flow behavior in self-affine aperture-based rock fractures based on the fractal theory and Navier–Stokes equations. Based on the fractal theory, the rough-walled fracture is reconstructed by fractional Brownian motion, and the surface roughness and aperture distribution can be characterized by fractal dimension. For rough-walled fractures with different surface roughness and the standard deviation of the aperture, the nonlinear relationship between flow rate and pressure gradient can be well described by the Forchheimer equation. The mathematical relationship between the nonlinear coefficient, fractal dimensions, and the standard deviation of the aperture is also quantified.

2. Geometrical Model of Rough Fractures

Previous studies [18, 25–27] have shown that the fracture roughness can be described by self-affinity, which can be simulated using fractional Brownian motion (fBm). The height of rough fracture surfaces is described by a continuous and single random function Z(x). The stationary increment [Z(x)−Z(x + Δ)] over the distance Δ follows a Gaussian distribution with mean zero and variance δ². The self-affinity relating to fBm obeys the following expressions:$$\begin{array}{c}\text{(1)}& Zx-Zx+\lambda \Delta =0,{\delta }_{\lambda \Delta }^2={\lambda }^{2H}{\delta }_{\Delta }^2\Rightarrow {\delta }_{\lambda \Delta }={\lambda }^H{\delta }_{\Delta },\\ {\delta }_{\lambda \Delta }^2={Zx-Zx+\lambda \Delta }^2,\\ {\delta }_{\Delta }^2={Zx-Zx+\Delta }^2,\end{array}$$where <·> is the mathematical expectation, x donates the coordinate component, and H represents the Hurst exponent varying from 0 to 1 and associated with the fractal dimension D by D = 3-H. λ is a constant, ${\delta }_{\lambda \mathrm{\Delta }}^2$ and ${\delta }_{\mathrm{\Delta }}^2$ are the variance corresponding to the height variation of fracture surface with the distance of λΔ and Δ, respectively, and δ is the standard deviation of the aperture.

To construct the geometrical model of rough-walled fractures, in present study, the successive random addition method (SRAM) [28–30] is used to generate the fracture surfaces. The generated square area is shown in Figure 1, and the specific steps are as follows:

(1) In the given single square region, the initial random values of the four corners, which are labeled as number 1, satisfy the Gaussian distribution $N0,{\delta }_0^2$

(2) The center point of the square and midpoints of each side are marked by number 2, and their height are the average value of the four corners initial height and the average value of two points of each side, respectively; at the same time, the random values from $N0,{\delta }_1^2$ should be added into all points, in which$$\begin{array}{}\text{(2)}& {\delta }_1^2=\frac{{\delta }_0^2}{2^{2H}}1-2^{2H-2}.\end{array}$$

(3) Step (2) is repeated for each newly generated square, and the random values from $N0,{\delta }_n^2$ should be added to all heights as well until $2^n\times 2^n$ squares are created at iteration$$\begin{array}{}\text{(3)}& {\delta }_n^2=\frac{{\delta }_{n-1}^2}{2^{2H}}=\frac{{\delta }_0^2}{{2^{2H}}^n}1-2^{2H-2}.\end{array}$$

(4) No new square is inserted again, but we keep adding random values from $N0,{\delta }_j^2$ to all points, in which

$$\begin{array}{}\text{(4)}& {\delta }_j^2=\frac{{\delta }_{j-1}^2}{2^{2H}}=\frac{{\delta }_0^2}{{2^{2H}}^j}1-2^{2H-2},\end{array}$$where j = n + 1, n + 2, …, NM; NM is large enough, so δ_NM/δ₀ is negligible

[figures omitted; refer to PDF]

The upper and lower surfaces of rough fracture are generated using the SRAM. The height of low surface is Z₁(x,y), and the height of the upper surface can be calculated by the following formula:$$\begin{array}{}\text{(5)}& Z_{2}x,y=Z_{1}x+\Delta x,y+\Delta y+u,\end{array}$$where u is the mean aperture between the upper and lower surfaces. The distribution function of the fracture aperture can be expressed as$$\begin{array}{}\text{(6)}& bx,y=\begin{array}{cc}{Z}_{2}x,y-Z_{1}x,y,& \text{if\hspace{0.17em}}{Z}_{2}x,y>Z_{1}x,y,\\ 0,& \text{else}.\end{array}\end{array}$$

To demonstrate the reliability of the algorithm used to characterize the fracture surface geometry, Figure 2 shows the fracture aperture distribution diagram of four groups with different roughness. The side length of square L = 40 mm, δ = 0.15 mm, and fractal dimensions are 2.1, 2.2, 2.3, and 2.4, respectively. It can be seen that, with the increase of D, the maximum values of the aperture increase and the minimum values of the aperture decrease; the aperture distribution is more scattered, the variation of adjacent apertures is more fluctuated, and the degree of roughness is greater.

[figures omitted; refer to PDF]

3. Theory of Flow Properties in Fractures

3.1. Cubic Law

Based on the smooth parallel plate model, the flow behaviour in a single fracture should obey the cubic law [10]:$$\begin{array}{}\text{(7)}& Q=-\frac{w{b}_h^3}{12\mu }\nabla p,\end{array}$$where Q is the volumetric flow rate per unit time, $b_h$ is the hydraulic aperture of the fracture, $w$ is the fracture width, μ is the dynamic viscosity, and $\nabla p=p_{\text{in}}-p_{\text{out}}/l$ is the pressure gradient between the inlet and outlet.

3.2. Forchheimer Equation

For larger flow rate, the flow behaviour of rough fractures is nonlinear due to the larger inertial effect. Thus, instead of equation (7), the Forchheimer equation [14, 15, 17, 19] is commonly applied to describe the nonlinear flow properties:$$\begin{array}{c}\text{(8)}& -\nabla P=AQ+B{Q}^2,\text{(9)}& A=\frac{\mu }{k{A}_h}=\frac{12\mu }{w{b}_h^3},B=\frac{\beta \rho }{w^2{b}_h^2},\end{array}$$where A and B are the linear coefficient and nonlinear coefficient, respectively; k is the fracture permeability; β is the non-Darcy flow inertial coefficient which depends on the geometric characteristics of fracture surfaces [16, 31]; and ρ is the fluid density. When the flow rate is small, the inertia force is much less than the viscosity force; in other words, the quadratic term (BQ²) is much less than the linear term (AQ), and equation (8) can be reduced to the cubic law.

3.3. Distinguishing the Flow Regime

In order to quantify the nonlinear flow behaviour of fractures, the Reynolds number Re is calculated as follows:$$\begin{array}{}\text{(10)}& \mathrm{Re}=\frac{\rho v{b}_h}{\mu }=\frac{\rho Q}{\mu w},\end{array}$$where $v$ is the average velocity of the fracture inlet. The Reynolds number Re represents the ratio of inertial force to viscous force in fracture flow, and the inertial effect is stronger with larger Re; thus, it is easier to enter the nonlinear flow regime.

Simultaneously, the non-Darcy effect factor E is defined according to the Forchheimer equation:$$\begin{array}{}\text{(11)}& E=\frac{{\text{BQ}}^2}{\text{AQ}+{\text{BQ}}^2}.\end{array}$$

The physical meaning of E is the proportion of pressure gradient caused by nonlinear flow in the total pressure gradient, and it is a dimensionless coefficient ranging from 0 to 1 representing the degree of nonlinear flow. The greater the E is, the stronger the nonlinear effect is. At present, E = 0.1 [15, 23, 32], and combined with equations (10) and (11), the critical Reynolds number Re_c for the transition from linear to nonlinear flow of fracture flow yields$$\begin{array}{}\text{(12)}& R{e}_c=\frac{A\rho }{9B\mu w}.\end{array}$$

The smaller the Re_c is, the more significant the inertia effect is and the easier it is to be the nonlinear flow state.

4. Numerical Simulation

4.1. Governing Equation

For the incompressible Newtonian fluid with a constant viscosity coefficient flowing in the rough fractures, the fluid motion is governed by the N–S equation and the continuity equation:$$\begin{array}{}\text{(13)}& \rho \mathbf{u}\cdot \nabla \mathbf{u}-\mu {\nabla }^2\mathbf{u}=-\nabla P,\text{(14)}& \nabla \mathbf{u}=0,\end{array}$$where u and $\nabla$ are the velocity vector and Hamiltonian operator, respectively. In this study, the fluid density is 997.1 (kg/m³), and the dynamic viscosity is 0.894 ×10⁻³ (Pa s) for water at 25°C. Since the flow rate is small and generally does not exceed the boundary of laminar flow, the laminar interface in finite element software COMSOL Multiphysics is selected for the solution [33].

4.2. Numerical Procedure

As shown in Figure 3, at first, three-dimensional rough fracture surfaces are generated based on the SRAM, and then, the solid fracture model is constructed by the Boolean operation. Next, the solid model is meshing to determine the microstructure of the computing model. Finally, the software COMSOL is employed to obtain a series of data of flow rate and pressure gradient which is fitted by the Forchheimer equation.

[figure omitted; refer to PDF]

Considering the solution accuracy, calculation cost, and errors caused by different element sizes, the tetrahedral element size is set as 0.25 mm, and the number of grid elements ranges from 400 thousands to 600 thousands. The size of the fracture model is 40 mm × 40 mm, and the mean aperture u is 0.8 mm; the geometric parameters including the standard deviation of the aperture and fractal dimension are shown in Table 1. The left side of the fracture model is specified as inflow boundary, the range of Q is 3.586 ×10⁻⁸∼3.228 × 10⁻⁶·m³/s, the corresponding Re ranges from 1 to 90 according to equation (10). The right side of the fracture model is defined as outflow boundary and the pressure is set as zero, and the rest of the boundaries are impervious.

Table 1

The values for linear coefficient A and nonlinear coefficient B.

Standard deviation of aperture (mm)	D = 2.1		D = 2.2		D = 2.3		D = 2.4		D = 2.5
	A	B	A	B	A	B	A	B	A	B
δ = 0.09	5.5111	1.1802	5.4536	1.2779	5.5118	1.3642	5.4648	1.4179	5.5147	1.5342
δ = 0.12	5.6875	1.2289	5.6428	1.3294	5.6168	1.4194	5.6985	1.5478	5.6647	1.6495
δ = 0.15	5.9625	1.2942	5.9763	1.4131	5.9316	1.5536	5.8904	1.6491	5.8862	1.8116
δ = 0.18	6.2710	1.3749	6.1622	1.5075	6.2558	1.6758	6.0919	1.7666	6.1335	2.0007
δ = 0.21	6.5353	1.4817	6.4933	1.5996	6.4495	1.7724	6.4531	1.9774	6.4546	2.2197

Notes: A (10⁸ kg·s⁻¹·m⁻⁵); B (10¹⁴ kg·m⁻⁸).

5. Results and Discussion

5.1. Relationship between Flow Rate and Pressure Gradient

The relationships between flow rate and pressure gradient of five groups are shown in Figure 4. The fitting linear coefficient A and nonlinear coefficient B of the Forchheimer equation are presented in Table 1, and the coefficients of determination R² are both greater than 0.99, which indicates the nonlinear relationship between the flow rate and pressure gradient in self-affine aperture-based fractures can be well described by the Forchheimer equation.

[figures omitted; refer to PDF]

In Figure 4, with the increment of flowrate, the deviation between Forchheimer curves and cubic law increases drastically. This deviation is also strengthened with the increment of D and δ, which indicates that flow resistance will be greater for the rougher surface. Zeng and Grigg [32] suggested that the Forchheimer-flow properties were mainly caused by inertial force; in other words, the inertial term $\rho \mathbf{u}\cdot \nabla \mathbf{u}$ in equation (13) was the main factor for the nonlinear behavior.

In order to analyze the primary factor of nonlinear flow behavior, Figure 5 compared the five velocity sections along the x direction when D = 2.5 and Q = 5.3796 × 10⁻⁷·m³/s of two fracture models with different δ (0.09 mm and 0.21 mm). It can be seen that the velocity distribution is more scattered in the fracture with more heterogeneous aperture distribution, and the local velocity becomes relatively larger. As a result, the inertia effect of flow is enhanced, and the local energy dissipation is increased, which leads to the fluid entering nonlinear state.

[figures omitted; refer to PDF]

5.2. Analysis of Permeability of Rough Fractures

Based on equation (9), the linear coefficient A is negatively correlated with the permeability of rough fractures. As shown in Table 1, A increases with the δ increase, indicating that the permeability will be lower when the distribution of the aperture is more discrete and irrelevant for a larger δ. The same phenomenon can also be found in Figure 6, the hydraulic aperture $b_h$ decreases apparently with the δ increases, and it is always less than the mean aperture 0.8 mm; therefore, the permeability of rough fractures is less than that of smooth fractures. Moreover, with the increment of δ, the flow paths become more tortuous and the permeability becomes smaller.

[figure omitted; refer to PDF]

However, D has no obvious effect on A except for a slight fluctuation as shown in Table 1. In other words, when the flow is small, the fractal dimension has little effect on the flow capacity of rough cracks The section diagrams of streamline distribution are presented in Figure 7, for δ = 0.21 mm and Q = 2.152 × 10⁻⁶·m³/s with D being 2.1, 2.3, and 2.5, y = 20 mm, and x ranging from 25 mm to 40 mm. For same δ, the streamline distribution is similar while the roughness of the fracture surfaces is greater for larger D, and there are some blank areas at the rough bulge. This is because the fluid flow in the fracture tends to the region with lower resistance and will bypass the high resistance area to form a dominant channel. Therefore, the effective flow space is substantially equal and the flow capacity does not change significantly. When the flow rate continues to increase, eddy flow will appear in these blank areas, resulting in accelerated energy consumption and the permeability of fractures being reduced [34].

[figures omitted; refer to PDF]

5.3. Analysis of Forchheimer-Flow Characteristics

The B and β represent the degree of evolution of Forchheimer-flow properties, and the nonlinear flow is stronger with a greater value of B and β. In Table 1, the values of B increase with larger δ and D, which indicates that the degree of flow nonlinearity is more drastic with the rougher surfaces and more heterogeneous aperture distribution.

According to equation (9), B is proportional to β but inversely proportional to the square of $b_h$, so it is difficult for B to fully reflect the impact of fracture morphology on flow properties of rough fractures. Instead, Figure 8 shows the relationship between δ, D, and β. It is observed that β increases linearly with δ and D. Therefore, the two-parameter expression is proposed:$$\begin{array}{}\text{(15)}& \beta =\eta \delta D+m\delta +nD,\end{array}$$where η, m, and n are empirical parameters. On the basis of nonlinear Levenberg–Marquardt algorithm, the fitting parameters η, m, and n are 517.5, −967, and 49.39, respectively, and the coefficient of determination R² is 0.9879, indicating that this expression could well explain the effect of morphology parameters on Forchheimer-flow properties. Combining with equation (12), the empirical model of nonlinear coefficient B is obtained:$$\begin{array}{}\text{(16)}& B=\frac{\rho }{w^2{b}_h^2}517.5\delta D-967\delta +49.39D.\end{array}$$

[figure omitted; refer to PDF]

It should be noted that when the value of $b_h$ decreases with the increment of δ, the growth of B is more significant. Therefore, the effect of δ on the nonlinear flow behavior of the rough fractures is greater than that of D.

5.4. Critical Reynolds Number

Re_c can be used to evaluate the flow regime of fractures, representing the Reynolds number of the fluid flow going from the linear state to nonlinear state. When Re < Re_c, the flow behavior satisfied the cubic law on account of the inertia effect being very small and the viscous effect having a dominant role of fluid flow in fractures; when Re > Re_c, the flow behavior was controlled by the inertia effect, appearing Forchheimer-flow properties. As shown in Figure 9, for the same condition, Re_c has a decreasing trend with the increment of δ and D, indicating that the fluid flow can more easily develop into nonlinear Forchheimer-flow resulting from the rougher surfaces and more heterogeneous aperture distribution.

[figure omitted; refer to PDF]

Previous studies [13, 14, 35–38] on Re_c of flow in rock fractures are listed in the Table 2. In this study, for rough fractures of u = 0.8 mm, δ being 0.09∼0.21 mm and D ranging from 2.1 to 2.5, the calculated values of Re_c are between 9.01 and 14.47, which are almost consistent with the previous researches.

Table 2

Literature summaries about critical Reynolds number of flow in rock fractures.

Authors	Materials	Aperture ($\mu m$)	Rec	Study method
Konzuk and Kueper	Tensile fractures of limestone	Mean = 381	2.8∼14.3	Laboratory test
Javadi et al.	Tensile fractures of granite	—	0.001∼25	Laboratory test
Zhang and Nemcik	Tensile fractures of sandstone	6.14∼18.95	3.5∼24.8	Laboratory test
Zimmerman et al.	3D fractures	Mean = 149	1∼10.20	Numerical simulation
Chen et al.	Sedimentary and intrusive rocks	—	25∼66	Field test
Qian et al.	Artificial parallel plate model	1000∼2500	245∼759	Laboratory test
This paper	3D rough fractures	Mean = 800	9.01∼14.47	Numerical simulation

6. Conclusions

In order to estimate the relationships between the nonlinear flow properties and self-affine fracture geometry, a systematic procedure with respect to the nonlinear flow analysis for the three-dimensional rough-walled fractures model is proposed based on the fractal characteristics, in which the roughness of fracture surfaces is characterized by the fractal dimension and the fractal model is generated by the SRAM. The numerical simulation of nonlinear flow analysis in self-affine aperture-based fractures is presented based on the fractal theory and N–S equation, in which the spatial variation of fracture geometry including roughness and aperture distribution is characterized by the fractal dimension and N–S equation is solved by the software COMSOL. The main conclusions are as follows:

(1) The larger the fractal dimension, the greater the fracture surface roughness and the lower the correlation of the local aperture. The larger the standard deviation of the aperture, the greater the fluctuation of fracture surface and the more heterogeneous the aperture distribution.

(2) The nonlinear relationship between flow rate and pressure gradient in self-affine aperture-based fractures can be well described by the Forchheimer equation. The flow rate, standard deviation of the aperture, and fractal dimension can improve the deviation of pressure gradient from the cubic law.

(3) With the increase of the standard deviation of the aperture, the linear coefficient of the Forchheimer equation is larger, and the hydraulic aperture becomes smaller. When the flow rate is small, fractal dimension has little effect on the permeability of fracture.

(4) The nonlinear coefficient of the Forchheimer equation increases with the increment of the standard deviation of the aperture and fractal dimension, and the empirical expression between the nonlinear coefficient, the standard deviation of aperture, and the fractal dimension is proposed. The standard deviation of the aperture has greater impact on the nonlinear flow behavior than the fractal dimension. The critical Reynolds number decreases with the increase of standard deviation of the aperture and fractal dimension, and its measured range is 9.01∼14.47, which is almost consistent with the previous results.

Acknowledgments

The financial supports from the National Natural Science Foundation of China (Nos. 42077243 and 51709207), Natural Science Foundation of Hubei Province (No. 2018CFB631), and Visiting Researcher Fund Program of State Key Laboratory of Water Resources and Hydropower Engineering Science (2019SGG04) are gratefully acknowledged.

References

[1] C. Yao, Y. Shao, J. Yang, F. Huang, C. He, Q. Jiang, "Effects of non-darcy flow on heat-flow coupling process in complex fractured rock masses," Journal of Natural Gas Science and Engineering, vol. 83,DOI: 10.1016/j.jngse.2020.103536, 2020.

[2] C. Yao, Y. Shao, J. Yang, "Effects of fracture density, roughness, and percolation of fracture network on heat-flow coupling in hot rock masses with embedded three-dimensional fracture network," Geothermics, vol. 87,DOI: 10.1016/j.geothermics.2020.101846, 2020.

[3] Z. Ye, H. Qin, Y. Chen, Q. Fan, "An equivalent pipe network model for free surface flow in porous media," Applied Mathematical Modelling, vol. 87, pp. 389-403, DOI: 10.1016/j.apm.2020.06.017, 2020.

[4] Z. Ye, Q. Fan, S. Huang, A. Cheng, "A one-dimensional line element model for transient free surface flow in porous media," Applied Mathematics and Computation, vol. 392,DOI: 10.1016/j.amc.2020.125747, 2021.

[5] I Faoro, D Elsworth, T Candela, "Evolution of the transport properties of fractures subject to thermally and mechanically activated mineral alteration and redistribution," Geofluids, vol. 16 no. 3, pp. 396-407, DOI: 10.1111/gfl.12157, 2016.

[6] S. Huang, Z. Lu, Z. Ye, Z. Xin, "An elastoplastic model of frost deformation for the porous rock under freeze-thaw," Engineering Geology, vol. 278,DOI: 10.1016/j.enggeo.2020.105820, 2020.

[7] S. Huang, Y. Ye, X. Cui, A. Cheng, G. Liu, "Theoretical and experimental study of the frost heaving characteristics of the saturated sandstone under low temperature," Cold Regions Science and Technology, vol. 174,DOI: 10.1016/j.coldregions.2020.103036, 2020.

[8] Z. C. Tang, Q. Z. Zhang, "Elliptical Hertz-based general closure model for rock joints," Rock Mechanics and Rock Engineering, vol. 54 no. 1, pp. 477-486, DOI: 10.1007/s00603-020-02275-0, 2021.

[9] Z. C. Tang, Q. Z. Zhang, J. Peng, "Effect of thermal treatment on the basic friction angle of rock joint," Rock Mechanics and Rock Engineering, vol. 53 no. 4, pp. 1973-1990, DOI: 10.1007/s00603-019-02026-w, 2020.

[10] D. T. Snow, "A parallel plate model of fractured permeable media," 1965. M.Sc. thesis

[11] Z. Tang, Y. Jiao, "Choosing appropriate appraisal to describe peak spatial features of rock joint profiles," International Journal of Geomechanics, vol. 20 no. 4,DOI: 10.1061/(asce)gm.1943-5622.0001641, 2020.

[12] Zhi Cheng Tang, "Experimental investigation on temperature-dependent shear behavior of granite discontinuity," Rock Mechanics and Rock Engineering, vol. 53 no. 9, pp. 4043-4060, DOI: 10.1007/s00603-020-02149-5, 2020.

[13] R. W. Zimmerman, A. Al-yaarubi, C. C. Pain, "Non-linear regimes of fluid flow in rock fractures," International Journal of Rock Mechanics & Mining Sciences, vol. 41 no. 3, pp. 163-169, DOI: 10.1016/j.ijrmms.2004.03.036, 2004.

[14] Z. Y. Zhang, J. Nemcik, "Fluid flow regimes and nonlinear flow characteristics in deformable rock fractures," Journal of Hydrology, vol. 477 no. 16, pp. 139-151, DOI: 10.1016/j.jhydrol.2012.11.024, 2013.

[15] J. Q. Zhou, S. H. Hu, S. Fang, "Nonlinear flow behavior at low Reynolds numbers through rough-walled fractures subjected to normal compressive loading," International Journal of Rock Mechanics & Mining Sciences, vol. 80, pp. 202-218, DOI: 10.1016/j.ijrmms.2015.09.027, 2015.

[16] Y. F. Chen, J. Q. Zhou, S. H. Hu, R. Hu, C. B. Zhou, "Evaluation of Forchheimer equation coefficients for non–Darcy flow in deformable rough–walled fractures," Journal of Hydrology, vol. 529, pp. 993-1006, DOI: 10.1016/j.jhydrol.2015.09.021, 2015.

[17] F. Xiong, Q. H. Jiang, Z. Y. Ye, "Nonlinear flow behavior through rough-walled rock fractures: The effect of contact area," Computers & Geotechnics, vol. 102 no. OCT, pp. 179-195, DOI: 10.1016/j.compgeo.2018.06.006, 2018.

[18] M. Wang, Y. F. Chen, G. W. Ma, J. Q. Zhou, C. B. Zhou, "Influence of surface roughness on nonlinear flow behaviors in 3D self-affine rough fractures: lattice Boltzmann simulations," Advances in Water Resources, vol. 96, pp. 373-388, DOI: 10.1016/j.advwatres.2016.08.006, 2016.

[19] Q. Yin, G. W. Ma, H. W. Jing, "Hydraulic properties of 3D rough-walled fractures during shearing: An experimental study," Journal of Hydrology, vol. 555, pp. 169-184, DOI: 10.1016/j.jhydrol.2017.10.019, 2017.

[20] X. Zhang, H. Chen, C. Yao, "Seepage Characteristics of Triaxial Compression-Induced Fractured Rocks under Varying Confining Pressures," International Journal of Geomechanics, vol. 20 no. 9,DOI: 10.1061/(asce)gm.1943-5622.0001796, 2020.

[21] F. Huang, C. Yao, J. Yang, "Connectivity evaluation of fracture networks considering the correlation between trace length and aperture," Applied Mathematical Modelling, vol. 88, pp. 870-887, DOI: 10.1016/j.apm.2020.07.011, 2020.

[22] G. Zhang, Y. Tian, L. Ying-Jun, "Numerical study on the mechanism of fluid flow through single rough fractures with different JRC," Science China Physics, Mechanics & Astronomy, vol. 49 no. 01, pp. 30-39, DOI: 10.1360/sspma2018-00247, 2019.

[23] C. C. Xia, X. Qian, P. Lin, "Experimental investigation of nonlinear flow characteristics of real rock joints under different contact conditions," Journal of Hydraulic Engineering, vol. 143 no. 3,DOI: 10.1061/(asce)hy.1943-7900.0001238, 2016.

[24] Y. W. Tsang, "The effect of tortuosity on fluid flow through a single fracture," Water Resources Research, vol. 20 no. 9, pp. 1209-1215, DOI: 10.1029/wr020i009p01209, 1984.

[25] S. R. Brown, "Fluid flow through rock joints: the effect of surface roughness," Journal of Geophysical Research: Solid Earth, vol. 92 no. B2, pp. 1337-1347, DOI: 10.1029/jb092ib02p01337, 1987.

[26] T. Liang, X. Liu, S. Wang, E. Wang, Q. Li, "Study on the fractal characteristics of fracture network evolution induced by mining," Advances in civil engineering, vol. 2018,DOI: 10.1155/2018/9589364, 2018.

[27] R. C. Liu, M. He, N. Huang, "Three-dimensional double-rough-walled modeling of fluid flow through self-affine shear fractures," Journal of Rock Mechanics and Geotechnical Engineering, vol. 12 no. 1, pp. 41-49, DOI: 10.1016/j.jrmge.2019.09.002, 2020.

[28] H. H. Liu, G. S. Bodvarsson, S. Lu, "A corrected and generalized successive random additions algorithm for simulating fractional levy motions," Mathematical Geology, vol. 36 no. 3, pp. 361-378, DOI: 10.1023/b:matg.0000028442.71929.26, 2004.

[29] Z. Y. Ye, H. H. Liu, Q. H. Jiang, "Two-phase flow properties of a horizontal fracture: The effect of aperture distribution," Advances in Water Resources, vol. 76, pp. 43-54, DOI: 10.1016/j.advwatres.2014.12.001, 2015.

[30] N. Huang, R. C. Liu, Y. J. Jiang, "Numerical study of the geometrical and hydraulic characteristics of 3D self-affine rough fractures during shear," Journal of Natural Gas Science and Engineering, vol. 45, pp. 127-142, DOI: 10.1016/j.jngse.2017.05.018, 2017.

[31] G. Rong, J. Tan, H. B. Zhan, "Quantitative evaluation of fracture geometry influence on nonlinear flow in a single rock fracture," Journal of Hydrology, vol. 589,DOI: 10.1016/j.jhydrol.2020.125162, 2020.

[32] Z. Zeng, R. Grigg, "A criterion for non-darcy flow in porous media," Transport in Porous Media, vol. 63 no. 1, pp. 57-59, DOI: 10.1007/s11242-005-2720-3, 2006.

[33] L. C. Zou, L. R. Jing, V. Cvetkovic, "Shear-enhanced nonlinear flow in rough-walled rock fractures," International Journal of Rock Mechanics & Mining Sciences, vol. 97, pp. 33-45, DOI: 10.1016/j.ijrmms.2017.06.001, 2017.

[34] Z. L. Wang, L. F. Shen, Z. M. Xu, "Influence of roughness of rock fracture on seepage characteristics," Chinese Journal of Geotechnical Engineering, vol. 38 no. 7, pp. 1262-1268, 2016.

[35] J. S. Konzuk, B. H. Kueper, "Evaluation of cubic law based models describing single-phase flow through a rough-walled fracture," Water Resources Research, vol. 40 no. 2,DOI: 10.1029/2003wr002356, 2004.

[36] M. Javadi, M. Sharifzadeh, K. Shahriar, "Critical Reynolds number for nonlinear flow through rough-walled fractures: the role of shear processes," Water Resources Research, vol. 50 no. 2, pp. 1789-1804, DOI: 10.1002/2013wr014610, 2014.

[37] Y. F. Chen, S. H. Hu, R. Hu, "Estimating hydraulic conductivity of fractured rocks from high-pressure packer tests with an Izbash’s law-based empirical model," Water Resources Research, vol. 51 no. 4, pp. 2096-2118, DOI: 10.1002/2014wr016458, 2015.

[38] J. Qian, H. Zhan, S. Luo, "Experimental evidence of scale-dependent hydraulic conductivity for fully developed turbulent flow in a single fracture," Journal of Hydrology, vol. 339 no. 3-4, pp. 206-215, DOI: 10.1016/j.jhydrol.2007.03.015, 2007.