The Effect of Fracture Surface Roughness and Mechanical Aperture on the Onset of Nonlinear Flow Behaviors in 3D Self-Affine Rough Fractures

,


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, CO 2 geological storage, etc. [1][2][3]. 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 [4][5][6]. 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 [12][13][14]. 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 [15][16][17]. 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 Jc that can judge the onset of nonlinear fluid flow, and the prediction method of coefficients (A and B) involved in Forchheimer's law was constructed simultaneously.

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 aðx i , y j Þ is used to describe the undulations of rough fracture surfaces, which is defined as ½aðx + Δx, y + ΔyÞ − aðx, yÞ and obeys the Gaussian normal distribution with mean zero and variance σ 2 at distance lðl = ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi ðΔxÞ 2 + ðΔyÞ 2 q Þ. The self-affine surface is defined as: where h·i is the mathematical expectation, r is a constant, H is the Hurst exponent, and σ 2 is the variance calculated by: Using the SRA method, rough surfaces with JRCs of 2.5, 7.5, 12.5, and 17.5 were generated and are shown in Figure 1.

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 b 0 is the initial mechanical pore size and Δb n and Δb s are the changes in mechanical aperture due to normal and shear stress, respectively. A widely used nonlinear model between fracture aperture change Δb n and normal stress σ n is the hyperbolic relationship [35]: where Δb m is the maximum possible closure of the fracture and k n is the normal stiffness.

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 2 Geofluids 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 JRCs 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 4 × 10 −4 m/s, increasing from 4 × 10 −4 m/s to 16 × 10 −4 m/s; medium-speed flow, with a gradient of 2 × 10 −3 m/s, increasing from 2 × 10 −4 m/s to 18 × 10 −4 m/s; high-speed flow, with 2 × 10 −2 m/s as the gradient, increasing from 2 × 10 −2 m/s to 100 × 10 −2 m/s.

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, U is the velocity vector, P is the hydraulic pressure, μ is the viscous coefficient, and f 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 Q is the flow rate; w is the width of the model; e is the hydraulic aperture, which is usually obtained by 3 Geofluids measurement in the test and/or field; μ is the fluid viscosity; ΔP 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 A ′ and B ′ are two coefficients. In this study, the hydraulic gradient J, a dimensionless parameter, is applied to quantify the hydraulic pressure difference between the inlet and outlet, which is defined by J = ∇P/ðρgÞ. Then, Equation (7) can be written as: To quantitatively determine the nonlinear effect in fluid flow, parameter E is proposed, which is defined as [42,43]: Substituting Equations (7) and (8) into Equation (9) yields: Conventionally, the critical value of E 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].

Roughness Surface Characterization.
The joint roughness coefficient (JRC), 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 Z 2 is the root mean square of the first derivative of the profile and can be obtained by: where x i and δ i are the coordinates of the i-th sampling point, and N t 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.  Figure 4 shows that the pressure gradient increases nonlinearly with the increase of flow rate for different fracture apertures with the same JRC. 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  7 Geofluids 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 per-meability 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 B of the fitting equation decrease, which means that the nonlinear effect is   Geofluids weakened. Further study of the relationship between coefficients A and B and the mechanical aperture e m will be discussed in a later section.

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 JRCs 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 20 × 10 −7 m 3 /s , 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. 9 Geofluids permeability of a rock fracture can be obtained by backcalculating by Darcy's law:

The Relationship between Permeability and Hydraulic
where S 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 K decreases sharply with the increase of J.
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.

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 k 0 is the permeability of the rock fracture in laminar flow and k 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 J, δ keeps increasing as J 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 J increases. As J 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.   Table 1. As the mechanical aperture of the fracture increases, both A and B shown in Table 1 decrease, which is consistent with the negative index of the aperture in Equation (15). Although the value of B is 5 orders of magnitude larger than the value of A, for the increase in hydraulic aperture from 0.2 mm to 0.8 mm, the changes in A and B are both 2 orders of magnitude. The only factor affecting coefficient A is the equivalent fracture aperture. The factors affecting coefficient B are related to the spatial distribution of the surface roughness and other factors. Therefore, the change in coefficient B is more obvious than that of A. In addition, the critical hydraulic gradient Jc obtained in the last section is summa-rized in Table 1. To establish a model for evaluating the Forchheimer coefficients A and B and the critical hydraulic gradient Jc, a multiple regression algorithm was used to fit the simulated data, and the best fit formula is: The prediction values of A, B, and Jc 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

Geofluids
A and B, as well as the critical hydraulic gradient Jc of rock fractures.

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 J 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 J 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 Jc for the condition that nonlinear flow cannot be ignored was defined, and the prediction model of A, B, and Jc 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 powerexponential function decrease. When the hydraulic gradient is less than Jc, the cubic law can be used for its simplicity and convenience. When the hydraulic gradient is greater than Jc, 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.