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.

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 [15]. 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 [611] 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 [1216]. 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: where is the pressure gradient between the inlet and outlet of the fracture, is the flow rate through the fracture, and and are the coefficients of viscosity and inertia, respectively. Zimmerman et al. [17] observed the Forchheimer flow phenomena of rough fractures when the Reynolds number 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 and 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 part and a nonlinear part . When the flow rate is low, a cubic law can be used to describe the relationship of the flow rate and pressure. Hence, can be expressed as where is the dynamic coefficient of viscosity of the fluid and is the width of the fracture. is the hydraulic aperture, which can be calculated as .

The coefficient represents the degree that the seepage curve deviates from that in the linear stage. Schrauf and Evans [28] put forward a form of using dimensional analysis: where is the fluid density and 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 : where and 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): where , , and are fitting parameters. is the hydrological curvature, which is defined as the ratio of the actual length of the seepage path to the horizontal length along the pressure gradient of the fracture. 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 measure and measurement scale is as follows: where is a parameter related to the fractal dimension and is in the range 1-3. is the apparent measurement. Based on this, the relationship between the fracture surface area and a square mesh with dimension is as follows: where is the fractal dimension of the fracture area. For a square fracture, when is equal to and equation (7) is substituted into , can be obtained as follows:

Hence, can be obtained:

In addition, the fractal relationship of the length of the seepage path and measure is as follows: where is the fractal dimension of the seepage path. When is equal to and equation (10) is substituted into , can be obtained as

Hence, can be obtained:

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:

Jin et al. [29] considered as equal to . Hence, can be obtained by substituting equation (9) and equation (12) into equation (5):

When and are , equation (1) becomes

Hence, a new model of parameter in Forchheimer’s law can be obtained: where and 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 and .

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 into 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 (Figure 1). 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 (, , ) 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.

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 Figure 2 are listed in Table 1.

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 of the profile as a variable and the interval distance 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 into grids. 2D profiles are divided from the fracture surface data in accordance with the directions . The height data of a directional line is then calculated. The height data of the coordinate is found on the 15-degree directional line . 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. where is the height of the point within a radius of 1 mm from point and is the distance from point to point .

Step 2. Calculation of the fractal dimension of all of the directional lines. The variogram function is defined as where is the semivariogram, and are the heights of the 2D profile from the baseline, and is the number of pairs of at a lag distance between them. can be simplified as a power-law function in the self-affine profile as : where is a proportionality constant and is the Hurst exponent, which is related to the fractal dimension by . However, equation (18) and equation (19) cannot be used to calculate directly. should be written in the logarithmic form so that can be obtained by linear regression analysis. At least seven variance functions at different intervals are calculated for each profile line, and fractal dimension is obtained by fitting equation (20). The fractal dimension of all of the profiles in one direction is averaged, and the fractal dimension in that direction is obtained.

To make the anisotropic characteristics of fracture surface roughness more intuitive, Table 2 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.

It is noted that the fractal dimension 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 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 under different flow rates was recorded during the test process; flow rates were in the range 0–100 ml/s during the tests.

5.2. Nonlinear Flow Behavior

Figures 4 and 5 show the relationship between the flow rate and the pressure gradient of each fracture in the 0° and 90° directions. When the flow rate is small (), 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].

In 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. where is the pressure gradient in the 90° direction and 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.

Normalized transmissivity () is determined by Zimmerman et al. [17] in the following form: where is the apparent transmissivity and 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 and is plotted in Figure 7. It can be seen from the figure that as increases, linearly increases. And is about 12 times than , which is consistent with Zimmerman et al. [17].

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

In order to solve the undetermined constants and , 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 and 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 Figure 8. 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.

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:

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:

By substituting equation (16) and equation (23) into equation (24), a new critical Reynolds number can be obtained:

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

Data Availability

The data are all available and have been explained in this article; readers can access the data supporting the conclusions of the study.

Conflicts of Interest

The authors declare that they have no conflicts of interest

Authors’ Contributions

Chun Zhu and Yun Lin contributed equally to this work.

Acknowledgments

This work was supported by the Key Research and Development Project of Zhejiang Province (Grant No. 2019C03104) and the Special Funds of Fundamental Research Funds for Central Universities (2015QB02). The first author is grateful to the Chinese Scholarship Council and the University of Adelaide for providing an opportunity to conduct this research as a joint Ph.D. student.