# Experimental Investigation on Nonlinear Flow Anisotropy Behavior in Fracture Media

Authors: Chun Zhu; Xiaoding Xu; Xiuting Wang; Feng Xiong; Zhigang Tao; Yun Lin; Jing Chen
Journal: Geofluids
 (2019)

---

## Abstract

A series of flow experiments were performed on matched fractures to study the problem of non-Darcy flow in fractured media. Five rock fractures of different roughness were generated using indirect tensile tests, and their surface topographies were measured using a stereo topometric scanning system. The fracture was assumed to be a self-affine surface, and its roughness and anisotropy were quantified by the fractal dimension. According to the flow tortuosity effect, the nonlinear flow was characterized by hydraulic tortuosity and surface tortuosity power law relationships based on Forchheimer’s law. Fracture seepage experiments conducted with two injection directions (0° and 90°) showed that Forchheimer’s law described the nonlinear flow well. Both the proposed model and Chen’s double-parameter model gave similar results to the experiment, but the match was closer with the proposed model. On this basis, a new formula for the critical Reynolds number is proposed, which serves to distinguish linear flow and Forchheimer flow.


---

## Body

## 1. Introduction

A long history of geological and human activities has caused most rock masses to be cut by a large number of faults and fractures [1–5]. These discontinuities form the main channels for groundwater flow, which control the permeability characteristics of the rock mass. In the study of rock mass hydrology, discontinuities are usually generalized into two smooth parallel plates, and the famous cubic law is hence obtained through theory and experiment. A variety of correction models [6–11] has been proposed to account for fracture roughness, contact, or filling.Some engineering projects involve a high hydraulic gradient, for example, dam construction in the deep weak overburden of a river valley, exploitation of low-permeability oil and gas wells, and coal mine gas outbursts [12–16]. Under this condition, fluid flow through fractures is not linear, and the use of the cubic law or related modified models would cause large deviations. The well-known Forchheimer law is used to describe this flow behavior:
(1)∇P=AQ+BQ2,where ∇P is the pressure gradient between the inlet and outlet of the fracture, Q is the flow rate through the fracture, and A and B are the coefficients of viscosity and inertia, respectively. Zimmerman et al. [17] observed the Forchheimer flow phenomena of rough fractures when the Reynolds number Re>20 through experiments and numerical methods. Zhang and Nemcik [18] discussed the linear and nonlinear flow characteristics of rough fractures under different confining pressures. Zhou et al. [19] explained the physical significance of the Forchheimer flow coefficients A and B and the internal transition mechanism on the basis of water pressure tests under different confining pressures. However, the effects of fracture roughness on flow were not explained in detail. Jin et al. (2013) pointed out that the influence of a rough geometry on fracture flow is manifested in three aspects: the frictional effect in the fluid, the tortuous effect of the fracture surface, and a local roughness effect. Tsang [20] considered that the roughness of the fracture surface would lead to flow tortuosity, and Xiao et al. [21] introduced a tortuosity factor to describe the tortuosity of flow. Watanabe et al. [22] carried out fluid flow experiments in fractures with shear displacements and found that the nonlinear flow effect decreases as shear increases.Fractal geometry was first put forward by B.B. Mandelbrot. Xie and Wang [23] introduced it into the description of fracture roughness and then used it to describe fluid flow characteristics at rock fracture surfaces. Murata and Saito [24] studied the influence of fractal parameters on the tortuosity effect, and Wang et al. [25] put forward a flow model using fractal parameters. Ju et al. [26] carried out flow experiments on rough single fractures with different fractal dimensions. These flow tests clarified the influence of a rough structure on seepage flow. Develi and Babadagli [27] carried out saturated seepage tests on seven kinds of artificial tensional fracture surfaces, describing the roughness of the fractures by means of the fractal dimension and discussing the influences of roughness, anisotropy, and normal stress on seepage characteristics.In this paper, a nonlinear fractal model for rough fractures is deduced based on the tortuous effect and the self-affine fractal characteristics of the fracture surface. The law and anisotropy of Forchheimer flow are analyzed, and the new model is verified by seepage tests.


## 2. Nonlinear Fractal Model for Rough-Walled Rock Fractures

Forchheimer’s law is composed of a linear partAQ and a nonlinear part BQ2. When the flow rate is low, a cubic law can be used to describe the relationship of the flow rate and pressure. Hence, A can be expressed as
(2)A=12μweh3,where μ is the dynamic coefficient of viscosity of the fluid and w is the width of the fracture. eh is the hydraulic aperture, which can be calculated as eh=12μQ/w/▽P1/3.The coefficientB represents the degree that the seepage curve deviates from that in the linear stage. Schrauf and Evans [28] put forward a form of B using dimensional analysis:
(3)B=bDρeh3w2,where ρ is the fluid density and bD is a parameter related to the roughness of the fracture surface. Chen et al. [6] used the peak asperity height to describe fracture roughness and obtained a two-parameter model for bD:
(4)bD=aξ2ehb,where a and b are fitting parameters, respectively. However, the peak asperity height does not account for flow tortuosity and anisotropy. Chen et al. [6] also used numerical simulation to study non-Darcy behavior in fracture flow. The results showed that the rougher the surface was, the more tortuous the flow would be, and eddy currents and recirculation would occur at high velocity, which would increase the inertial resistance. In order to characterize the effect of flow tortuosity, the following power law relations are proposed by Murata and Saito [24] and Ji et al. (2015):
(5)bD=cτaτsb,where a, b, and c are fitting parameters. τ is the hydrological curvature, which is defined as the ratio of the actual length Lt of the seepage path to the horizontal length Lc along the pressure gradient of the fracture. τs is the curvature of the surface, which is defined as the ratio of the surface area to the projection area of the fracture surface.For fractal fractures, the relationship between measureFδ and measurement scale δ is as follows:
(6)Fδ=F0δα,where α is a parameter related to the fractal dimension D and D is in the range 1-3. F0 is the apparent measurement. Based on this, the relationship between the fracture surface area As and a square mesh with dimension δ1 is as follows:
(7)As=F1δ12−Ds,where Ds is the fractal dimension of the fracture area. For a square fracture, when δ1 is equal to Lc and equation (7) is substituted into AsLc=Ac, F1 can be obtained as follows:
(8)F1=AcDs/2.Hence,τs can be obtained:
(9)τs=Asδ1Ac=δ12Ac2−Ds.In addition, the fractal relationship of the length of the seepage pathLt and measure δ2 is as follows:
(10)Ltδ2=F2δ21−DT,where DT is the fractal dimension of the seepage path. When δ2 is equal to Lc and equation (10) is substituted into LtLc=Lc, F2 can be obtained as
(11)F2=LcDT.Hence,τ can be obtained:
(12)τ=Ltδ2Lc=δ2Lc1−DT.Mandelbrot, the founder of fractal theory, suggested that the fractal dimension of the surface could be calculated by adding the dimension of a surface profile to 1. Therefore, the relationship between the fractal dimension of the profile length and the fractal dimension of the area is as follows:(13)Ds=Dl+1.Jin et al. [29] considered Dl as equal to DT. Hence, bD can be obtained by substituting equation (9) and equation (12) into equation (5):
(14)bD=cδ2Lca1−Dlδ12Acb1−Dl.Whenδ1 and δ2 are eh, equation (1) becomes
(15)bD=cehLca+b1−Dl.Hence, a new model of parameterB in Forchheimer’s law can be obtained:
(16)B=aehLcb1−Dlρeh3w2,where a and b are constants that can be determined with fracture permeability tests. Firstly, the fractal dimension is calculated from fracture surface data. The curve relating flow to pressure gradient is then obtained through fracture permeability tests, and this is used to obtain a and b.


## 3. Fracture Surface Measurement

### 3.1. Rock Fracture Preparation

The effect of fracture surface roughness on fluid flow was studied by way of saturated seepage tests of rock fracture surfaces of different roughness. A natural fracture surface is difficult to obtain, so artificial tension fracture specimens were used to study the characteristics of fluid flow in fractures. In this study, natural granite selected from a quarry in Sichuan Province was processed into150mm×150mm×150mm square specimens in the laboratory, and then, the specimens were split using the Brazilian splitting test method to obtain artificial tensional joint specimens. Finally, five groups of fracture surfaces (F1, F2, F3, F4, and F5) with different roughness were prepared.


### 3.2. Measurement Procedure

A portable 3D optical three-dimensional scanning system was used to measure the three-dimensional morphology of the fracture surface (Figure1). The system broadly comprises a scanning control unit, scanning lens, turn table, and tripod. The scanning lens is placed on the tripod, which can rotate freely and adjust its position conveniently. The other components are connected via USB. The system acquires a visible grating image projected onto the surface of the object then accurately determines the spatial coordinates (X, Y, Z) of each point using the phase and triangulation methods according to the shape of fringe curvature change, forming a three-dimensional point cloud data. This approach has the benefits of being fast, high precision (measuring accuracy 25 μm), and allowing noncontact measurement.Figure 1
Stereo topometric scanning system.In the actual measurement process, features of the measurement environment (light, dust, etc.) and the measurement methods will have an impact on the accuracy of the three-dimensional topographic data. Therefore, after acquiring the three-dimensional data for a fracture surface, the point cloud data were processed by using the self-contained software CloudForm to reduce noise, remove irrelevant points, filter data ripples, and patch. Additionally, the original point cloud of the fracture surface is composed of hundreds of thousands of discrete points with an average spacing of 0.025 mm, which amount to a huge amount of disordered data. Data analysis was facilitated by compiling the program with MATLAB software to delete and reorder the measured points. The newly obtained fracture surface has a total of 22801 points at an average spacing of 1 mm. The measured topographic parameters of the fracture surfaces shown in Figure2 are listed in Table 1.Figure 2
Surfaces topographies of fractures (lower surfaces).
(a)
F1
(b)
F2
(c)
F3
(d)
F4
(e)
F5Table 1
Geometrical details of fractures.
Specimen
LengthL (mm)
Widthw (mm)
Peak asperity heightξ (mm)
VarianceRrms (mm2)
AverageRm (mm)
JRC
F1
149.9
150.1
3.74
2.06
1.58
11.3
F2
150.1
150.0
3.65
2.07
1.55
11.0
F3
150.1
150.1
3.97
4.14
3.52
12.8
F4
150.2
149.9
4.89
5.16
4.39
16.9
F5
150.1
149.9
3.27
3.87
3.29
13.9


## 3.1. Rock Fracture Preparation

The effect of fracture surface roughness on fluid flow was studied by way of saturated seepage tests of rock fracture surfaces of different roughness. A natural fracture surface is difficult to obtain, so artificial tension fracture specimens were used to study the characteristics of fluid flow in fractures. In this study, natural granite selected from a quarry in Sichuan Province was processed into150mm×150mm×150mm square specimens in the laboratory, and then, the specimens were split using the Brazilian splitting test method to obtain artificial tensional joint specimens. Finally, five groups of fracture surfaces (F1, F2, F3, F4, and F5) with different roughness were prepared.


## 3.2. Measurement Procedure

A portable 3D optical three-dimensional scanning system was used to measure the three-dimensional morphology of the fracture surface (Figure1). The system broadly comprises a scanning control unit, scanning lens, turn table, and tripod. The scanning lens is placed on the tripod, which can rotate freely and adjust its position conveniently. The other components are connected via USB. The system acquires a visible grating image projected onto the surface of the object then accurately determines the spatial coordinates (X, Y, Z) of each point using the phase and triangulation methods according to the shape of fringe curvature change, forming a three-dimensional point cloud data. This approach has the benefits of being fast, high precision (measuring accuracy 25 μm), and allowing noncontact measurement.Figure 1
Stereo topometric scanning system.In the actual measurement process, features of the measurement environment (light, dust, etc.) and the measurement methods will have an impact on the accuracy of the three-dimensional topographic data. Therefore, after acquiring the three-dimensional data for a fracture surface, the point cloud data were processed by using the self-contained software CloudForm to reduce noise, remove irrelevant points, filter data ripples, and patch. Additionally, the original point cloud of the fracture surface is composed of hundreds of thousands of discrete points with an average spacing of 0.025 mm, which amount to a huge amount of disordered data. Data analysis was facilitated by compiling the program with MATLAB software to delete and reorder the measured points. The newly obtained fracture surface has a total of 22801 points at an average spacing of 1 mm. The measured topographic parameters of the fracture surfaces shown in Figure2 are listed in Table 1.Figure 2
Surfaces topographies of fractures (lower surfaces).
(a)
F1
(b)
F2
(c)
F3
(d)
F4
(e)
F5Table 1
Geometrical details of fractures.
Specimen
LengthL (mm)
Widthw (mm)
Peak asperity heightξ (mm)
VarianceRrms (mm2)
AverageRm (mm)
JRC
F1
149.9
150.1
3.74
2.06
1.58
11.3
F2
150.1
150.0
3.65
2.07
1.55
11.0
F3
150.1
150.1
3.97
4.14
3.52
12.8
F4
150.2
149.9
4.89
5.16
4.39
16.9
F5
150.1
149.9
3.27
3.87
3.29
13.9


## 4. Calculation of Fractal Dimension

There are many methods for calculating the fractal dimension of a rough fracture surface. Clarke [30] first proposed the triangular prism surface area method, which takes the spatial surface area as the variable and a square grid as the measure scale. Later, Xie and Wang [23] proposed a projection coverage method. These two methods can be categorized as driver methods. The dimension calculated by these two methods is a similar fractal dimension, not a geometric fractal dimension. Zhou et al. [19] proposed box-counting methods for calculating the fractal dimension of a three-dimensional surface, including a cube coverage method and an improved cube coverage method. The above computational theories are based on statistical self-similarity. However, Brown [31] and Kulatilake et al. [32, 33] argue that rough rock fracture surfaces conform to the characteristics of a self-affine model.Kulatilake et al. [32] put forward a variogram method as a self-affine model to determine the fractal dimension. This takes the variogram function 2γx,h of the profile as a variable and the interval distance h as the measure scale.The detailed method is as follows:Step 1.
Generation of two-dimensional contour data in different directions. Firstly, the fracture surface data are meshed into1mm×1mm grids. 2D profiles are divided from the fracture surface data in accordance with the directions θ=15×kk=1,2,⋯,11. The height data Zx,y of a directional line is then calculated. The height data ZQ of the coordinate Qx,y is found on the 15-degree directional line P. The height of the coordinate radius within a 1 mm range is calculated according to equation (17). The profile data are then obtained cyclically. The next set of profile data with a distance of 10 mm is obtained by the same method. Finally, the profile data in each direction are obtained by repeated cyclic calculation.
(17)Zx,y=∑i=1nZi/ri∑i=1n1/ri,where Zi is the height of the point within a radius of 1 mm from point Q and ri is the distance from point i to point Q.Step 2.
Calculation of the fractal dimension of all of the directional lines. The variogram function is defined as(18)2γh=1N∑i=1NZi+1−Zi2,where γhL2 is the semivariogram, ZiL and Zi+1L are the heights of the 2D profile from the baseline, and N is the number of pairs of Z at a lag distance h between them. γh can be simplified as a power-law function in the self-affine profile as h→0:
(19)2γh=Kvh2H,where Kv is a proportionality constant and H is the Hurst exponent, which is related to the fractal dimension by Dv=2−H. However, equation (18) and equation (19) cannot be used to calculate Dv directly. Dv should be written in the logarithmic form
(20)log2γh=22−Dvlogh+logKv,so that Dv can be obtained by linear regression analysis. At least seven variance functions at different intervals h are calculated for each profile line, and fractal dimension D is obtained by fitting equation (20). The fractal dimension D of all of the profiles in one direction is averaged, and the fractal dimension D in that direction is obtained.To make the anisotropic characteristics of fracture surface roughness more intuitive, Table2 presents rose diagrams of the fractal dimensions of five different fracture surfaces (F1–F5). The fractal dimension is randomly distributed in all directions, and the fracture surfaces are characterized by an irregular anisotropic roughness structure. The fractal dimension of the fracture surface of F4 is larger than that of the others and is the highest, 1.60, in the 90° direction. Therefore, F4 has the greatest roughness of the fracture surfaces.Table 2
The fractal dimensions at different directions.
Directions (°)
F_1
F_2
F_3
F_4
F_5
0
1.422
1.497
1.421
1.418
1.475
30
1.398
1.486
1.391
1.426
1.426
60
1.431
1.477
1.417
1.486
1.435
90
1.542
1.503
1.434
1.600
1.436
120
1.478
1.466
1.471
1.489
1.389
150
1.425
1.459
1.370
1.424
1.407
180
1.422
1.497
1.421
1.418
1.475
210
1.398
1.486
1.391
1.426
1.426
240
1.431
1.477
1.417
1.486
1.435
270
1.542
1.503
1.434
1.600
1.436
300
1.478
1.466
1.471
1.489
1.389
330
1.425
1.459
1.370
1.424
1.407It is noted that the fractal dimensionD does not take the difference between forward and backward on the 2D profile into consideration. For example, the fractal dimension of the 0° profile is the same as that of 180°. Hence, equation (20) is correct when assuming there is the same nonlinear flow in the θ and θ+180° directions.


## 5. Nonlinear Flow Behavior of Rough-Walled Rock Fractures

### 5.1. Seepage Tests of Rough-Walled Rock Fractures

The five groups of rock fracture surfaces (F1–F5) mentioned in the previous section were prepared for flow testing with a self-designed device. A detailed description of the device is given in Xiong et al. [7]. Two flow directions (0° and 90°) were tested for each group of fractures, as shown in Figure 3. The pressure difference P under different flow rates Q was recorded during the test process; flow rates were in the range 0–100 ml/s during the tests.Figure 3
Different directions of flow (0° and 90°) in a rock fracture.


### 5.2. Nonlinear Flow Behavior

Figures4 and 5 show the relationship between the flow rate Q and the pressure gradient ∇P of each fracture in the 0° and 90° directions. When the flow rate is small (Q<10ml/s), the pressure gradient increases linearly with the flow rate. When the flow rate increases, the pressure gradient increases nonlinearly, showing an increase in the inertia effect. In order to describe this relationship, the Forchheimer formula was used to fit the test data for the 0° and 90° directions. The fitting results are shown in Table 3. These results indicate that Forchheimer’s law (equation (1)) is able to quantitatively describe the nonlinear flow behavior, which is consistent with Zimmerman et al. [17].Figure 4
Relationship between pressure gradient and flowrate in the 0° direction.Figure 5
Relationship between pressure gradient and flowrate in the 90° direction.Table 3
Fractal dimension and fitted parameters of Forchheimer’s law.
Specimen
A×106 (Kg·s-1·m-5)
B×1010 (Kg·m-8)
eh (mm)
Rec
D
bD
β
F1
0°
17.62
31.47
1.656
33.90
1.4221
0.00322
0.0027
90°
15.97
20.50
1.711
56.18
1.542
0.0232
0.002
F2
0°
11.04
15.01
1.935
43.53
1.497
0.0245
0.0021
90°
8.78
14.21
2.089
42.97
1.503
0.0292
0.0023
F3
0°
27.56
63.18
1.427
35.88
1.421
0.0414
0.0033
90°
27.47
56.93
1.428
38.03
1.434
0.0374
0.003
F4
0°
6.30
15.10
2.333
28.99
1.418
0.0433
0.0035
90°
5.09
6.09
2.504
57.82
1.600
0.0216
0.0018
F5
0°
13.07
20.08
1.829
40.84
1.475
0.0277
0.0022
90°
18.96
2.24
1.890
34.12
1.436
0.0034
0.00026In order to analyze the anisotropic characteristics of flow in a rock fracture, a new parameter describing anisotropy is proposed: the ratio of the difference between the 90° and 0° pressure gradients and the 90° pressure gradient.(21)anisotropy=∇P90°−∇P0°∇P90°,where ∇P90° is the pressure gradient in the 90° direction and ∇P0° is the pressure gradient in the 0° direction. Figure 6 shows the variation of anisotropy with flow rate. The anisotropy values differ between the different groups, which indicate that the anisotropy of fracture flow exists and is related to the fracture morphology and aperture distribution.Figure 6
Pressure anisotropy of fracture samples F1–F5.Normalized transmissivity (Ta/T0) is determined by Zimmerman et al. [17] in the following form:
(22)TaT0=11+βRe,where Ta is the apparent transmissivity and T0 is a special apparent transmissivity in Darcy’s flow state and is typically called intrinsic transmissivity. According to experimental data, the values of β are listed in Table 2. The relationship of bD and β is plotted in Figure 7. It can be seen from the figure that as bD increases, β linearly increases. And β is about 12 times than bD, which is consistent with Zimmerman et al. [17].Figure 7
The relationships ofbD and β of fracture samples F1–F5.


### 5.3. Verification of the Nonlinear Flow Model Based on Fractal Theory

In order to solve the undetermined constantsa and b, fracture morphology data were first obtained for F1, F2, and F3, and the fractal dimension was then calculated according to the method detailed in Section 4. Then, the seepage test results were fitted according to the Levenberg-Marquardt method, and the parameters a and b were found to be 0.246 and -0.964, respectively.In order to verify the model, the nonlinear fractal model is compared with the seepage test data and Chen’s two-parameter model [6].For fractures F4 and F5, pressure gradients were calculated according to the proposed model and the Chen model, respectively, and the results were compared with the experimental values, as shown in Figure8. It can be seen that the results calculated with the nonlinear fractal model are close to the measured values, and the relative errors are mostly within 20%. This shows that the nonlinear fractal model gives a better description of nonlinear seepage in fractured media than does Chen’s model. The error is larger in Chen’s model because absolute roughness is not sufficient to quantify the effect of surface topography on fracture flow.Figure 8
Comparison of results from equation (21), equation (22), and experimental measurements (open point is flow test along the 0° direction and solid point is flow test along the 90° direction; dashed line and solid line are proposed and the Chen model in 0° and 90° directions, respectively).
(a)
F4
(b)
F5


## 5.1. Seepage Tests of Rough-Walled Rock Fractures

The five groups of rock fracture surfaces (F1–F5) mentioned in the previous section were prepared for flow testing with a self-designed device. A detailed description of the device is given in Xiong et al. [7]. Two flow directions (0° and 90°) were tested for each group of fractures, as shown in Figure 3. The pressure difference P under different flow rates Q was recorded during the test process; flow rates were in the range 0–100 ml/s during the tests.Figure 3
Different directions of flow (0° and 90°) in a rock fracture.


## 5.2. Nonlinear Flow Behavior

Figures4 and 5 show the relationship between the flow rate Q and the pressure gradient ∇P of each fracture in the 0° and 90° directions. When the flow rate is small (Q<10ml/s), the pressure gradient increases linearly with the flow rate. When the flow rate increases, the pressure gradient increases nonlinearly, showing an increase in the inertia effect. In order to describe this relationship, the Forchheimer formula was used to fit the test data for the 0° and 90° directions. The fitting results are shown in Table 3. These results indicate that Forchheimer’s law (equation (1)) is able to quantitatively describe the nonlinear flow behavior, which is consistent with Zimmerman et al. [17].Figure 4
Relationship between pressure gradient and flowrate in the 0° direction.Figure 5
Relationship between pressure gradient and flowrate in the 90° direction.Table 3
Fractal dimension and fitted parameters of Forchheimer’s law.
Specimen
A×106 (Kg·s-1·m-5)
B×1010 (Kg·m-8)
eh (mm)
Rec
D
bD
β
F1
0°
17.62
31.47
1.656
33.90
1.4221
0.00322
0.0027
90°
15.97
20.50
1.711
56.18
1.542
0.0232
0.002
F2
0°
11.04
15.01
1.935
43.53
1.497
0.0245
0.0021
90°
8.78
14.21
2.089
42.97
1.503
0.0292
0.0023
F3
0°
27.56
63.18
1.427
35.88
1.421
0.0414
0.0033
90°
27.47
56.93
1.428
38.03
1.434
0.0374
0.003
F4
0°
6.30
15.10
2.333
28.99
1.418
0.0433
0.0035
90°
5.09
6.09
2.504
57.82
1.600
0.0216
0.0018
F5
0°
13.07
20.08
1.829
40.84
1.475
0.0277
0.0022
90°
18.96
2.24
1.890
34.12
1.436
0.0034
0.00026In order to analyze the anisotropic characteristics of flow in a rock fracture, a new parameter describing anisotropy is proposed: the ratio of the difference between the 90° and 0° pressure gradients and the 90° pressure gradient.(21)anisotropy=∇P90°−∇P0°∇P90°,where ∇P90° is the pressure gradient in the 90° direction and ∇P0° is the pressure gradient in the 0° direction. Figure 6 shows the variation of anisotropy with flow rate. The anisotropy values differ between the different groups, which indicate that the anisotropy of fracture flow exists and is related to the fracture morphology and aperture distribution.Figure 6
Pressure anisotropy of fracture samples F1–F5.Normalized transmissivity (Ta/T0) is determined by Zimmerman et al. [17] in the following form:
(22)TaT0=11+βRe,where Ta is the apparent transmissivity and T0 is a special apparent transmissivity in Darcy’s flow state and is typically called intrinsic transmissivity. According to experimental data, the values of β are listed in Table 2. The relationship of bD and β is plotted in Figure 7. It can be seen from the figure that as bD increases, β linearly increases. And β is about 12 times than bD, which is consistent with Zimmerman et al. [17].Figure 7
The relationships ofbD and β of fracture samples F1–F5.


## 5.3. Verification of the Nonlinear Flow Model Based on Fractal Theory

In order to solve the undetermined constantsa and b, fracture morphology data were first obtained for F1, F2, and F3, and the fractal dimension was then calculated according to the method detailed in Section 4. Then, the seepage test results were fitted according to the Levenberg-Marquardt method, and the parameters a and b were found to be 0.246 and -0.964, respectively.In order to verify the model, the nonlinear fractal model is compared with the seepage test data and Chen’s two-parameter model [6].For fractures F4 and F5, pressure gradients were calculated according to the proposed model and the Chen model, respectively, and the results were compared with the experimental values, as shown in Figure8. It can be seen that the results calculated with the nonlinear fractal model are close to the measured values, and the relative errors are mostly within 20%. This shows that the nonlinear fractal model gives a better description of nonlinear seepage in fractured media than does Chen’s model. The error is larger in Chen’s model because absolute roughness is not sufficient to quantify the effect of surface topography on fracture flow.Figure 8
Comparison of results from equation (21), equation (22), and experimental measurements (open point is flow test along the 0° direction and solid point is flow test along the 90° direction; dashed line and solid line are proposed and the Chen model in 0° and 90° directions, respectively).
(a)
F4
(b)
F5


## 6. Discussion

Forchheimer’s law has been widely applied for nonlinear seepage flow in fractured media, but the mechanism of the transition from linear flow to nonlinear flow needs to be further discussed.Zimmerman et al. (2014) considered the linear to nonlinear transition process in fracture fluid flow, distinguishing the two flow regimes by the nonlinear Darcy effect factorα with a value of 0.1:
(23)α=BQ2AQ+BQ2.Therefore, many scholars have used the critical Reynolds number to describe the transition mechanism. For fractured media, the Reynolds number can be expressed as follows:(24)Re=ρQμw.By substituting equation (16) and equation (23) into equation (24), a new critical Reynolds number can be obtained:
(25)Rec=5.42ehLc−0.964Dl−1.Equation (25) shows that the critical Reynolds number is closely related to hydraulic aperture, the fractal dimension of the fracture surface, and the flow direction. The smaller the hydraulic aperture and the rougher the fracture surface, the smaller the critical Reynolds number. The critical Reynolds numbers calculated by this method are shown in Table 2. Its value ranges between 30 and 60, much smaller than the 2300 value considered by Wang et al. (2015). The nonuniform distributions of the fracture aperture and the rough surface result in eddy currents and reflux flow, which make the flow tortuous and increase the inertial resistance. This leads to nonlinear flow at a low Reynolds number.


## 7. Conclusions

This paper discusses the effect of roughness on nonlinear flow in a rock fracture based on previous research and analysis of physical laboratory experiments. The main conclusions are as follows:(1)
A new nonlinear seepage model for rough fractures, equation (16), is proposed according to flow tortuosity in the fracture and the fractal characteristics of the fracture
(2)
The 3D optical three-dimensional scanning system was used to acquire point cloud data from fracture surfaces. The self-affine fractal dimension calculation method proposed by Kulatilake et al. [32] was used to analyze the anisotropic characteristics of the roughness of the fracture surface
(3)
Five different kinds of rough fractures were tested in seepage experiments in the 0° and 90° directions. The results show that fracture flow conforms to Forchheimer’s law and has clear isotropic characteristics. The new model generates results that are closer to those from the experiment than does Chen’s two-parameter model
(4)
According to the new model, a new expression of the critical Reynolds number (equation (25)) for distinguishing Darcy flow from Forchheimer flow is proposed. It shows that the critical Reynolds number is closely related to hydraulic aperture, the roughness of the fracture surface, and the flow direction



---

*Source: 5874849-2019-06-24.xml*