Equivalent hydraulic aperture and fracture surface roughness are two significant factors affecting the fluid flow behaviors in rock fractures. To understand the role of fracture surface roughness and aperture in the fluid flow through 3D self-affine rough fractures, roughness fracture surfaces with joint roughness coefficients equal to 2.5, 7.5, 12.5, and 17.5 were established, and the Navier–Stokes equation was used to compute the fluid flow in these 3D self-affine rough fractures with a mechanical aperture increase from 0.2 mm to 0.8 mm with a gradient of 0.2 mm. The results show that when the fracture mechanical aperture is 0.2 mm, the impact of fracture surface roughness on fluid flow is considerable, while this effect decreases obviously with the increase of fracture mechanical aperture. Comparing the permeability obtained by the Navier–Stokes equation with the cube law under different hydraulic gradients, we found that their deviation increased with the increase of hydraulic gradients. This allows for the definition of a critical hydraulic gradient (), below which the permeability can be properly predicted using the cubic law for its simplicity, and above which Forchheimer’s equation should be applied, and Forchheimer’s coefficients , , and can be calculated by the prediction equations established in this study.

1. Introduction

Groundwater circulation and contaminant migration in fractured rock have attracted substantial attention in many geoscience and engineering disciplines related to nuclear waste disposal, geothermal resource development, oil and gas extraction, CO2 geological storage, etc. [13]. Compared with rock fracture, the permeability of tight matrix is negligibly due to the permeability of a fracture conventionally being several higher magnitudes than that of rock matrix, and the fracture controls the predominant pathways of fluid flow and solute transport through fractured rock masses [46]. Individual fractures are the basic element that constitutes the fracture network, and quantitative estimation of the permeability of single fractures is the foundation for understanding the fluid flow through complex natural fractured rock masses [7, 8]. Natural fracture surfaces are typically rough with 3D self-affine distribution and display complex hydraulic behavior coming from fracture surface roughness [9]. Fluid flow through individual fractures is affected by multiple factors, such as fracture aperture, surface roughness, and inertial effects, among which aperture and surface roughness may be the two most critical parameters that significantly affect fluid flow behavior [10, 11].

The fracture is typically represented by two smooth parallel plates separated by a small distance, and fluid flow through a single fracture is assumed to be analogous to laminar flow [1214]. This derivative of the so-called “cubic law,” in which the flow flux is proportional to the cubic fracture aperture and the pressure head reduces due to surface morphology, is neglected [1517]. The parallel plate model is a qualitative description of fluid flow through real fractures. In fact, the rough fracture surface due to the spotty asperity height makes the fluid flow difficult to conform to the laminar regime [18, 19]. The fluid flow in the real fracture is tortuous, which deviates from the expectations of the cube law, and the preferred fluid flow channel has been observed in field and laboratory investigations, in which more than 90% of the total flow only occupies 5%–20% of the fracture plane [20, 21]. In addition, the mechanical aperture of a fracture is usually heterogeneous due to the irregular distribution of aperture asperities [22, 23]. Alternatively, the hydraulic aperture is proposed to describe the hydraulic characteristic of fractures with rough surfaces. Fractures have different mechanical apertures caused by surface morphology, and they may have the same hydraulic aperture [24, 25]. Many efforts have been made to establish the correlations between hydraulic aperture and mechanical aperture, in which the tortuosity factor, contact ratio, the standard deviation of mechanical apertures, and the most widely applied parameter called joint roughness coefficients are proposed [26].

To improve the understanding of the fluid flow behavior in rock fractures, quantitative and systematic analyses of the fluid flow in rock fractures with different fracture apertures and surface morphologies is necessary. The fluid flow behaviors of single fractures have been studied extensively using theoretical modeling, laboratory experiments, and numerical methods [27, 28]. Numerical simulation is of great advantage in handling various complex conditions, such as different fracture apertures with the same surface roughness, which is technically challenging in the laboratory [29]. Considering the aperture heterogeneity in individual fractures, there exists a humongous number of meshes after the model geometry is discretized, and catching the fluid flow behavior in rock fractures through computation of the Navier–Stokes (NS) equations that are composed of a set of coupled nonlinear partial derivatives is burdensome [30]. Some simplified form, cubic law as an example, is presented for computational convenience, but this equation can only be characteristic laminar flow through fractures with smooth surfaces [31]. Therefore, a criterion is needed to judge whether the simplified form should be applied.

In this study, a rough fracture surface was generated using a modified successive random addition (SRA) algorithm, and numerical simulations of fluid flow were performed on a series of 3D self-affine rough fractures with different fracture apertures and JRCs. Then, the effects of surface roughness and mechanical aperture on the evolution of permeability in different hydraulic gradients were analyzed. Finally, multivariable regressions were used to establish the mathematical expressions of critical hydraulic gradient that can judge the onset of nonlinear fluid flow, and the prediction method of coefficients ( and ) involved in Forchheimer’s law was constructed simultaneously.

2. Model Generation

2.1. The Method for Rough Fracture Surface Generation

The surface morphology of the structural surface in the rock mass is rough and conforms to the fractal distribution, so the fractal algorithm can be used to generate the rough surface of the fracture [32]. The successive random additions algorithm (SRA), which has been widely used in previous studies due to its high efficiency, was also used in this study [33]. In two-dimensional SRA, a single-valued continuous function is used to describe the undulations of rough fracture surfaces, which is defined as and obeys the Gaussian normal distribution with mean zero and variance at distance . The self-affine surface is defined as: where is the mathematical expectation, is a constant, is the Hurst exponent, and is the variance calculated by:

Using the SRA method, rough surfaces with s of 2.5, 7.5, 12.5, and 17.5 were generated and are shown in Figure 1.

2.2. Stress Response of Rock Fractures

The effects of stress/deformation on the deformation of fractured rock are one of the major concerns for performing fluid flow through fractures. The influence of stress on fractured rock has been widely investigated, and the results indicate that the hydraulic conductivity and flow patterns will change under different stress conditions [34]. Figure 2 shows that when stress is applied to fractured rock, the fracture aperture will increase or decrease.

The change in fracture aperture caused by stress can be calculated by the following formula: where is the initial mechanical pore size and and are the changes in mechanical aperture due to normal and shear stress, respectively. A widely used nonlinear model between fracture aperture change and normal stress is the hyperbolic relationship [35]: where is the maximum possible closure of the fracture and is the normal stiffness.

2.3. Performed Fluid Flow in Rock Fracture

COMSOL Multiphysics, a multiphysics coupling simulation software based on finite elements, is widely used in the calculation of seepage in rock fractures. A widely used method is employed in this study to construct the rough surface of a fracture. The rough surface generated by the SRA algorithm mentioned was first imported into COMSOL and then raised a small distance to establish a fracture cavity model [35]. It is notable that the stress and displacement of the fracture surface are not considered because the stress is not the factor that this study focused on. By raising the fracture surface with different distances, fracture apertures of 0.2 mm, 0.4 mm, 0.6 mm. and 0.8 mm are constructed. The role of fracture morphology is investigated by importing surfaces with different roughness surface generation in Section 2.1 based on fractal theory [36]. In total, the fracture surfaces of four different s were imported to analyze the effect of fracture surface roughness on fluid flow, and four distances of fracture (in Figure 3) were set to investigate the role of mechanical aperture in fluid flow. In total, sixteen models were studied.

To perform the fluid flow simulation in the geometric model, a velocity boundary is set on the left side of the model, the outlet on the right boundary is set to a constant pressure, and the other four faces are set as impervious boundaries. To study the different flow states in the fracture, the inlet velocity was set as follows: low velocity flow, with a gradient of , increasing from to ; medium-speed flow, with a gradient of , increasing from to ; high-speed flow, with as the gradient, increasing from to .

3. Theoretical Background

3.1. Governing Equations of Fluid Flow

Numerical simulation calculations are performed in the COMSOL numerical software based on the finite element method by directly solving the Navier–Stokes (N–S) equations. Incompressible fluid flow in a fracture is governed by [37, 38]: where is the fluid density, is the velocity vector, is the hydraulic pressure, is the viscous coefficient, and is the body force.

To overcome the difficulty in solving the N–S equation caused by the nonlinear term, a simplified model is proposed in which the fracture is regarded as two smooth parallel plates separated by a small distance in which the fluid flow at low velocity follows the cubic law [39]: where is the flow rate; is the width of the model; is the hydraulic aperture, which is usually obtained by measurement in the test and/or field; is the fluid viscosity; is the hydraulic pressure difference between the inlet and outlet. When the fluid flow velocity increases, the fluid velocity and hydraulic gradient no longer obey the linear law, and the cubic law no longer applies. The most widely used mathematical description of nonlinear flow relationships in fractured and porous media is Forchheimer’s law [40, 41], which is a zero-intercept quadratic equation: where and are two coefficients. In this study, the hydraulic gradient , a dimensionless parameter, is applied to quantify the hydraulic pressure difference between the inlet and outlet, which is defined by . Then, Equation (7) can be written as:

To quantitatively determine the nonlinear effect in fluid flow, parameter is proposed, which is defined as [42, 43]:

Substituting Equations (7) and (8) into Equation (9) yields:

Conventionally, the critical value of is defined as 0.1, which means that when the pressure drop caused by nonlinear fluid flow accounts for more than 10% of the total pressure drop, the nonlinear fluid flow should be considered in the computation [44, 45].

3.2. Roughness Surface Characterization

The joint roughness coefficient (), a widely used parameter usually in the range of 0 to 20, was proposed by Barton and Choubey [46] for quantifying the roughness of a curve, which was calculated by [47]: where is the root mean square of the first derivative of the profile and can be obtained by: where and are the coordinates of the -th sampling point, and is the number of sampling points.

For the three-dimensional rough fracture surface, this study selects 20 curves along the direction of fluid flow at equal intervals as the samples, and the average value of their JRC was calculated to characterize the roughness of the fracture surface.

4. Fluid Flow in Discrete Fracture Networks

4.1. The Relationship between the Hydraulic Gradient and Flow Flux under Different JRCs and Mechanical Apertures ()

In the process of numerical simulation, flow fluid behavior in different regimes is realized by changing the inflow velocity. Solving the N–S equation is time-consuming. Taking equal to 17.5 as an example, when the fracture aperture was set to 0.2 mm, 0.4 mm, 0.6 mm, and 0.8 mm, the number of meshes with an average quality greater than 0.7 was 5 119 468, 6 657 317, 7 350 285, and 8 645 096, and the corresponding computation times were 42 h, 56 h, 78 h, and 91 h, respectively. Figure 4 shows that the pressure gradient increases nonlinearly with the increase of flow rate for different fracture apertures with the same . The curvature of the curve in Figure 4 gradually increases, which indicates that the nonlinearity of fluid flow is increasing, and the nonlinear relationship between the pressure gradient and flow rate is becoming increasingly obvious. The permeability of the fracture was calculated based on the cubic law simultaneously; that is, the pressure loss represented by the nonlinear term in the Forchheimer equation is ignored. The results show that there is a nonnegligible deviation between the results obtained by the cube law and the N–S equation. In the fracture with an aperture of 0.2 mm, when the flow rate is , the hydraulic gradient obtained by the N–S equation is 40, and the value obtained by the cube law is 100, which is approximately 2.5 times that of the former. This indicates that when the flow velocity (pressure) is above a certain value, the permeability of the fracture cannot be obtained by the cubic law. As the mechanical aperture increases from 0.2 mm to 0.8 mm, the quadratic term coefficients of the fitting equation decrease, which means that the nonlinear effect is weakened. Further study of the relationship between coefficients and and the mechanical aperture will be discussed in a later section.

4.2. The Relationship between the Hydraulic Gradient and Flow Flux under Different Fracture Surface Roughness

Figure 5 shows the relationship between the flow rate and hydraulic gradient at different s with a mechanical aperture of 0.2 mm. With the increase in fracture surface roughness, the pressure gradient increases to achieve the same flow flux; for example, when the flow flux is , the pressure gradient is 55, 65, 80, and 100. This means that the greater the roughness of the fracture surface is, the greater the resistance of the fluid flow through the fracture. The increase in the quadratic term coefficient of the fitted curves also indicates that the increase in fracture surface roughness leads to more obvious nonlinear effects of the fluid in the fracture. Comparing the hydraulic gradients calculated with the NS equation and the cube law when achieving the same flow rate, the deviation between them increases as the fracture surface roughness increases or the flow flux increases. That is, the cube law can only be used at low flow rates and a relatively smooth fracture surface. The applicable conditions of the cube law and Forchheimer’s law will be discussed quantitatively in later sections.

4.3. The Relationship between Permeability and Hydraulic Gradient under Different and JRC

The equivalent permeability of a rock fracture can be obtained by back-calculating by Darcy’s law: where is the area of the outlet in the three-dimensional model and represents the length of the outlet boundary in the two-dimensional model.

By calculating the permeability using Equation (13), the evolution of the permeability of rock fractures with different fracture surface roughness and different mechanical apertures in terms of the hydraulic gradient is shown in Figure 6. In Figure 6, when the fluid flow velocity is small, the fluid flow state is laminar flow, which conforms to Darcy’s law, and the permeability coefficient remains constant. When the fluid flow velocity gradually increases, nonlinear flow appears, and the flow regime changes into weak nonlinearity. As the fluid flow velocity continues to increase, the nonlinear flow in the fluid intensifies, the flow transforms to a strongly nonlinear state, the hydraulic gradient and flow rate conform to Forchheimer’s law, and the permeability decreases sharply with the increase of .

To analyze the fluid flow regime in rough fractures, Figure 7 plots the flow velocity distribution diagrams under four mechanical apertures. The fluid flow in different rock fractures has the same pattern when flowing at a low speed, which belongs to Darcy flow without forming a vortex. As the flow rate increases, vortices begin to appear, energy loss occurs in the flow process, and the equivalent permeability begins to decrease. The flow velocity continues to increase, a large number of eddies appear, the flow in the fracture cavity is more complex, and the nonlinear effect is more intense. This analysis is consistent with the variation trend of the model permeability in Figure 6.

4.4. The Deviation of Permeability Calculated by the N–S Equation and Cubic Law

To quantify the impact of nonlinear effects in fluid flow on model permeability, a parameter that is characteristic of the deviation of permeability calculated by the N–S equation and cubic law is proposed, which is defined as: where is the permeability of the rock fracture in laminar flow and is the permeability of the rock fracture at an arbitrary hydraulic gradient.

Figure 8 shows the changing trend of permeability deviation in terms of the variety of hydraulic gradient , keeps increasing as increases. Flow flux and hydraulic gradient have a linear relationship when a small hydraulic gradient (below 0.1) is applied to the fracture model, fluid flow conforms to the Darcy flow, and permeability remains constant, although increases. As increases (approximately 0.1 to 0.5), the appearance of eddies results in a nonlinear relationship between the flow flux and hydraulic gradients, and fluid flow transforms into a weakly nonlinear regime. As the hydraulic gradient continues to increase, the nonlinearity of fluid flow becomes more pronounced and shows a stronger nonlinear flow phenomenon, and the permeability deviation increases sharply. In addition, under the same hydraulic gradient and surface roughness, the permeability deviation also increases as the mechanical aperture increases. This means that the mechanical aperture and fracture surface roughness are worth considering in the investigation of the effect of hydraulic gradient on the fluid flow regime and permeability of rock fractures.

4.5. A Prediction Model for Critical Hydraulic Gradient and Forchheimer’s Coefficients and

The Forchheimer coefficients and obtained by fitting with different mechanical apertures and surface roughness are listed in Table 1. As the mechanical aperture of the fracture increases, both and shown in Table 1 decrease, which is consistent with the negative index of the aperture in Equation (15). Although the value of is 5 orders of magnitude larger than the value of , for the increase in hydraulic aperture from 0.2 mm to 0.8 mm, the changes in and are both 2 orders of magnitude. The only factor affecting coefficient is the equivalent fracture aperture. The factors affecting coefficient are related to the spatial distribution of the surface roughness and other factors. Therefore, the change in coefficient is more obvious than that of . In addition, the critical hydraulic gradient obtained in the last section is summarized in Table 1. To establish a model for evaluating the Forchheimer coefficients and and the critical hydraulic gradient , a multiple regression algorithm was used to fit the simulated data, and the best fit formula is:

The prediction values of , , and calculated by Equation (15) are listed in Table 1, and their correlation coefficients are larger than 0.99, indicating that the prediction model can reasonably forecast the Forchheimer coefficients and , as well as the critical hydraulic gradient of rock fractures.

5. Conclusions

The primary goal of this study is to investigate the fluid flow characteristics in rock fractures and the model for judging the critical condition of nonlinear flow. The research carried out a numerical study on the fluid flow characteristics in the rock fracture by establishing geometric models of different mechanical apertures and rough fracture surfaces. First, the fracture surface of different roughness is generated based on SRA algorithms and geometric models with different apertures are constructed in COMSOL. Next, fluid dynamics calculations were performed by directly solving the N–S equation, and the permeability calculated by the cubic law was considered simultaneous for comparison. The results show that when the hydraulic gradient is small, there is a linear relationship between the flow and the hydraulic gradient. When the hydraulic gradient increases, there is a nonlinear relationship between the flow rate and the hydraulic gradient, at which point the equivalent permeability decreases, resulting in a significant deviation between using the N–S equation and the cube law. When is small, the flow velocity in models with different mechanical apertures has the same distribution form, and no vortex exists and conforms to the Darcy flow. When is larger, the flow velocity increases, and the eddy current is generated and causes energy loss in the flow process, which in turn reduces the flow rate and the equivalent permeability, which transforms into strong nonlinear flow. The quantitative analysis of the factors affecting the flow state in fractured rock masses was carried out first. Then, the critical hydraulic gradient for the condition that nonlinear flow cannot be ignored was defined, and the prediction model of , , and involving the fracture aperture and surface roughness was obtained using the multiple regression algorithm. With the increase in the mechanical aperture, the area generated by the vortex increases, and the critical hydraulic gradient and the coefficients in the Forchheimer equation all show a power-exponential function decrease. When the hydraulic gradient is less than , the cubic law can be used for its simplicity and convenience. When the hydraulic gradient is greater than , the Forchheimer equation should be applied, and the coefficient can be calculated according to the empirical formula proposed in this study.

Although this research quantitatively studies the nonlinear flow in rock fractures, the effect of model size on fluid flow characteristics is not considered, and the relationship between stress, fracture aperture, and fluid flow under loading conditions is not realized. In future work, we will conduct further research on the above defects.

Data Availability

The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.


This study was financially supported by the National Natural Science Foundation of China (Grant No. 51674047 and 51911530152).