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

Authors: Xianshan Liu; Man Li; Nandou Zeng; Tao Li
Journal: Geofluids
 (2020)

---

## Abstract

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.


---

## Body

## 1. 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 3-dimensional 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–29] is not applicable to the case of fracture 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 numerically, 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 3-dimensional rough fractures considering the space distribution has not been deeply investigated. Taking a reservoir sandstone as a case study, fractures with different roughness are formed and scanned the point cloud data to obtain the key parameters describing the rough fracture. And then, flow tests have been conducted on these fracture models and the corresponding numerical simulations have been carried out by solving the Navier-Stokes equations based on software COMSOL Multiphysics to characterize the nonlinear relationship and flow distribution behavior. And also, groups of 3-dimensional fracture intersections have been studied with particular focus on the nonlinear flow behavior. And semiempirical expressions for the nonlinear term in the Forchheimer law have been proposed for different cases of rough fractures and fracture intersections, revealing the flow nonlinearity of 3-D rough fractures. The achievements can help to verify all kinds of seepage theories, especially promoting the general application in corresponding engineering.


## 2. Flow Tests of Rock Rough Fractures

### 2.1. Preparation of Rock Rough Fractures

The testing sandstones taken from western Xindianzi anticline in Chongqing City are selected for the experimental study considering different seepage designs. It is observed from Figure1(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.Figure 1
Arrangement of rock rough fractures.
(a)
Sandstones samples(b)
Brazilian test(c)
Fracture surface after the Brazilian testing


### 2.2. 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).Figure 2
3-dimensional morphology reconstruction of rough fracture of specimens S1-S4.
(a)
Specimen S1:Δ=0.58 mm(b)
Specimen S2:Δ=0.86 mm(c)
Specimen S3:Δ=1.01 mm(d)
Specimen S4:Δ=1.59 mmIn addition, the above point cloud data is changed into Cartesian coordinates, supposeZT,ZB 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.(1)ZT=ZT11⋯ZT1n⋮⋱⋮ZTm1⋯ZTmn,(2)ZB=ZB11⋯ZB1n⋮⋱⋮ZBm1⋯ZBmn,where, ZTij,ZBij 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 (Fr) can be written by the following equations.(3)Fr=ZT−ZB,(4)Fr=b11⋯b1n⋮⋱⋮bm1⋯bmn,where bij is the fracture aperture at nodes (i,j), which is the subtraction of ZTij,ZBijshown in the following equation.(5)bij=ZTij−ZBij≥0.So the average aperturebmof the 3-dimensional fracture surface can be written by(6)bm=1m×n∑i=1m∑j=1nbij.And then, based on the point cloud data calculated by Equations (1) and (2), the average surface Z coordinates (ZT) 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.(7)Z¯T=1m+1n+1∑i=0m∑j=0nZTij,(8)ΔT=1m+1n+1∑i=0m∑j=0nZTij−Z¯T.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 S1, S2, S3, and S4 are shown in Table 1.(9)Δ=ΔT+ΔB2.Table 1
Geometric parameters of the rough fractures.
Specimen no.LengthL (mm)DiameterD (mm)Average aperturebm (mm)Absolute roughnessΔ (mm)Relative roughness (Δ/bm)S1100.0450.080.870.580.67S299.9450.041.260.860.68S3100.0850.010.641.011.58S4100.0249.980.991.591.61


### 2.3. The Procedure of the Flow Test about the Rock Rough Fractures

All tests will be carried out on a rock servocontrolled triaxial equipment named Rock 600-50HT PLUS manufactured by TOP INDUSTRIE in France plotted in Figure3. 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.Figure 3
Triaxial multifield coupling test system.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.


## 2.1. Preparation of Rock Rough Fractures

The testing sandstones taken from western Xindianzi anticline in Chongqing City are selected for the experimental study considering different seepage designs. It is observed from Figure1(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.Figure 1
Arrangement of rock rough fractures.
(a)
Sandstones samples(b)
Brazilian test(c)
Fracture surface after the Brazilian testing


## 2.2. 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).Figure 2
3-dimensional morphology reconstruction of rough fracture of specimens S1-S4.
(a)
Specimen S1:Δ=0.58 mm(b)
Specimen S2:Δ=0.86 mm(c)
Specimen S3:Δ=1.01 mm(d)
Specimen S4:Δ=1.59 mmIn addition, the above point cloud data is changed into Cartesian coordinates, supposeZT,ZB 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.(1)ZT=ZT11⋯ZT1n⋮⋱⋮ZTm1⋯ZTmn,(2)ZB=ZB11⋯ZB1n⋮⋱⋮ZBm1⋯ZBmn,where, ZTij,ZBij 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 (Fr) can be written by the following equations.(3)Fr=ZT−ZB,(4)Fr=b11⋯b1n⋮⋱⋮bm1⋯bmn,where bij is the fracture aperture at nodes (i,j), which is the subtraction of ZTij,ZBijshown in the following equation.(5)bij=ZTij−ZBij≥0.So the average aperturebmof the 3-dimensional fracture surface can be written by(6)bm=1m×n∑i=1m∑j=1nbij.And then, based on the point cloud data calculated by Equations (1) and (2), the average surface Z coordinates (ZT) 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.(7)Z¯T=1m+1n+1∑i=0m∑j=0nZTij,(8)ΔT=1m+1n+1∑i=0m∑j=0nZTij−Z¯T.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 S1, S2, S3, and S4 are shown in Table 1.(9)Δ=ΔT+ΔB2.Table 1
Geometric parameters of the rough fractures.
Specimen no.LengthL (mm)DiameterD (mm)Average aperturebm (mm)Absolute roughnessΔ (mm)Relative roughness (Δ/bm)S1100.0450.080.870.580.67S299.9450.041.260.860.68S3100.0850.010.641.011.58S4100.0249.980.991.591.61


## 2.3. The Procedure of the Flow Test about the Rock Rough Fractures

All tests will be carried out on a rock servocontrolled triaxial equipment named Rock 600-50HT PLUS manufactured by TOP INDUSTRIE in France plotted in Figure3. 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.Figure 3
Triaxial multifield coupling test system.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.


## 3. Analysis of Seepage Characteristic of 3-Dimensional Rough Fractures

### 3.1. 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.Figure 4
Curves of∇P−Q for specimens S1-S4.
(a)
Specimen S1:Δ=0.58 mm(b)
Specimen S2:Δ=0.86 mm(c)
Specimen S3:Δ=1.01 mm(d)
Specimen S4:Δ=1.59 mmIt can be observed from the comparisons plotted in Figures4(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 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.


### 3.2. 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) 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:(10)−dPdL=AQ+BQ2,where AQ is the linear item by viscous forces, BQ2 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/(m5·s)) and nonlinear coefficients (kg/m8), A=12μ/wb3, B=βρ/w2b2, ρ is the fluid density (kg/m3), 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.(11)Re=ρuDhμ=2ρQμw,where Dhis the hydraulic diameter (m), Dh=2b.For determining the influence of the nonlinear item based on Forchheimer formula, non-Darcy coefficientEwritten by Equation (12) proposed by Zeng and Grigg [26] is introduced to explain the mechanism of nonlinear flow in the rough fractures.(12)E=BQ2AQ+BQ2,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 Rec based on the Equations (11) and (12) can be derived as follows [31]:(13)Rec=2AρEBμw1−E.Therefore, non-Darcy coefficients considering different testing conditions can be calculated based on fitting coefficientsA,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 Recof 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 Rec; 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.Figure 5
Variation of non-Darcy coefficient of fracture surface for samples S1-S4.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 coefficientf 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.(14)−∇P=f1Dhρu22=fρQ24b3w2,(15)f=−4∇Pb3w2ρQ2,(16)f=4b3w2ρQ2AQ+BQ2.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 fand 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; fexceeds 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.Figure 6
Wall friction coefficient of rough fractures for specimens S1-S4.


## 3.1. 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.Figure 4
Curves of∇P−Q for specimens S1-S4.
(a)
Specimen S1:Δ=0.58 mm(b)
Specimen S2:Δ=0.86 mm(c)
Specimen S3:Δ=1.01 mm(d)
Specimen S4:Δ=1.59 mmIt can be observed from the comparisons plotted in Figures4(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 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.


## 3.2. 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) 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:(10)−dPdL=AQ+BQ2,where AQ is the linear item by viscous forces, BQ2 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/(m5·s)) and nonlinear coefficients (kg/m8), A=12μ/wb3, B=βρ/w2b2, ρ is the fluid density (kg/m3), 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.(11)Re=ρuDhμ=2ρQμw,where Dhis the hydraulic diameter (m), Dh=2b.For determining the influence of the nonlinear item based on Forchheimer formula, non-Darcy coefficientEwritten by Equation (12) proposed by Zeng and Grigg [26] is introduced to explain the mechanism of nonlinear flow in the rough fractures.(12)E=BQ2AQ+BQ2,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 Rec based on the Equations (11) and (12) can be derived as follows [31]:(13)Rec=2AρEBμw1−E.Therefore, non-Darcy coefficients considering different testing conditions can be calculated based on fitting coefficientsA,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 Recof 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 Rec; 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.Figure 5
Variation of non-Darcy coefficient of fracture surface for samples S1-S4.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 coefficientf 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.(14)−∇P=f1Dhρu22=fρQ24b3w2,(15)f=−4∇Pb3w2ρQ2,(16)f=4b3w2ρQ2AQ+BQ2.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 fand 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; fexceeds 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.Figure 6
Wall friction coefficient of rough fractures for specimens S1-S4.


## 4. The Numerical Simulation of Flow Behavior through 3-Dimensional Rough Fractures

For better describing the flow variation in rough fractures, 3-dimensional 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 COMSOL Multiphysics is proposed to reconstruct the 3-dimensional 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 then analyzes the seepage characteristics in detail considering different conditions.


### 4.1. The Model Setup and Corresponding Calibration of the Seepage Model through the 3-D Rough Fractures

#### 4.1.1. The Model Setup for the Fracture Seepage

The software named COMSOL Multiphysics is a popular numerical software solving the multiphysics coupling problems in many fields of solid mechanics, fluid mechanics, and thermal mechanics. For the model setup of the flow in rough fractures, the scanning data about the 3-D fracture surface is introduced into the software to form the 3-D rough models shown in Figure7(a), and the complete specimens, corresponding upper surface and lower surface, are plotted in Figures 7(b) and 7(c). And also, the mesh of the fracture models is shown in Figures 8(a) and 8(b).Figure 7
Physical model of rough fractured for specimen S1.
(a)
Fracture model(b)
Upper surface(c)
Lower surfaceFigure 8
3-dimensional numerical model of sample S1 with rough surface.
(a)
3-dimensional mesh of the complete model(b)
3-dimensional mesh of the fracturesAs for 3-D seepage calculation, suppose the incompressible fluid, constant viscosity and fluid density, and the laminar flow, the continuity equation and Navier-Stokes equation are included in Equations (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:(17)∇u=0.The Navier-Stokes equation is shown as follows:(18)ρu⋅∇u=μ∇2u−∇P,where ρ is the fluid density, valued 998.2 kg/m3, μ is the dynamic viscosity coefficient, valued 0.001Pa⋅s, P is the fluid pressure (Pa), uis the flow velocity (m/s), ρu⋅∇u represents the item describing the fluid inertial forces, and μ∇2u represents the item describing the fluid viscous forces.


#### 4.1.2. 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 rateQ 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.Table 2
Comparison of flow rate between the experiments and numerical simulation of specimen S1.
Pressure gradient∇P(kPa/m)Flow rate at outletQ(×105 m3/s)Experimental valueSimulated valueRelative error10.2890.23917.3%20.4990.44910.0%30.6560.6333.5%40.8030.7941.1%50.9560.9381.9%71.1881.1890.1%91.4801.4025.3%111.6421.5893.2%131.8211.7563.6%


### 4.2. Analysis of the Seepage Behavior through 3-D Rough Fractures

#### 4.2.1. The Description of the Flow Behavior

The flow behavior of the specimens S1-S4 is deeply investigated based on the calibrated seepage model. And the pressure distribution under different pressure gradients is plotted in Figure9, and the pressure decreases from the inlet to the outlet and the pressure lines are uneven along the flow direction, indicating that rough fracture wall blocks the fluid flow forward and the curves are bended obviously to spread to the left and right while increasing the pressure gradient. Meanwhile, the pressure p0=200Pa is supposed to get the pressure distribution on the fracture surface plotted in Figures 10(a)–10(d), and different pressures are obvious for same points under same conditions, illuminating the fracture roughness plays an important role on the seepage characteristics.Figure 9
Distribution of seepage pressure on the fracture surface of specimen S1.Figure 10
Comparison of surface pressure of specimens S1-S4.
(a)
Specimen S1(b)
Specimen S2(c)
Specimen S3(d)
Specimen S4For deeply investigating the flow characteristics, cloud figures plotted in Figures11(a)–11(d) and the velocity distribution on section X=25mm of the specimens S1-S4 plotted in Figures 12(a)–12(e) are treated to describe the flow velocity on the different sections. It can be observed from the figures mentioned above that the smallest flow velocity is close to zero near the upper wall and lower wall, otherwise the greatest in the middle of the fracture, which is in accidence with the velocity variation described by the Darcy-Weisbach equation. In order to display the flow variation along the fracture height, the intersected line on the section Y=50mm and X=25mm of four specimens is selected to show the curves variation plotted in 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.Figure 11
Velocity distribution on fracture surface of specimens S1-S4.
(a)
Specimen S1(b)
Specimen S2(c)
Specimen S3(d)
Specimen S4Figure 12
Velocity distribution along the lines (X=25mm) of specimens S1-S4.
(a)
The diagram of the researched line(b)
Specimen S1(c)
Specimen S2(d)
Specimen S3(e)
Specimen S4Figure 13
Velocity distribution along the fracture height (Z direction) of specimens S1-S4 (X=25mm interested Y=50mm).


#### 4.2.2. Analysis of the Seepage Behavior considering the Aperture Variation

The seepage behavior considering the aperture variation has been deeply researched in this section. Considering the aperture with 0.50 mm, 0.75 mm, 1.00 mm, and 1.25 mm and same intersected lines on sectionX=25mm and Y=50mm, the curves of the flow velocity along the fracture height representing the aperture direction (Z direction) are compared shown in Figures 14(a)–14(d). 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.Figure 14
Comparison of flow velocity of fractures in samples S1-S4 considering varied aperture.
(a)
Specimen S1(b)
Specimen S2(c)
Specimen S3(d)
Specimen S4


## 4.1. The Model Setup and Corresponding Calibration of the Seepage Model through the 3-D Rough Fractures

### 4.1.1. The Model Setup for the Fracture Seepage

The software named COMSOL Multiphysics is a popular numerical software solving the multiphysics coupling problems in many fields of solid mechanics, fluid mechanics, and thermal mechanics. For the model setup of the flow in rough fractures, the scanning data about the 3-D fracture surface is introduced into the software to form the 3-D rough models shown in Figure7(a), and the complete specimens, corresponding upper surface and lower surface, are plotted in Figures 7(b) and 7(c). And also, the mesh of the fracture models is shown in Figures 8(a) and 8(b).Figure 7
Physical model of rough fractured for specimen S1.
(a)
Fracture model(b)
Upper surface(c)
Lower surfaceFigure 8
3-dimensional numerical model of sample S1 with rough surface.
(a)
3-dimensional mesh of the complete model(b)
3-dimensional mesh of the fracturesAs for 3-D seepage calculation, suppose the incompressible fluid, constant viscosity and fluid density, and the laminar flow, the continuity equation and Navier-Stokes equation are included in Equations (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:(17)∇u=0.The Navier-Stokes equation is shown as follows:(18)ρu⋅∇u=μ∇2u−∇P,where ρ is the fluid density, valued 998.2 kg/m3, μ is the dynamic viscosity coefficient, valued 0.001Pa⋅s, P is the fluid pressure (Pa), uis the flow velocity (m/s), ρu⋅∇u represents the item describing the fluid inertial forces, and μ∇2u represents the item describing the fluid viscous forces.


### 4.1.2. 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 rateQ 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.Table 2
Comparison of flow rate between the experiments and numerical simulation of specimen S1.
Pressure gradient∇P(kPa/m)Flow rate at outletQ(×105 m3/s)Experimental valueSimulated valueRelative error10.2890.23917.3%20.4990.44910.0%30.6560.6333.5%40.8030.7941.1%50.9560.9381.9%71.1881.1890.1%91.4801.4025.3%111.6421.5893.2%131.8211.7563.6%


## 4.1.1. The Model Setup for the Fracture Seepage

The software named COMSOL Multiphysics is a popular numerical software solving the multiphysics coupling problems in many fields of solid mechanics, fluid mechanics, and thermal mechanics. For the model setup of the flow in rough fractures, the scanning data about the 3-D fracture surface is introduced into the software to form the 3-D rough models shown in Figure7(a), and the complete specimens, corresponding upper surface and lower surface, are plotted in Figures 7(b) and 7(c). And also, the mesh of the fracture models is shown in Figures 8(a) and 8(b).Figure 7
Physical model of rough fractured for specimen S1.
(a)
Fracture model(b)
Upper surface(c)
Lower surfaceFigure 8
3-dimensional numerical model of sample S1 with rough surface.
(a)
3-dimensional mesh of the complete model(b)
3-dimensional mesh of the fracturesAs for 3-D seepage calculation, suppose the incompressible fluid, constant viscosity and fluid density, and the laminar flow, the continuity equation and Navier-Stokes equation are included in Equations (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:(17)∇u=0.The Navier-Stokes equation is shown as follows:(18)ρu⋅∇u=μ∇2u−∇P,where ρ is the fluid density, valued 998.2 kg/m3, μ is the dynamic viscosity coefficient, valued 0.001Pa⋅s, P is the fluid pressure (Pa), uis the flow velocity (m/s), ρu⋅∇u represents the item describing the fluid inertial forces, and μ∇2u represents the item describing the fluid viscous forces.


## 4.1.2. 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 rateQ 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.Table 2
Comparison of flow rate between the experiments and numerical simulation of specimen S1.
Pressure gradient∇P(kPa/m)Flow rate at outletQ(×105 m3/s)Experimental valueSimulated valueRelative error10.2890.23917.3%20.4990.44910.0%30.6560.6333.5%40.8030.7941.1%50.9560.9381.9%71.1881.1890.1%91.4801.4025.3%111.6421.5893.2%131.8211.7563.6%


## 4.2. Analysis of the Seepage Behavior through 3-D Rough Fractures

### 4.2.1. The Description of the Flow Behavior

The flow behavior of the specimens S1-S4 is deeply investigated based on the calibrated seepage model. And the pressure distribution under different pressure gradients is plotted in Figure9, and the pressure decreases from the inlet to the outlet and the pressure lines are uneven along the flow direction, indicating that rough fracture wall blocks the fluid flow forward and the curves are bended obviously to spread to the left and right while increasing the pressure gradient. Meanwhile, the pressure p0=200Pa is supposed to get the pressure distribution on the fracture surface plotted in Figures 10(a)–10(d), and different pressures are obvious for same points under same conditions, illuminating the fracture roughness plays an important role on the seepage characteristics.Figure 9
Distribution of seepage pressure on the fracture surface of specimen S1.Figure 10
Comparison of surface pressure of specimens S1-S4.
(a)
Specimen S1(b)
Specimen S2(c)
Specimen S3(d)
Specimen S4For deeply investigating the flow characteristics, cloud figures plotted in Figures11(a)–11(d) and the velocity distribution on section X=25mm of the specimens S1-S4 plotted in Figures 12(a)–12(e) are treated to describe the flow velocity on the different sections. It can be observed from the figures mentioned above that the smallest flow velocity is close to zero near the upper wall and lower wall, otherwise the greatest in the middle of the fracture, which is in accidence with the velocity variation described by the Darcy-Weisbach equation. In order to display the flow variation along the fracture height, the intersected line on the section Y=50mm and X=25mm of four specimens is selected to show the curves variation plotted in 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.Figure 11
Velocity distribution on fracture surface of specimens S1-S4.
(a)
Specimen S1(b)
Specimen S2(c)
Specimen S3(d)
Specimen S4Figure 12
Velocity distribution along the lines (X=25mm) of specimens S1-S4.
(a)
The diagram of the researched line(b)
Specimen S1(c)
Specimen S2(d)
Specimen S3(e)
Specimen S4Figure 13
Velocity distribution along the fracture height (Z direction) of specimens S1-S4 (X=25mm interested Y=50mm).


### 4.2.2. Analysis of the Seepage Behavior considering the Aperture Variation

The seepage behavior considering the aperture variation has been deeply researched in this section. Considering the aperture with 0.50 mm, 0.75 mm, 1.00 mm, and 1.25 mm and same intersected lines on sectionX=25mm and Y=50mm, the curves of the flow velocity along the fracture height representing the aperture direction (Z direction) are compared shown in Figures 14(a)–14(d). 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.Figure 14
Comparison of flow velocity of fractures in samples S1-S4 considering varied aperture.
(a)
Specimen S1(b)
Specimen S2(c)
Specimen S3(d)
Specimen S4


## 4.2.1. The Description of the Flow Behavior

The flow behavior of the specimens S1-S4 is deeply investigated based on the calibrated seepage model. And the pressure distribution under different pressure gradients is plotted in Figure9, and the pressure decreases from the inlet to the outlet and the pressure lines are uneven along the flow direction, indicating that rough fracture wall blocks the fluid flow forward and the curves are bended obviously to spread to the left and right while increasing the pressure gradient. Meanwhile, the pressure p0=200Pa is supposed to get the pressure distribution on the fracture surface plotted in Figures 10(a)–10(d), and different pressures are obvious for same points under same conditions, illuminating the fracture roughness plays an important role on the seepage characteristics.Figure 9
Distribution of seepage pressure on the fracture surface of specimen S1.Figure 10
Comparison of surface pressure of specimens S1-S4.
(a)
Specimen S1(b)
Specimen S2(c)
Specimen S3(d)
Specimen S4For deeply investigating the flow characteristics, cloud figures plotted in Figures11(a)–11(d) and the velocity distribution on section X=25mm of the specimens S1-S4 plotted in Figures 12(a)–12(e) are treated to describe the flow velocity on the different sections. It can be observed from the figures mentioned above that the smallest flow velocity is close to zero near the upper wall and lower wall, otherwise the greatest in the middle of the fracture, which is in accidence with the velocity variation described by the Darcy-Weisbach equation. In order to display the flow variation along the fracture height, the intersected line on the section Y=50mm and X=25mm of four specimens is selected to show the curves variation plotted in 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.Figure 11
Velocity distribution on fracture surface of specimens S1-S4.
(a)
Specimen S1(b)
Specimen S2(c)
Specimen S3(d)
Specimen S4Figure 12
Velocity distribution along the lines (X=25mm) of specimens S1-S4.
(a)
The diagram of the researched line(b)
Specimen S1(c)
Specimen S2(d)
Specimen S3(e)
Specimen S4Figure 13
Velocity distribution along the fracture height (Z direction) of specimens S1-S4 (X=25mm interested Y=50mm).


## 4.2.2. Analysis of the Seepage Behavior considering the Aperture Variation

The seepage behavior considering the aperture variation has been deeply researched in this section. Considering the aperture with 0.50 mm, 0.75 mm, 1.00 mm, and 1.25 mm and same intersected lines on sectionX=25mm and Y=50mm, the curves of the flow velocity along the fracture height representing the aperture direction (Z direction) are compared shown in Figures 14(a)–14(d). 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.Figure 14
Comparison of flow velocity of fractures in samples S1-S4 considering varied aperture.
(a)
Specimen S1(b)
Specimen S2(c)
Specimen S3(d)
Specimen S4


## 5. 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.


### 5.1. 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.Figure 15
Curves of∇P and Q considering the variation of the hydroaperture for samples.
(a)
Fracture apertureb=0.50 mm(b)
Fracture apertureb=0.75 mm(c)
Fracture apertureb=1.00 mm(d)
Fracture apertureb=1.25 mmTherefore, 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 coefficientsA,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.Table 3
Fitting coefficientsA and B in the Forchheimer equation and corresponding critical Reynolds numbers.
Fracture apertureb (mm)Specimen S1Specimen S2Specimen S3Specimen S4A (108)B (1013)fA (108)B (1013)fA (108)B (1013)fA (108)B (1013)f0.5020.0859.6885.0620.45610.6685.5720.59111.4595.9820.8212.5736.560.754.863.7331.955.1173.9532.065.3144.1842.185.5884.612.411.001.4621.8950.991.5121.9821.031.5742.1461.121.6272.2851.191.250.4151.2010.630.4371.2820.670.4761.4100.740.5031.5110.79Figure 16
Fitting curves of the coefficientsA and B based on Forchheimer equation and fracture aperture for samples.
(a)
Specimen S1(b)
Specimen S2(c)
Specimen S3(d)
Specimen S4


### 5.2. Triggering Mechanism of the Nonlinear Flow Behavior through 3-D Rough Fractures

It can be concluded from the seepage analysis mentioned in Section3 that Recwill decrease with roughness increasing considering same pressure and aperture, showing that nonlinear flow will occur with smaller Reynolds number; otherwise, Recwill 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 coefficientf, 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.(19a)f=96Re+4bβ,(19b)4bβ=a1Δbc1,where a1 and c1 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, which indicates a strong relationship. The best regression values of the coefficients in Equation (20) are a1=0.067 and c1=−1.128 with coefficient of determination R2=0.944.(20)f=96Re+0.067Δb−1.128.Figure 17
Relationship off~Re and relative roughness.In addition, the wall friction coefficients considering aperture variation of specimens S1-S4 are calculated to build the curves off and Re plotted in Figures 18(a)–18(d). It can be seen that for some specimen, the curves of f and Re are linear with smaller Reynolds number and then show nonlinear characteristics and deviate the curve considering the parallel plate model. And the conclusion can be done that greater relative roughness causes obvious flow nonlinearity, and f is controlled by the Reynolds number with smaller flow velocity; otherwise, f is controlled by multifactors with greater flow velocity.Figure 18
Curves of wall friction factor and Reynolds.
(a)
b=0.50 mm(b)
b=0.75 mm(c)
b=1.00 mm(d)
b=1.25 mmIn addition, comparison of critical wall friction coefficients is listed in Table3. 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.


### 5.3. 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.Figure 19
The intersection of rough fractures.
(a)
The diagram of the interested fractures(b)
The intersected fractures of S1 with angle 30°(c)
The intersected fractures of S4 with angle 30°


#### 5.3.1. Seepage Flow considering Different Interested Angles of the Rough Fractures

The experimental results of curves∇P−Q for the specimens S1 and S4 are shown in Figures 20(a)–20(b). These results suggest that pressure gradient increases linearly initially and then nonlinearly as flow rate increases. The data in the linear stage is used to calculate hydraulic apertures. The aperture of the corresponding numerical model uses the same value of hydraulic aperture obtained by experimental data. Simulations are also plotted in Figures 21 and 22.Figure 20
Curves of pressure gradient and flowrate.
(a)
Specimen S1(b)
Specimen S4Figure 21
Velocity distribution of specimen S1.
(a)
Interested angle of 30°(b)
Interested angle of 60°(c)
Interested angle of 90°(d)
Interested angle of 120°(e)
Interested angle of 150°Figure 22
Velocity distribution of specimen S4.
(a)
Interested angle of 30°(b)
Interested angle of 60°(c)
Interested angle of 90°(d)
Interested angle of 120°(e)
Interested angle of 150°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 R2 (>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 BQ2 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 enhanced as flow rate increases. According to the definition of critical Reynolds number in Equation (13), the Rec is ranging from 85.94 to 133.42 listed in Table 4.Table 4
The seepage parameters considering different intersected angles of rough fractures.
CFIα/°Specimen S1Specimen S4B (kg/m8/1014)RecβB (kg/m8/1014)Recβ302.309133.420.0202.396128.610.021602.806109.810.0243.043101.230.026902.969103.780.0263.28493.830.0281203.022101.930.0263.51187.760.0301503.36991.440.0293.58585.940.031Figures21 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 Reis 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°, 60°, 90°, 120°, and 150°. The numerical simulations by solving the Navier-Stokes equations were also performed by COMSOL 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=βρ/w2b2) 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.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.


#### 5.3.2. 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 Figure23. 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.Figure 23
Curves of∇P~Q for rough fractures with different interested angles.
(a)
Interested angle of 30°(b)
Interested angle of 90°(c)
Interested angle of 150°Figure 24
The comparison of nonlinear coefficients for rough fractures with different interested angles.


## 5.1. 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.Figure 15
Curves of∇P and Q considering the variation of the hydroaperture for samples.
(a)
Fracture apertureb=0.50 mm(b)
Fracture apertureb=0.75 mm(c)
Fracture apertureb=1.00 mm(d)
Fracture apertureb=1.25 mmTherefore, 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 coefficientsA,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.Table 3
Fitting coefficientsA and B in the Forchheimer equation and corresponding critical Reynolds numbers.
Fracture apertureb (mm)Specimen S1Specimen S2Specimen S3Specimen S4A (108)B (1013)fA (108)B (1013)fA (108)B (1013)fA (108)B (1013)f0.5020.0859.6885.0620.45610.6685.5720.59111.4595.9820.8212.5736.560.754.863.7331.955.1173.9532.065.3144.1842.185.5884.612.411.001.4621.8950.991.5121.9821.031.5742.1461.121.6272.2851.191.250.4151.2010.630.4371.2820.670.4761.4100.740.5031.5110.79Figure 16
Fitting curves of the coefficientsA and B based on Forchheimer equation and fracture aperture for samples.
(a)
Specimen S1(b)
Specimen S2(c)
Specimen S3(d)
Specimen S4


## 5.2. Triggering Mechanism of the Nonlinear Flow Behavior through 3-D Rough Fractures

It can be concluded from the seepage analysis mentioned in Section3 that Recwill decrease with roughness increasing considering same pressure and aperture, showing that nonlinear flow will occur with smaller Reynolds number; otherwise, Recwill 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 coefficientf, 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.(19a)f=96Re+4bβ,(19b)4bβ=a1Δbc1,where a1 and c1 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, which indicates a strong relationship. The best regression values of the coefficients in Equation (20) are a1=0.067 and c1=−1.128 with coefficient of determination R2=0.944.(20)f=96Re+0.067Δb−1.128.Figure 17
Relationship off~Re and relative roughness.In addition, the wall friction coefficients considering aperture variation of specimens S1-S4 are calculated to build the curves off and Re plotted in Figures 18(a)–18(d). It can be seen that for some specimen, the curves of f and Re are linear with smaller Reynolds number and then show nonlinear characteristics and deviate the curve considering the parallel plate model. And the conclusion can be done that greater relative roughness causes obvious flow nonlinearity, and f is controlled by the Reynolds number with smaller flow velocity; otherwise, f is controlled by multifactors with greater flow velocity.Figure 18
Curves of wall friction factor and Reynolds.
(a)
b=0.50 mm(b)
b=0.75 mm(c)
b=1.00 mm(d)
b=1.25 mmIn addition, comparison of critical wall friction coefficients is listed in Table3. 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.


## 5.3. 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.Figure 19
The intersection of rough fractures.
(a)
The diagram of the interested fractures(b)
The intersected fractures of S1 with angle 30°(c)
The intersected fractures of S4 with angle 30°


### 5.3.1. Seepage Flow considering Different Interested Angles of the Rough Fractures

The experimental results of curves∇P−Q for the specimens S1 and S4 are shown in Figures 20(a)–20(b). These results suggest that pressure gradient increases linearly initially and then nonlinearly as flow rate increases. The data in the linear stage is used to calculate hydraulic apertures. The aperture of the corresponding numerical model uses the same value of hydraulic aperture obtained by experimental data. Simulations are also plotted in Figures 21 and 22.Figure 20
Curves of pressure gradient and flowrate.
(a)
Specimen S1(b)
Specimen S4Figure 21
Velocity distribution of specimen S1.
(a)
Interested angle of 30°(b)
Interested angle of 60°(c)
Interested angle of 90°(d)
Interested angle of 120°(e)
Interested angle of 150°Figure 22
Velocity distribution of specimen S4.
(a)
Interested angle of 30°(b)
Interested angle of 60°(c)
Interested angle of 90°(d)
Interested angle of 120°(e)
Interested angle of 150°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 R2 (>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 BQ2 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 enhanced as flow rate increases. According to the definition of critical Reynolds number in Equation (13), the Rec is ranging from 85.94 to 133.42 listed in Table 4.Table 4
The seepage parameters considering different intersected angles of rough fractures.
CFIα/°Specimen S1Specimen S4B (kg/m8/1014)RecβB (kg/m8/1014)Recβ302.309133.420.0202.396128.610.021602.806109.810.0243.043101.230.026902.969103.780.0263.28493.830.0281203.022101.930.0263.51187.760.0301503.36991.440.0293.58585.940.031Figures21 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 Reis 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°, 60°, 90°, 120°, and 150°. The numerical simulations by solving the Navier-Stokes equations were also performed by COMSOL 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=βρ/w2b2) 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.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.


### 5.3.2. 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 Figure23. 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.Figure 23
Curves of∇P~Q for rough fractures with different interested angles.
(a)
Interested angle of 30°(b)
Interested angle of 90°(c)
Interested angle of 150°Figure 24
The comparison of nonlinear coefficients for rough fractures with different interested angles.


## 5.3.1. Seepage Flow considering Different Interested Angles of the Rough Fractures

The experimental results of curves∇P−Q for the specimens S1 and S4 are shown in Figures 20(a)–20(b). These results suggest that pressure gradient increases linearly initially and then nonlinearly as flow rate increases. The data in the linear stage is used to calculate hydraulic apertures. The aperture of the corresponding numerical model uses the same value of hydraulic aperture obtained by experimental data. Simulations are also plotted in Figures 21 and 22.Figure 20
Curves of pressure gradient and flowrate.
(a)
Specimen S1(b)
Specimen S4Figure 21
Velocity distribution of specimen S1.
(a)
Interested angle of 30°(b)
Interested angle of 60°(c)
Interested angle of 90°(d)
Interested angle of 120°(e)
Interested angle of 150°Figure 22
Velocity distribution of specimen S4.
(a)
Interested angle of 30°(b)
Interested angle of 60°(c)
Interested angle of 90°(d)
Interested angle of 120°(e)
Interested angle of 150°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 R2 (>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 BQ2 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 enhanced as flow rate increases. According to the definition of critical Reynolds number in Equation (13), the Rec is ranging from 85.94 to 133.42 listed in Table 4.Table 4
The seepage parameters considering different intersected angles of rough fractures.
CFIα/°Specimen S1Specimen S4B (kg/m8/1014)RecβB (kg/m8/1014)Recβ302.309133.420.0202.396128.610.021602.806109.810.0243.043101.230.026902.969103.780.0263.28493.830.0281203.022101.930.0263.51187.760.0301503.36991.440.0293.58585.940.031Figures21 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 Reis 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°, 60°, 90°, 120°, and 150°. The numerical simulations by solving the Navier-Stokes equations were also performed by COMSOL 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=βρ/w2b2) 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.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.


## 5.3.2. 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 Figure23. 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.Figure 23
Curves of∇P~Q for rough fractures with different interested angles.
(a)
Interested angle of 30°(b)
Interested angle of 90°(c)
Interested angle of 150°Figure 24
The comparison of nonlinear coefficients for rough fractures with different interested angles.


## 6. 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 calculatef. 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 transformation points representing the occurring of the nonlinear flow(6)
The flow nonlinearity is strongly affected by the fracture intersected angle. The nonlinearity of fluid flow is mainly caused by the formation of eddies at fracture intersections, which are clearly observable in numerical models. And it can be observed from the numerical results that the critical Reynolds number decreases with increasing angle. Furthermore, the relationships based on simulations that relate the nonlinear coefficientB in Forchheimer law to intersected angle and surface roughness have been proposed to reveal the flow nonlinearity through the rough fractures



---

*Source: 8818749-2020-11-12.xml*