Effect of Lower Surface Roughness on Nonlinear Hydraulic Properties of Fractures

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 
 
 drop
 =
 
 
 10
 
 
 −
 4
 
 
 ~
 
 
 10
 
 
 5
 
 
  
 Pa
 /
 m
 
 by solving the Navier-Stokes equations based on rough fracture models with lower surface roughness varying from 
 
 JRC
 =
 1
 
 to 
 
 JRC
 =
 19
 
 . 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.


Introduction
Hydraulic properties of rock fractures are very important for engineering practices such as enhanced oil recovery [1], CO 2 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][5][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][8][9].
The rough surface of fractures gives rise to the increase in the flow paths, decreasing the permeability/transmissivity [10][11][12][13]. Many methods have been used to characterize fracture surface roughness, such as root-mean-square of first derivative of asperity height (Z 2 ), Hurst exponent (H), fractal dimension (D f ), and joint roughness coefficient (JRC) [14][15][16][17]. JRC is widely used in rock mechanics and rock engineer-ing 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]: JRC = 32:2 + 32:47 log Z 2 , where x i and z i are the coordinates of the fracture profile along length and height directions, respectively, and M is the number of sampling points along the length direction; x i − x i−1 is the interval of sampling points along the length direction and z i − z i−1 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 Z 2 = 0:5 is 50~70% of that of smooth fractures with Z 2 = 0, 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 JRC = 6~18, 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][25][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][29][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 a and b, the relationships between normalized transmissivity and Reynolds number, and the variations in critical Reynolds number versus lower surface roughness are systematically analyzed and discussed.

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 u is the flow velocity tensor, ρ is the fluid density, p is the hydraulic pressure, T is the shear stress tensor, t is the time, and f is the body force tensor. The convective acceleration terms, ðu ⋅ ∇Þu, 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 (Re) that is defined as the ratio of viscous force to inertial force can be calculated according to the following [23,34]: where Q is the flow rate, μ is the dynamic viscosity, and w is the width of fractures. For 2D fractures, it is assumed that w = 1 by default.
When fluid flows with a low Re, the inertial force is negligibly small with respect to viscous force. In such a case, the nonlinear terms, ðu ⋅ ∇Þu, can be deleted. Thus, Equation (2) reduces to the cubic law by neglecting the effect of fracture surface roughness, written as follows [35]: where e is the hydraulic aperture, ∇P is the hydraulic pressure drop that equals to ΔP/ΔL, and L is the fracture length. Equation (4) implies that Q is linearly proportional to ∇P, which is also linearly correlated to the cube of e. Equation (4) is applicable for modeling fluid flow at a low Re through parallel-plate models. Due to its simplification, Equation (4) has been widely used for estimating hydraulic properties of fractured rock masses. When the Re is larger than a critical value, fluid flows into the nonlinear flow regime, in which Q is nonlinearly correlated with ∇P. In such a case, Forchheimer's law can be employed, written as follows [36,37]: where a and b are the linear coefficient and nonlinear coefficient, respectively. The linear coefficient a can be expressed as follows [38]: Equation (4) can be rewritten as follows: where T is the transmissivity and is correlated to the cube of e, 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 T 0 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 T/T 0 = 0:9, the corresponding Re is the critical Reynolds number (Re c ) that quantifies the onset of nonlinear flow of fluids [39,40]. T/T 0 = 0:9 implies that the nonlinear flow-induced transmissivity decrease occupies 10% of the initial transmissivity. When the applied Re is smaller than Re c , the fluid flow is in the linear flow regime and Equation (4) can be employed. When the applied Re is larger than Re c , the fluid flow is in the nonlinear flow regime and Equation (2) should be solved. Substituting T/T 0 = 0:9 into Equation (11) yields the following:

Verification of the Numerical Method
To verify the validity of the numerical method, three parallelplate models with a length of 100 mm (L = 100 mm) and mean mechanical apertures (E) varying from 2.52 mm to 5.51 mm are established, as shown 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.

Numerical Models and Streamline Distributions
In the present study, the distributions of local apertures (E local ) along fracture length direction are fixed as shown in Figure 3, which are extracted from cutting lines of sheared 3D rough fractures. The minimum E local is 0.91 mm, which is larger than 0 guaranteeing that fluid can flow through the model. The maximum E local is 3.64 mm. The average E local 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 meshes = 0:1 mm and number of iterations = 2000, are adopted for the analysis. Five rough lower surfaces of fractures with JRC ðlowerÞ = 1, 3, 9, 11, and 19 are borrowed from Barton profiles proposed by Barton and Choubey [41]. For comparison, a smooth lower surface with JRCðlowerÞ = 0 is added, as shown in Figure 4. The height of the upper surface (h upper ) is the summation of the height of the lower surface (h lower ) and the height of the local aperture (E local ), written as follows: Note that the E local 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.     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 JRCðlowerÞ = 0, the eddies exist in the place where E local changes significantly and in the place where E local 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.

Nonlinear Relationship between Hydraulic Pressure Drop and Flow
Rate. For the six models as shown in Figure 4, −∇ P = 10 −4~1 0 5 Pa/m spanning ten orders of magnitude are applied and the macroscopic flow rate Q 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 Q and −∇P are presented in Figure 7. The Q~−∇P curves can be described by quadratic functions with a zero intercept. With the increment of Q, −∇P 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 a and b can be calculated as shown in Figure 8. With the increment of JRC(lower), a decreases and b increases, both following exponential functions, written as follows: where a has a unit of Pa·s·m -4 and b has a unit of Pa·s 2 ·m -7 . 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

Geofluids
contributing to the increase in the nonlinear terms and increasing the nonlinear coefficient b. Although the variation trends of a and b are different, the values of b varying from 2:15 × 10 6 Pa · s · m −4 to 4:72 × 10 6 Pa · s 2 · m −7 are approximately 3 orders of magnitude larger than those of a.

Variations in Normalized Transmissivity and Critical
Reynolds Number. The variations in T/T 0 versus Re with JRC(lower) varying from 0 to 19 are presented in Figure 9. When the Re is small (i.e., less than 10), the T/T 0 varies negligibly small approximately equaling to 1, indicating that fluid flow is in the linear regime and the cubic law is applicable. When Re increases from 10 to 100, T/T 0 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 Re corresponding to T/T 0 = 0:9 decreases. With continuously increasing Re (i.e., Re > 100), the T/T 0 continuously decreases with an increasing rate. The T/T 0~R e 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) The values of β are in the reasonable magnitudes as reported in the literature, such as β = 0:00477 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 Re c can be calculated. As shown in Figure 10, as JRC(lower) increases, the Re c decreases, following an exponential function, written as follows:

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 parallelplate 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., 10 5 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.