Abstract
This study investigates the effect of fracture lower surface roughness on the nonlinear flow behaviors of fluids through fractures when the aperture fields are fixed. The flow is modeled with hydraulic pressure by solving the Navier-Stokes equations based on rough fracture models with lower surface roughness varying from to . Here, JRC represents joint roughness coefficient. The results show that the proposed numerical method is valid by comparisons between numerically calculated results with theoretical values of three parallel-plate models. With the increment of hydraulic pressure drop from 10-4 to 105 Pa/m spanning ten orders of magnitude, the flow rate increases with an increasing rate. The nonlinear relationships between flow rate and hydraulic pressure drop follow Forchheimer’s law. With increasing the JRC of lower surfaces from 1 to 19, the linear Forchheimer coefficient decreases, whereas the nonlinear Forchheimer coefficient increases, both following exponential functions. However, the nonlinear Forchheimer coefficient is approximately three orders of magnitude larger than the linear Forchheimer coefficient. With the increase in Reynolds number, the normalized transmissivity changes from constant values to decreasing values, indicating that fluid flow transits from linear flow regimes to nonlinear flow regimes. The critical Reynolds number that quantifies the onset of nonlinear fluid flow ranges from 21.79 to 185.19.
1. Introduction
Hydraulic properties of rock fractures are very important for engineering practices such as enhanced oil recovery [1], CO2 sequestration [2], and geothermal energy development [3]. The permeability/transmissivity is commonly calculated/predicted based on Darcy’s law using parallel-plate models, neglecting the effects of fracture surface roughness and inertial force [4–6]. However, for fluid flow in karst systems and in high-pressure pump tests, the fluid flow enters the nonlinear flow regime and the effect of surface roughness of fractures should be considered [7–9].
The rough surface of fractures gives rise to the increase in the flow paths, decreasing the permeability/transmissivity [10–13]. Many methods have been used to characterize fracture surface roughness, such as root-mean-square of first derivative of asperity height (), Hurst exponent (), fractal dimension (), and joint roughness coefficient (JRC) [14–17]. JRC is widely used in rock mechanics and rock engineering due to its simplicity [18]. However, the determination of JRC is subjective, significantly depending on personal experiences [19]. Therefore, some empirical functions have been proposed to quantitatively determine JRC, such as the following [17, 20]: where and are the coordinates of the fracture profile along length and height directions, respectively, and is the number of sampling points along the length direction; is the interval of sampling points along the length direction and is the corresponding height variation. It is generally accepted that with the increase in fracture surface roughness, the permeability/transmissivity decreases. For example, Liu et al. [21] reported that the transmissivity of rough fractures with is 50~70% of that of smooth fractures with , when the ratio of mechanical aperture to the length of fractures varies from 0.01 to 0.10. Huang et al. [22] concluded that when , the contact area accounts for 5~27% of the total fracture plane during shear, and the permeability of rough fracture models is approximately 26~80% of that of smooth fracture models. However, the apertures are shear-induced or assigned following Gaussian distributions with the lower surface to be a smooth plane [22, 23]. In the nature, the lower surface of fractures is rough and should be taken into consideration.
In the early studies, the cubic law is used to predict the flow behaviors of fluids through rock fractures, which is suitable for parallel-plate models [24–26]. Later, to extend the model to natural fractures with rough surfaces, the cubic law is modified by investigating the relationships between hydraulic aperture and mechanical aperture [18, 27]. Although the modified cubic law can be used to characterize fluid flow through rough fractures, the model is also simplified, in which the lower surface and the upper surface are well-mated deviating from the natural fracture profiles. Recently, the fractures with different lower and upper profiles are established and fluid flow is modeled by solving the Navier-Stokes equations [28–30]. Thus, the nonlinear flow characterizations can be well understood. However, when addressing the effect of fracture surface roughness, the profiles of both lower surface and upper surface are changed [31], and the previous studies did not estimate the effect of lower surface roughness of fracture on the nonlinear flow behaviors of fluids when the aperture fields are fixed.
The present study is aimed at studying the nonlinear hydraulic properties of rough fractures with varying roughness of lower surfaces. First, three parallel-plate models are established and the Navier-Stokes equations are solved to model fluid flow. The results are compared with theoretical values, verifying the validity of proposed numerical method. Then, six rough models with the joint roughness coefficient of the lower surface varying from 0 to 19 are utilized to estimate the nonlinear hydraulic properties. Finally, the streamline distributions at different hydraulic pressure drops, the nonlinear relationships between hydraulic pressure drop and flow rate, the evolutions of Forchheimer coefficients and , the relationships between normalized transmissivity and Reynolds number, and the variations in critical Reynolds number versus lower surface roughness are systematically analyzed and discussed.
2. Theoretical Background
Assuming that the water is a kind of incompressible Newtonian fluid, the fluid flow can be governed by the Navier-Stokes equations, written as follows [24, 32, 33]: where is the flow velocity tensor, is the fluid density, is the hydraulic pressure, is the shear stress tensor, is the time, and is the body force tensor. The convective acceleration terms, , take into consideration the effect of inertial forces and are the source of nonlinear relationships between flow rate and hydraulic pressure drop. The nonlinearity of fluid flow is generally affected by the variations in local aperture and surface roughness.
The Reynolds number () that is defined as the ratio of viscous force to inertial force can be calculated according to the following [23, 34]: where is the flow rate, is the dynamic viscosity, and is the width of fractures. For 2D fractures, it is assumed that by default.
When fluid flows with a low , the inertial force is negligibly small with respect to viscous force. In such a case, the nonlinear terms, , can be deleted. Thus, Equation (2) reduces to the cubic law by neglecting the effect of fracture surface roughness, written as follows [35]: where is the hydraulic aperture, is the hydraulic pressure drop that equals to , and is the fracture length.
Equation (4) implies that is linearly proportional to , which is also linearly correlated to the cube of . Equation (4) is applicable for modeling fluid flow at a low through parallel-plate models. Due to its simplification, Equation (4) has been widely used for estimating hydraulic properties of fractured rock masses. When the is larger than a critical value, fluid flows into the nonlinear flow regime, in which is nonlinearly correlated with . In such a case, Forchheimer’s law can be employed, written as follows [36, 37]: where and are the linear coefficient and nonlinear coefficient, respectively. The linear coefficient can be expressed as follows [38]:
Equation (4) can be rewritten as follows: where is the transmissivity and is correlated to the cube of , written as follows:
Thus, Equations (4) and (5) yield the following Equations (9) and (10), representing the transmissivity in the linear and nonlinear flow regimes, respectively [39]. where is the initial transmissivity that equals to the transmissivity in the linear flow regime.
Substituting Equation (3) into Equation (10) gives the following [39]: where is a coefficient that can be expressed as follows:
It is assumed that when , the corresponding is the critical Reynolds number () that quantifies the onset of nonlinear flow of fluids [39, 40]. implies that the nonlinear flow-induced transmissivity decrease occupies 10% of the initial transmissivity. When the applied is smaller than , the fluid flow is in the linear flow regime and Equation (4) can be employed. When the applied is larger than , the fluid flow is in the nonlinear flow regime and Equation (2) should be solved. Substituting into Equation (11) yields the following:
3. Verification of the Numerical Method
To verify the validity of the numerical method, three parallel-plate models with a length of 100 mm () and mean mechanical apertures () varying from 2.52 mm to 5.51 mm are established, as shown in Figure 1. The fluid is injected into the model through the left side and flows out of the model through the right side. The upper and lower surfaces are impermeable. The geometry of the fractures is plotted in AutoCAD and exported as SAT files. The SAT files are imported into ANSYS ICEM for meshing. The quadrilateral meshes are adopted with a maximum side length of 0.1 mm. Thus, there are more than 25 layers along the aperture direction. We have also checked the effect of number of iterations and found that the calculated results are stable after 2000 iterations. So, the number of iterations is determined as 2000. The meshed geometries of the models are saved as MESH files and then imported into ANSYS FLUENT for calculation. The water is at a temperature of 25°, which has a density () of 998.2 m3/kg and a viscosity () of 0.001 Pa·s. To guarantee a linear flow, a sufficiently small hydraulic pressure drop () of 10-4 Pa/m is applied between inlet and outlet. By solving the Navier-Stokes equations, the flow rate () through the model can be calculated. For the three models, , , and , respectively. The theoretical values of are derived according to Equation (4), and the calculated , , and , respectively. The relative errors of all cases are less than 0.8%. The numerically calculated and theoretical results are presented in Figure 2, which agree well with each other verifying the validity of the proposed numerical method. Thus, the proposed numerical method is adopted for the following analysis.


4. Numerical Models and Streamline Distributions
In the present study, the distributions of local apertures () along fracture length direction are fixed as shown in Figure 3, which are extracted from cutting lines of sheared 3D rough fractures. The minimum is 0.91 mm, which is larger than 0 guaranteeing that fluid can flow through the model. The maximum is 3.64 mm. The average is 2.52 mm, which is the same as that shown in Figure 1(a). Therefore, the same parameters used for the model shown in Figure 1(a), such as the maximum side length of and number of , are adopted for the analysis. Five rough lower surfaces of fractures with are borrowed from Barton profiles proposed by Barton and Choubey [41]. For comparison, a smooth lower surface with is added, as shown in Figure 4. The height of the upper surface () is the summation of the height of the lower surface () and the height of the local aperture (), written as follows:


(a)

(b)

(c)

(d)

(e)

(f)
Note that the at different locations for the models with different rough lower surfaces is the same. This guarantees that the mean mechanical aperture is the same and the surface roughness of lower surfaces is the only variable. Thus, the effect of lower surface roughness can be investigated.
5. Nonlinear Hydraulic Properties of Fractures
5.1. Streamline Distributions
To visually observe the flow paths, a number of particles are injected at the inlet and the streamlines are recorded according to the local flow velocity tensors. When the hydraulic pressure drop is sufficiently small, i.e., , the flow rate is relatively small neglecting the effect of inertial forces and the fluid flow is in the linear regime. When the hydraulic pressure drop is large, i.e., , the fluid flow enters the nonlinear flow regime and the effect of inertial forces cannot be negligible with respect to viscous forces. In this study, and are chosen to represent the linear and nonlinear flow regimes and the corresponding streamline distributions are presented in Figures 5 and 6, respectively. In the linear flow regime, the particles smoothly flow through the void spaces formed by the tortuous lower and upper surfaces. Since the viscous force is much larger than the inertial force, no eddies are formed. In the nonlinear flow regime, the inertial force cannot be negligible with respect to the viscous force. Many eddies located at different locations with different sizes and shapes. These eddies give rise to energy losses, decreasing the transmissivity/permeability of fractures. When , the eddies exist in the place where changes significantly and in the place where does not change robustly, there are almost no eddies (i.e., the right part of Figure 6(a)). Whereas when JRC(lower) is large (i.e., = 19), the eddies are distributed within the total aperture fields, due to the influences of local aperture variations and rough surfaces of lower and upper walls. Therefore, the energy losses more significantly with a larger JRC(lower), resulting in a smaller transmissivity/permeability.

(a)

(b)

(c)

(d)

(e)

(f)

(a)

(b)

(c)

(d)

(e)

(f)
5.2. Nonlinear Relationship between Hydraulic Pressure Drop and Flow Rate
For the six models as shown in Figure 4, spanning ten orders of magnitude are applied and the macroscopic flow rate is calculated, covering the linear flow regime, weak nonlinear flow regime, and strong nonlinear flow regime. A total of 60 fluid flow simulations are performed. The relationships between and are presented in Figure 7. The curves can be described by quadratic functions with a zero intercept. With the increment of , increases with an increasing rate, following Forchheimer’s law as shown in Equation (5). As JRC(lower) increases, the curve moves leftwards, indicating that the transmissivity decreases due to the rough surfaces. Since the curves are fitted using Equation (5), the variations in coefficients and can be calculated as shown in Figure 8. With the increment of JRC(lower), decreases and increases, both following exponential functions, written as follows: where has a unit of Pa·s·m-4 and has a unit of Pa·s2·m-7.


(a) Relationships between and JRC(lower)

(b) Relationships between and JRC(lower)
Figure 8 indicates that with the increment of roughness of lower surfaces, the linear coefficient decreases due to the increase in flow lengths induced by the increased tortuous length, yet the tortuous surface will induce energy losses contributing to the increase in the nonlinear terms and increasing the nonlinear coefficient . Although the variation trends of and are different, the values of varying from to are approximately 3 orders of magnitude larger than those of .
5.3. Variations in Normalized Transmissivity and Critical Reynolds Number
The variations in versus with JRC(lower) varying from 0 to 19 are presented in Figure 9. When the is small (i.e., less than 10), the varies negligibly small approximately equaling to 1, indicating that fluid flow is in the linear regime and the cubic law is applicable. When increases from 10 to 100, decreases from values larger than 0.9 to values smaller than 0.9, meaning that fluid flow transits from linear flow regimes to nonlinear flow regimes. With the increment of JRC(lower), the corresponding to decreases. With continuously increasing (i.e., ), the continuously decreases with an increasing rate. The relationships can be well described using Equation (11), in which the coefficient can be determined. The results show that with increasing JRC(lower) from 0 to 19, the increases from 0.0006 to 0.0051, following an exponential function as shown in Equation (11), written as follows:

The values of are in the reasonable magnitudes as reported in the literature, such as and 0.00838 by Zimmerman et al. [39] and 0.00471 when the confining pressure is 0 by Yin et al. [42]. By substituting the values of into Equation (13), the can be calculated. As shown in Figure 10, as JRC(lower) increases, the decreases, following an exponential function, written as follows:

This indicates that the rougher surface of fractures gives rise to the onset of nonlinear flow at a smaller . When JRC(lower) increases from 0 to 1, decreases from 185.19 to 101.01 by a rate of 45.46%. When JRC(lower) increases from 1 to 19, decreases from 101.01 to 21.79 by a rate of 79.43%. The exhibits a significant decrease and then gradual decrease trend with increasing surface roughness of fractures, which is similar with that reported by Liu et al. [43].
6. Conclusions
The present study investigated the effect of lower surface roughness of fractures on the nonlinear hydraulic properties by solving the Navier-Stokes equations. The streamline distributions, nonlinear relationships between hydraulic pressure drop and flow rate, evolutions of normalized transmissivity and coefficient , and critical Reynolds number versus lower surface roughness are analyzed and discussed.
The results show that the proposed numerical method is valid by comparisons with theoretical result-based parallel-plate models with different apertures. At a low hydraulic pressure drop (i.e., 10-4 Pa/m), the fluid flow is in the linear flow regimes and no eddies are formed, whereas at a high hydraulic pressure drop (i.e., 105 Pa/m), the fluid flow is in the nonlinear flow regime and a number of eddies are modeled. The eddies occur at the places where apertures change significantly and near the rough surfaces. The hydraulic pressure drop has a quadratic function with flow rate, following Forchheimer’s law, when the lower surface roughness is in the range 0~19. With the increment of lower surface roughness, the linear coefficient in Forchheimer’s law decreases and the nonlinear coefficient in Forchheimer’s law increases, both following exponential functions. The values of nonlinear coefficient in Forchheimer’s law are approximately three orders of magnitude larger than the linear coefficient in Forchheimer’s law. As Reynolds number increases, the normalized transmissivity holds a constant value of 1 and then decreases with an increasing rate, indicating that fluid flow transits from a linear regime to a nonlinear regime. With the increment of lower surface roughness from 0 to 19, the coefficient describing the variations of normalized transmissivity increases from 0.0006 to 0.0051, following an exponential relationship as shown in Figure 11, which are in the similar magnitudes as reported in the literature. With increasing lower surface roughness from 0 to 19, the critical Reynolds number decreases from 185.19 to 21.79, indicating the fluid flow is easier to enter the nonlinear flow regimes for fractures having a rougher surface.

Future works will extent the current 2D models to 3D models to estimate the effect of fracture surface roughness on nonlinear fluid flow properties. Besides, the influences of aperture field and the ratio of aperture to fracture length will be taken into consideration.
Data Availability
The data can be obtained by contacting the corresponding author.
Conflicts of Interest
The authors declare that they have no conflict of interest.
Acknowledgments
This study has been partially funded by the Natural Science Foundation of China, China (Grant No. 51909142). The support is gratefully acknowledged.