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

#### 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 *δ*^{2}. The self-affinity relating to fBm obeys the following expressions: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, and 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 (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 should be added into all points, in which(3)Step (2) is repeated for each newly generated square, and the random values from should be added to all heights as well until squares are created at iteration(4)No new square is inserted again, but we keep adding random values from to all points, in whichwhere *j* = *n* + 1, *n* + 2, …, NM; NM is large enough, so *δ*_{NM}/*δ*_{0} is negligible

**(a)**

**(b)**

**(c)**

The upper and lower surfaces of rough fracture are generated using the SRAM. The height of low surface is *Z*_{1}(*x*,*y*), and the height of the upper surface can be calculated by the following formula:where *u* is the mean aperture between the upper and lower surfaces. The distribution function of the fracture aperture can be expressed as

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.

**(a)**

**(b)**

**(c)**

**(d)**

#### 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]:where *Q* is the volumetric flow rate per unit time, is the hydraulic aperture of the fracture, is the fracture width, *μ* is the dynamic viscosity, and 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: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^{2}) 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:where 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:

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

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:where **u** and are the velocity vector and Hamiltonian operator, respectively. In this study, the fluid density is 997.1 (kg/m^{3}), and the dynamic viscosity is 0.894 ×10^{−3} (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.

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^{−8}∼3.228 × 10^{−6}·m^{3}/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.

#### 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*^{2} 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.

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

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 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^{−7}·m^{3}/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.

**(a)**

**(b)**

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

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^{−6}·m^{3}/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].

**(a)**

**(b)**

**(c)**

##### 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 , 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: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*^{2} 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:

It should be noted that when the value of 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.

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.

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

#### Data Availability

The data used in the present study can be made available upon request from the authors.

#### Conflicts of Interest

The authors declare that they have no conflicts of interest.

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