Investigation on Nonlinear Flow Behavior through Rock Rough Fractures Based on Experiments and Proposed 3-Dimensional Numerical Simulation

Rock fractures as the main flow channels, their morphological features, and spatial characteristics deeply influence the seepage behavior. Reservoir sandstones as a case study, four splitting groups of fractures with different roughness were scanned to get the geometric features, and then the seepage experiments were taken to analyze the relationship of the pressure gradient ∇Pand flow rate Q, and the critical Reynolds number(Rec) and wall friction factor (f ) were determined to explain the translation of linear seepage to nonlinear seepage condition. Based on the scanning cloud data of different rough fractures, the fractures were reconstructed and introduced into the COMSOL Multiphysics software; a 3-dimensional seepage model for rough fractures was calibrated and simulated the seepage process and corresponding pressure distribution, and explained the asymmetry of flow velocity. And also, the seepage characteristics were researched considering aperture variation of different sample fractures; the results indicated that increasing aperture for same fracture decreased the relative roughness, the fitting coefficients by Forchheimer formula based on the data ∇P~Q decreased, and the figures about the coefficients and corresponding aperture described nonlinear condition of the above rough fractures. In addition, the expression of wall friction factor was derived, and relationship of f , Re, and relative roughness indicated that f increased with increasing fracture roughness considering the same aperture, resulting in nonlinear flow more easily, otherwise is not, showing that f could be used to describe the seepage condition and corresponding turning point. Finally, it can be seen from the numerical results that the nonlinearity of fluid flow is mainly caused by the formation of eddies at fracture intersections and the critical pressure gradient decreases with increasing angle. And also, analysis about the coefficient B in the Forchheimer law corresponding to fracture intersections considering the intersecting angle and surface roughness is proposed to reveal the flow nonlinearity. The above investigations give the theoretical support to understand and reveal the seepage mechanism of the rock rough fractures.


Introduction
Rock fractures are common in the field of the rock engineering. And the fractures as the controlled seepage channels play an important role on the groundwater assessment, oil and gas development, and mass transportation; how to solve the flow behavior through the complex fractures is key to above engineering [1].
And many researches indicate that natural fractures with complex morphology form tortuous seepage channels to influence the seepage behavior. Some researchers such as Babadagli et al. [2] and Singh et al. [3] have researched the cubic law to describe the low-velocity flow in the fractures considering the fracture morphology. Considering the flow velocity increasing, Yuan et al. [4] and Zimmerman et al. [5] have conducted flow tests and explained the flow nonlinearity, indicating that greater flow velocity in the rough fractures and wall friction coefficient enlarge the inertial forces to show nonlinear flow behavior; therefore, the cubic law considering the viscous forces cannot be applicable for the flow through the rough fractures, especially under higher flowing conditions.
In recent years, many researchers have conducted nonlinear flow tests, theories, and numerical simulations; Ranjith and Darlington [6], Qian et al. [7], Zhang and Nemcik [8], Zoorabadi et al. [9] Chen and Zhou [10], and Xiong et al. [11] have conducted flow tests considering different fracture roughness of a large amount of rocks to analyze the influence of the fracture roughness on the flow nonlinearity under loading and unloading conditions, revealing the seepage mechanism of the nonlinear flow behavior through the rough fractures. The testing achievements mentioned above give supports to set up a rational model to describe the flow nonlinearity. Therefore, some researchers have proposed a few theoretical models such as seepage cavity model by Weimin et al. [12], improved cubic model by Hongguang et al. [13], T model similar to the Forchheimer model by Javadi et al. [14], and a model based on corrected Reynolds equation by Lee et al. [15] to illuminate the seepage characteristics. However, the theoretical models with many simplifications cannot completely describe the nonlinear flow behavior through the rough fractures. Therefore, numerical simulations based on the experiments and theories have been developed for complex flow behavior through the rough fractures.
Wang et al. [16] and Zou et al. [17] have proposed some numerical models based on Boltzmann and Fluent software to analyze the flow nonlinearity in the rough fractures. Furthermore, the above tests, theoretical analysis, and numerical simulations have given many useful achievements for deep investigation of the flow behavior through the complex fractures. However, rougher fracture surfaces may influence on the friction loss, inertial forces, and seepage behavior, while the fluid flowing through rough fractures. Nazridoust et al. [18], Yang et al. [19], Zhang et al. [20], Zhou et al. [21], and Xiaopeng et al. [22] have conducted experiments and theories to propose the seepage model considering wall friction effect and reveal the evolution mechanism of the fluid flow in rough fractures considering the wall friction effect, which deeply explains that the seepage mechanism of the flow through rough fractures especially represents the transformation of the seepage behavior from the linearity to the nonlinearity while the fracture roughness and flow velocity change suddenly. However, above investigations are focused on flow behavior through the 2-dimensional fractures; advanced research on the flow characteristics in 3dimensional fractures should be deeply conducted.
In addition, when fluid is injected into a fracture intersection, the flow pressure and flow rate distribution are affected by many factors such as fluid properties [23], fracture intersection angle [24,25], fracture surface roughness [23], and pressure gradient [26]. The fluid flow behavior derived for a single fracture [27][28][29] is not applicable to the case of fracture   3 Geofluids intersections, especially when the nonlinear flow regime is considered. Numerical simulations have been performed by Liu et al. [30] and Xiong et al. [31] to estimate the effects of fracture intersections on the nonlinear flow behaviors in crossing fractures. It is found that a nonlinear flow regime occurs in which the inertial effects caused by the fracture intersections should not be neglected. Although many studies described have been conducted experimentally and numeri-cally, the influence mechanisms of fracture intersection angle and roughness on fluid flow behavior are rarely discussed. In addition, the corresponding nonlinear flow model is also rarely reported.
Although many studies described were conducted experimentally, theoretically, and numerically, the flow through 3dimensional rough fractures considering the space distribution has not been deeply investigated. Taking a reservoir    Figure 1(a) that the sandstone is gray and white, and the size of all tested sandstone specimens is cylindrical with 50 mm in diameter and 100 mm in length approximately according to the method by the International Society for Rock Mechanics (ISRM). And then, the Brazilian test plotted in (Figure 1(b)) is taken to obtain the artificial rough fractures. As for this testing process, the prepared samples are placed in the special clamp and then loaded with different loading rates to form fractures with obvious roughness plotted in Figure 1(c). And greater loading rate causes stress concentration to form fractures with smaller waviness and roughness; otherwise, the fractures with rougher surface are formed normally.

Acquisition and Analysis of 3-D Morphology of the Rough
Fractures. Four different groups of splitting sandstone fractures mentioned above with different roughness named S1, S2, S3, and S4 are selected for laser scanning system to get the geometric morphology. And the laser scanning system is composed by the 3-dimensional noncontact optical scanner named OKIO-B and corresponding scanning software, so it can be taken to scan the fracture surface contactlessly to obtain high precision and large amount of data and exhibit the real-time photographs to get the point cloud data of the 3-dimensional fracture surface. Furthermore, the point cloud data based on the scanning system are optimized by noise reduction and patching and then imported to the software Surfer to reconstruct the 3-dimensional fracture surface [24] plotted in Figures 2(a)-2(d) (unit: mm).
In addition, the above point cloud data is changed into Cartesian coordinates, suppose Z T , Z B as the coordinates of the upper surface and the lower surface in Z direction considering the upper surface and lower surface of the fractures in the same coordinates system, as in the following equations.
where, Z Tij , Z Bij are the coordinates in Z direction at nodes (i, j) in a common coordinate system for upper and lower fracture surfaces, respectively. So, the fracture aperture (F r ) can be written by the following equations.
where b ij is the fracture aperture at nodes (i, j), which is the subtraction of Z Tij , Z Bij shown in the following equation.
So the average aperture b m of the 3-dimensional fracture surface can be written by

Geofluids
And then, based on the point cloud data calculated by Equations (1) and (2), the average surface Z coordinates (Z T ) and absolute roughness (Δ T ) may be, respectively, obtained by Equations (7) and (8) for upper fracture surfaces. Similarly, the absolute roughness of the lower fracture surface (Δ B ) can be obtained.
Thus, the average absolute roughness of fracture surface (Δ)for upper and lower fracture surfaces can be calculated by Equation (9) .The average aperture, absolute average roughness, and relative roughness for the samples  7 Geofluids S1, S2, S3, and S4 are shown in Table 1.    Figure 6: Wall friction coefficient of rough fractures for specimens S1-S4. 8 Geofluids Figure 3. The triaxial cell is capable of performing triaxial compression tests at confining pressures (P2) up to 60 MPa, with increasing deviatoric stress (P1) up to 500 MPa, with increasing transducer that has a resolution of 0.01 MPa. The system can deal with the constant head, constant flow rate, and transient pulse permeability tests under low or high confining and water pressures, and a high-precision electronic balance is selected to measure the real-time flow rate [27]. The servocontrolled fluid pump can produce pore pressure up to 60 MPa (P3/P4). Furthermore, the upstream and downstream fluid pressures can be regulated with pore pressure pumps P3 and P4; as a result, the seepage tests can be performed on constant fluid pressure or constant volume condition according to the experimental target. As for the flow procedure, the sandstone specimens are firstly saturated, then vacuumized for 4 hours by the vacuum pump, wet pumped for 4 hours with distilled water, and finally soaked for 16 hours to ensure the specimens are full of water. And then, the specimens are enclosed in a 3 mm thick Viton rubber jacket and then placed in the sample assembly. Porous spacers are inserted on to the ends of the samples to ensure even distribution of pore pressure over the ends of the samples. Considering the influence of the temperature on deformation and seepage response, all the tests are all performed at room temperature (25 ± 2°C). The seepage tests under different load combinations can be performed as follows. The samples are firstly applied with the confining pressure of the desired value, and in this stage, the axial stress is proportionally increased to the value of the confining pressure, bringing the samples to an initial isotropic stress or zero deviatoric stress to ensure that there is no gap between the specimen and the rubber jacket to prevent the fluid leakage. Furthermore, the saturated specimens are conducted at a constant pore pressure indicating the balance of the upstream pressure (P3) and downstream pressure (P4) to ensure the fluid in a single phase in the process of the seepage testing. Afterwards, water is introduced into the fractures through the inlet and collected at the outlet. A series of flow pressures at the inlet were tested and their corresponding flow rates were recorded. Therefore, the relationships between flow rate and pressure gradient or hydraulic gradient for the rough fractures can be obtained.

Variation of the Pressure Gradient and Flow Rate of the Fractures.
The testing data about the pressure gradient ∇P and the flow rate Q at the outlet for the specimens S1-S4 has been obtained and fitted into ∇P~Q curves plotted in Figure 4. It can be observed from the comparisons plotted in Figures 4(a)-4(d) of ∇P~Q about the specimens S1-S4 (black curves) and the ∇P~Q about the parallel plate model (blue curves) that the flow behavior through the rough fractures has changed obviously. These results suggest that pressure gradient ∇Pincreases linearly initially and then nonlinearly inclined to the flow rate axis as the flow rate Q increases. It can be concluded that rough fractures with smaller pressure gradient, the fluid flow mainly shows viscous forces, then inertial forces when increasing the pressure gradient, especially the roughness of the fractures surface causing obvious friction effect increases the pressure loss to transform the  9 Geofluids linear flow to the nonlinear flow. In addition, it can be seen from the mentioned variation that in view of different average apertures and absolute roughness, the pressure gradient corresponding to the ∇P~Q deviating from the linear curves is also different, meaning that the critical pressure gradient is different for different rough fractures, indicating that the surface morphology is the key factor to influence on the seepage characteristics.

Nonlinear Characteristics of the Flow in Rough Fractures.
In can be seen from the curves ∇P~Q plotted in Figures 4(a)-4(d) that increasing inertial forces has controlled the seepage behavior while increasing the pressure gradient, so the relationship of the curves ∇P~Q is nonlinear and the cubic law only considering the viscous forces cannot be suitable to describe the flow nonlinear through the rough fractures. Therefore, the equation named Navier-Stokes (N-S)   10 Geofluids describing the complex flow is introduced to analyze the fracture flow. However, the solution about the complete N-S equation is very difficult; some researchers have proposed some rational models to quantitatively present the flow nonlinear caused by the flow inertial forces, and the very common model named the Forchheimer model [10,11] is as follows: where AQ is the linear item by viscous forces, BQ 2 is the nonlinear item by inertial forces, dP/dL is the pressure gradient (Pa/m), Q is the flow rate through the fractures, A, B represent the linear coefficients (kg/(m 5 ·s)) and nonlinear coeffi- , wis the fracture width (m), μ is the dynamic viscosity, and b is the equivalent hydraulic aperture (m). In addition, the Forchheimer model has been proposed to fit corresponding data of ∇P and Q for specimens S1-S4 based on flow tests; the fitting expressions with correlation coefficient 0.99 show a good agreement with the experimental results in Section 3.1, indicating that the Forchheimer model is used to represent the flow linearity considering viscous forces with smaller pressure gradient and flow nonlinearity considering inertial forces with greater pressure gradient, revealing the complete transformation of the linear flow to nonlinear flow in the rough fractures.
3.3. The Transformation of the Flow Behavior through the Rough Fractures. It can be seen that the curves ∇P~Q that the seepage behavior can be transformed from the linearity to the nonlinearity, so the transformation of the seepage behavior should be deeply investigated to apply in the related geotechnical engineering. Reynolds number Re written by Equation (11) is an important parameter to describe the ratio of the inertial forces and the viscous forces, and the parameter reveals that Re will increase with greater flow velocity causing greater inertial forces to reflect the variation of the fluid density, viscous parameters, flow velocity, and flow paths, representing the transformation points of the seepage behavior while flowing in the rough fractures.
where D h is the hydraulic diameter (m), D h = 2b.
For determining the influence of the nonlinear item based on Forchheimer formula, non-Darcy coefficient E written by Equation (12) proposed by Zeng and Grigg [26] is introduced to explain the mechanism of nonlinear flow in the rough fractures.
where E represents the intensity influenced by the flow nonlinearity. And in the geotechnical field, the inertial forces cannot be ignored while the parameter E reaches 10% [5], indicating the nonlinear flow behavior has been triggered. Thus, the critical Reynolds number Re c based on the Equations (11) and (12) can be derived as follows [31]:

Geofluids
Therefore, non-Darcy coefficients considering different testing conditions can be calculated based on fitting coefficients A, B, and the curves E~Q and E~Re are plotted in Figure 5. It can be seen from the figure that E increases with increasing Reynolds number or flow rate, and greater inertial forces cause greater loss of the nonlinear pressure gradient, indicating that nonlinear flow in rough fractures is more intense with smooth amplification, which shows a good agreement with the conclusions in Zhou et al.. [21]. Furthermore, the critical Re c of the specimens S1-S4 are also calculated as 48.8, 25.5, 37.9, and 20.3; it can be seen that E increases with rougher fractures surface considering some Reynolds number, indicating that greater loss of the pressure gradient will cause smaller critical Reynolds number Re c ; the reason is mainly that the rougher fractures caused larger tortuosity of the seepage paths and stronger inertial forces to trigger the nonlinear flow more easily. Therefore, different roughness and aperture will generate different seepage paths, resulting in different transformation points for the linear flow to nonlinear flow in the rough fractures.
It can be also concluded that nonlinear flow through the rough fractures is closely related to the wall morphology, and increasing pressure gradient causes greater inertial forces to strengthen the friction effect due to the contact of the fluid and rough wall. Then, the wall friction coefficient f can be written by Equations (15) and (16) based on Darcy-Weisbach equation shown by Equation (14), and mentioned by Forchheimer formula, the formulas can be applied to represent the influence of the fracture roughness on the flow behavior and reveal the mechanism of nonlinear flow through rough fractures.
And then, the experimental data of ∇P~Q are used for calculating the wall friction coefficients f based on Equations (15) and (16) to build curves of f~Re plotted in Figure 6. It can be seen from the figure that smaller flow velocity will cause negative relationship of wall friction coefficient f and Reynolds number Re in the initial stage. However, the nonlinear relationship of f and Re is obvious when increasing the flow rate. And the critical wall friction coefficient f can be obtained as 1.44, 1.65, 1.78, and 1.96, indicating that the seepage state is from the linear to nonlinear; f exceeds the above values, and variation of f is induced by the variation of the seepage state; however, the seepage state is caused by velocity variation resulting from the pressure gradient, so f is only an assessment index. And the results also show that rougher fracture surface will result in larger wall friction coefficient. Also, it can be observed from the local enlarged figure plotted in Figure 6 that different fracture roughness cause obvious variation of the curves f~Re, deeply illuminating that different surface morphologies cause tortuous seepage paths [23] resulting in stronger wall friction effect to change the seepage behavior.

The Numerical Simulation of Flow Behavior through 3-Dimensional Rough Fractures
For better describing the flow variation in rough fractures, 3dimensional numerical model is proposed to simulate the seepage behavior and reveal the variation of some seepage parameters such as pressure, flow velocity, and flow rate. According to the experimental and theoretical data of the rough fractures mentioned above, the software named COM-SOL Multiphysics is proposed to reconstruct the 3dimensional fracture surface of the specimens S1, S2, S3, and S4 to form numerical models considering different rough fractures. And then, the comparison has been conducted to prove the feasibility of the proposed numerical model and   (17) and (18). And then, initial boundary is set with no-flow and no-slip on the upper, lower, left, and right boundaries, and the pressure is loaded at the inlet and outlet to complete the seepage numerical model.
The continuity equation is written as follows: The Navier-Stokes equation is shown as follows: where ρ is the fluid density, valued 998.2 kg/m 3 , μ is the dynamic viscosity coefficient, valued 0.001Pa ⋅ s, P is the fluid

Calibration of the Numerical Model for Seepage in
Rough Fractures. In order to calibrate the proposed seepage model, the specimen S1 as a case study, the flow rate Q calculated by the formula Q = Ð udA under same conditions for the experiments has been obtained in Table 2. It can be seen from the comparisons that the simulated flow rate is close to the testing value (less than 5% discrepancy). And then, the simulated curves in red of ∇P and Q plotted in Figure 4 show a good agreement with the fitting curves by Forchheimer formula and smaller discrepancy with increased pressure gradient. Therefore, the proposed seepage model can be applicable to analyze the nonlinear flow through the rough fractures.  Figure 13 based on the software Tecplot, which illuminates the curve shape which obeys the parabola shape. And in addition, the curves of the flow velocity for four specimens in the middle of the fractures with maximum values [21] show unsmooth characteristics, even though in the same intersected lines, the main reason is that different roughness cause different wall friction effects resulting in the different velocity distributions, which emphasizes the ignorance of the fracture roughness and corresponding spatial distribution influencing on the flow through the rough fractures.    Observed from the velocity distribution, it can be concluded that for the same rough fracture, increasing fracture aperture will enlarge the seepage channel to decrease the relative roughness, so the wall friction coefficient decreases resulting in greater flow velocity. And also, the flow velocity shows asymmetric because of the fracture roughness changing the seepage paths, which is in accordance with the seepage characteristics through the rough fractures, deeply indicating the importance of the fracture roughness on correct description of the seepage characteristics.

Seepage Mechanism of the Rough Fractures considering Multifactor Combinations
In order to generally describe the flow behavior through the rough fractures, the seepage mechanism based on the proposed numerical model will be deeply investigated considering the combination of the fracture roughness and aperture variation, revealing the seepage mechanism of the rough fractures under the combination of multifactors.

Nonlinear Behavior of the Flow through 3-D Rough
Fractures. According the analysis mentioned above, once the fracture aperture for the same fracture changes causing the variation of the relative roughness, the curves ∇P~Q will change greatly. Considering different combinations of the roughness and aperture, the curves ∇P~Q of four specimens S1-S4 are plotted in Figures 15(a)-15(d). It can be observed from the curve variation that the curves show linear in the initial stage but nonlinear while increasing the pressure gradient, which is in accordance with the experimental curves. However, once the fracture aperture increases causing smaller relative roughness, corresponding flow rate and velocity will increase. Therefore, the seepage behavior based on combination of the roughness and aperture can be described by the Forchheimer model, and the fitting curves show good agreement with the simulated values. The calculated coefficients A, Bare listed in Table 3. The results indicate that considering the same aperture for the four specimens, increasing roughness will narrow the seepage channel to weaken the flow rate; the fitting coefficients are greater; otherwise, increasing aperture of same specimen causing smaller relative roughness causes smaller fitting coefficients and greater flow rate and velocity. And in addition, the relationship between the fitting coefficients A, Bis modeled by the power function plotted in Figures 16(a)-16(d), which shows a good agreement with the variation in the reference written by Chen et al. [10]. It can be concluded that the fitting coefficient based on the Forchheimer model considering the geometric parameters of the rough fractures can describe the nonlinear flow behavior through the rough fractures.

Triggering Mechanism of the Nonlinear Flow Behavior through 3-D Rough
Fractures. It can be concluded from the seepage analysis mentioned in Section 3 that Re c will decrease with roughness increasing considering same pressure and aperture, showing that nonlinear flow will occur with smaller Reynolds number; otherwise, Re c will be greater with aperture increasing considering same pressure and roughness, showing that nonlinear flow is difficult to occur because of smaller relative roughness. Therefore, the aperture variation causes the relative roughness (Δ/b) representing the fracture roughness changing greatly resulting in tortuous seepage paths and different interactions of the fracture wall and fluid.
Furthermore, according to the analysis about the wall friction coefficient f , another formula describing the wall friction coefficient can be obtained by Equations (19) and (20) based on the Darcy-Weisbach equation and Forchheimer formula, and it can be seen that the wall friction coefficient f is closely related to the Reynolds number, fracture aperture, and roughness.
where a 1 and c 1 are the fitting coefficients for calculating the wall friction coefficient. The best-fit curve of f and the variables Re and Δ/b is plotted in Figure 17 In addition, the wall friction coefficients considering aperture variation of specimens S1-S4 are calculated to build the curves of f and Re plotted in Figures 18(a)-18(d). It can  In addition, comparison of critical wall friction coefficients is listed in Table 3. It can be seen that f equals to 5.06, 5.57, 5.98, and 6.56 with aperture 0.50 mm, f equals to 1.95, 2.06, 2.18, and 2.41 with aperture 0.75 mm, f equals to 0.99, 1.03, 1.12, and 1.19 with aperture 1.0 mm, and f equals to 0.63, 0.67, 0.74, and 0.79 with aperture 1.25 mm. And the figure about f , Re , and relative roughness plotted in Figure 17 is setup to show that greater roughness caused the frictional resistance work and pressure loss increasing to get greater wall friction coefficient and smaller Reynolds number, and greater aperture caused smaller wall friction coefficient and greater Reynolds number. As for the former condition, the nonlinear flow occurs more easily; otherwise, the nonlinear flow for the latter condition may not occur probably. Therefore, the wall friction coefficients related to the fracture morphology, Reynolds number, and so on can explain the seepage characteristics of rough fractures and also describe the transformation points representing the occurrence of the nonlinear flow.

Nonlinear Flow Behavior of 3-D Intersected Fractures.
For deeply investigating the seepage behavior through the rough fractures, the central part of the specimens S1 and S4 (15 ≤ X ≤ 35, 30 ≤ Y ≤ 70) is selected to reveal the seepage behavior considering fractures intersection. This cooperation not only can save the computational cost but also the roughness of S1 and S4 is almost in good with the original specimens. Thus, the data based on the scanning system is introduced to the software COMSOL and forms the intersected fractures plotted in Figure 19; then, corresponding flow behavior is analyzed in detail in the following sections. The quadratic polynomial regression in the form of Equation (10) is used to fit the ∇P − Q curves in the cases of crossing fracture intersection (CFI). The high coefficients of determination of R 2 (>0.99) demonstrate that the Forchheimer law adequately describes the nonlinear flow behavior at fracture intersections. When flow rate is small (Figures 20(a)-20(b)), the inertial force is much less than the viscous force and pressure gradient increases linearly with increasing Q. In this condition, the nonlinear term BQ 2 could be ignored and the cubic law is introduced to describe the relationship between flow rate and pressure gradient. The nonlinearity and the deviations from the cubic law are    28 Geofluids enhanced as flow rate increases. According to the definition of critical Reynolds number in Equation (13), the Re c is ranging from 85.94 to 133.42 listed in Table 4. Figures 21 and 22 deeply researched on the flow characteristics at fracture intersections and their influences on the nonlinear flow behavior for two types of rough intersected fractures. For crossing angle α of 30°, the flow streamlines follow closely the curvatures of the fracture walls. For crossing angles α of 60°changing to 150°, two eddies emerge in lower parts of the fracture which cause the reduction in the effective area available for flow when Reynolds number Re is 100. This result is consistent with that of Liu et al. [30].
To systematically investigate the nonlinear flow behavior in single fracture intersections, many extended numerical models beyond the flow experimental cases were established. Models have been established with crossing angle α of 30°,   31 Geofluids 60°, 90°, 120°, and 150°. The numerical simulations by solving the Navier-Stokes equations were also performed by COM-SOL Multiphysics. In order to better understand the influence of fracture intersections on the nonlinear flow behavior, the Forchheimer coefficient B listed in Table 4 (mentioned in Section 3.2: B = βρ/w 2 b 2 ) is calculated and the values for different cases. The relationships between B and intersected angle are also listed in Table 4. Hence, it is also proved that the intersected angle has significant influences on nonlinear flow behavior of fracture intersections.

Geofluids
Therefore, the formation of eddies is related to the fracture intersected angle, which demonstrates that the intersected angle plays a significant role in the nonlinear behavior of flow through the intersected fractures.

The Fracture Roughness Influencing on the Seepage
Characteristics. In addition, in order to describe the flow characteristics through fractures with different roughness considering some intersected angles, the simulations are plotted in Figure 23. The results indicate that the curves of ∇P − Q increased linearly initially and then nonlinearly as flow rate increased, which show a good agreement with the variation in Figure 20. Observed from the curves considering same intersection angle plotted in Figure 24 described by histogram, the variation of the nonlinear fluid coefficients B increased with the greater roughness. Therefore, the fracture surface roughness and intersection angle are both important factors influencing on the flow behavior of rough fractures.

Conclusions
In this work, flow tests have been carried out on single rough fractures to study the nonlinear behavior of fluid flow through rough fractures. Corresponding extended numerical simulations have been conducted to supplement the analyses and also analysis about influence of the fracture intersections on seepage behavior transformation of the fluid flow. The main conclusions are drawn as follows: (1) The flow tests through the rough fractures of specimens S1-S4 show that increasing pressure gradient can enlarge the flow velocity and weaken the viscous forces, and then increasing inertial forces cause the transformation of the seepage behavior from the linearity to nonlinearity. And also, the non-Darcy coefficient is introduced to explain the transformation points of the flow behavior, and the wall friction coefficient is derived based on equation Darcy-Weisbach and Forchheimer formula to reveal the transformation mechanism of fluid flow through the rough fractures (2) The 3-D fractures are reconstructed based on the scanning point cloud data, and the seepage numerical model using the Navier-Stokes equation based on the software COMSOL Multiphysics is proposed considering different roughness about four specimens S1-S4. The comparison of the ∇P~Q calibrates the feasibility of the above model (3) It can be seen from the distribution of the pressure gradient and velocity of mentioned fractures based on the proposed model that pressure decreases along the flow direction; however, greater pressure causes uneven pressure distribution and obvious nonlinear flow behavior. And also, the velocity shape is close to the parabola shape along the fracture height direction (aperture direction); however, the velocity shape shows asymmetry because of nonlinear flow behavior through rough fractures (4) The curves of the ∇P~Q considering combination of the aperture variation indicate that greater absolute roughness weakens the flow rate to enlarge the Forchheimer coefficients; however, greater aperture causes larger velocity and flow rate to get smaller Forchheimer coefficients, which reveals the relationship of the coefficients A, B,aperture and the flow nonlinearity (5) The triggering mechanism of the transformation of the flow behavior from the linearity to the nonlinearity is deeply researched, and another expressions of the wall friction coefficient have been derived to calculate f . It can be concluded that greater roughness causes larger f and smaller Re considering same apertures resulting in nonlinear flow; otherwise, increasing aperture representing smaller relative causes smaller f and larger Re considering same roughness resulting in nonlinear flow more difficultly. The results indicate that the wall friction coefficients related to the fracture morphology, Reynolds number, and so on can explain the seepage characteristics of rough fractures and also describe the