Abstract

The mechanism and quantitative descriptions of nonlinear fluid flow through rock fractures are difficult issues of high concern in underground engineering fields. In order to study the effects of fracture geometry and loading conditions on nonlinear flow properties and normalized transmissivity through fracture networks, stress-dependent fluid flow tests were conducted on real rock fracture networks with different number of intersections (1, 4, 7, and 12) and subjected to various applied boundary loads (7, 14, 21, 28, and 35 kN). For all cases, the inlet hydraulic pressures ranged from 0 to 0.6 MPa. The test results show that Forchheimer’s law provides an excellent description of the nonlinear fluid flow in fracture networks. The linear coefficient and nonlinear coefficient in Forchheimer’s law generally decrease with the number of intersections but increase with the boundary load. The relationships between and can be well fitted with a power function. A nonlinear effect factor was used to quantitatively characterize the nonlinear behaviors of fluid flow through fracture networks. By defining a critical value of = 10%, the critical hydraulic gradient was calculated. The critical hydraulic gradient decreases with the number of intersections due to richer flowing paths but increases with the boundary load due to fracture closure. The transmissivity of fracture networks decreases with the hydraulic gradient, and the variation process can be estimated using an exponential function. A mathematical expression for decreased normalized transmissivity against the hydraulic gradient was established. When the hydraulic gradient is small, holds a constant value of 1.0. With increasing hydraulic gradient, the reduction rate of first increases and then decreases. The equivalent permeability of fracture networks decreases with the applied boundary load, and permeability changes at low load levels are more sensitive.

1. Introduction

Rock fracture networks constitute the main pathways of fluid flow and solute migration in deep underground projects, and during the past several decades, substantial efforts have been devoted to the estimation of fluid flow behavior and transmissivity of fractures in many geoengineering and geosciences such as underground tunneling [13], CO2 sequestration [4, 5], geothermal energy extraction [68], and hazardous wastes isolation [911]. The fluid flow in rock fractures is commonly assumed to follow the cubic law, in which the flow rate is linearly proportional to the pressure gradient [1214]. However, when the flow rate/hydraulic head difference is large, deviation from the linear Darcy law may occur. In such case, the conductivity of the fractures calculated using the cubic law will be overestimated [1517]. In addition, natural rock fractures are often characteristic of rough walls, intersections, and asperity contacts, which make the fluid flow process even more complex and difficult to accurately describe [1822]. Therefore, a thorough understanding of nonlinear flow properties within fracture networks is of great significance to ensure performance and safety of engineering activities.

The transmissivity , which is linked with volumetric flow rate, pressure gradient, and the fracture width, has been applied to characterize the onset of nonlinear flow regimes in rock fractures in some previous studies. Zimmerman et al. [23] conducted Navier-Stokes simulations and fluid flow tests in a natural 3D sandstone fracture and confirmed that weak inertia regimes exist for Reynolds number (Re) in the range of 1–10, where the fluid flow enters the nonlinear flow regime. The normalized transmissivity () remains constant when but shows a decreasing trend when (see Figure 1(a)). Zhang and Nemcik [24] experimentally investigated the fluid flow regimes through deformable rock fractures subjected to changing confining stress from 1.0 to 3.5 MPa. They found that for a nonmated sandstone fracture with the joint roughness coefficient of 9.2, the transmissivity is not a constant value but decreases with the hydraulic gradient (see Figure 1(b)). Liu et al. [14] calculated the - curves of single fractures with various ratios of mechanical aperture to the fracture length, in which the mechanical aperture is defined as the arithmetic average of the point-to-point distance between the two walls of a fracture. The results show that as increases, decreases. When , remains constant approaching 1.0; the fluid flow is linear. For , decreases dramatically; fluid flow enters a transitional regime. When , is sufficiently small and approaches 0; fluid flow enters a complete nonlinear regime (see Figure 1(c)). Yin et al. [25] estimated the relationships between and Re for rough-walled fractures during shearing and subjected to various normal loads from 7 to 35 kN. As an example, for one of their samples with the shear displacement of 9 mm, exhibits a decreasing trend with Re and the rate of its decrease diminishes. as a function of Re can be well described using a six-order polynomial function (see Figure 1(d)). The above analysis generally focused on variations in (normalized) transmissivity of single fractures. However, when it comes to complex fracture networks, the stress-dependent nonlinear flow mechanisms and transmissivity are still not fully understood.

In engineering practices, both nature activities and human perturbations entail significant changes in the effective stress field of underground rock masses, which would then impact the fracture aperture in fracture networks. The void space between opposing surfaces can vary due to normal stress-induced closures or opening [18, 24, 26] or shear stress-induced dilations [2729]. Hence, the fluid flow behavior and the transmissivity of rock fracture networks are stress-dependent, which has so far rarely been investigated.

The purpose of this paper is to investigate the stress-dependent nonlinear flow properties and normalized transmissivity of real rock fracture networks. First, plate granite specimens containing fracture networks with different number of intersections (1, 4, 7, and 12) were machined using a high-pressure water jet cutting system. Next, a series of high precision stress-dependent water flow tests with respect to different inlet hydraulic pressures ranging from 0 to 0.6 MPa and increasing applied boundary loads from 7 to 35 kN were conducted. The nonlinear flow behaviors, variations of the flow nonlinearity, normalized transmissivity and equivalent permeability of the fracture networks as a function of fracture networks geometry, and boundary load conditions were all examined.

2. Governing Equations of Fluid Flow in Fractures

Fluid flow through a single fracture is generally governed by the following Navier-Stokes (NS) equations, written in a tensor form as [15, 30, 31]:where is the flow velocity vector, is the hydraulic pressure, is the shear stress tensor, is the fluid density, is the time, is the body force acting on the fluid, and in gravity environment denotes the gravitational acceleration term.

For the case of incompressible and steady-state Newtonian flow, the terms involving time drop out and the NS equations can be expressed as [1, 17]where is the dynamic viscosity.

Equation (2) is composed of a set of coupled nonlinear partial derivatives of varying orders. In order to solve these equations, an additional equation, called the mass conservation equation which can be achieved through mass continuity, is employed and the equation takes the form [1, 15]

In (2), the convective acceleration term denotes the inertial forces acting on the fluid, which are the source of nonlinearity [14]. The complexity of the Navier-Stokes equations makes exact solutions extremely difficult to achieve, especially for natural rock fractures with complicated geometries.

For certain cases with very low Reynolds number or flow rate, by assuming that the inertial forces of fluid flow through the fractures are negligible compared to the viscous forces. The convective acceleration term in (2) vanishes, and the Navier-Stokes equations can be reduced to solvable forms as [12]where is the total volumetric flow rate, is the pressure gradient, is the fracture width, and is the hydraulic aperture. The flow rate is proportional to the cube of . This equation is commonly referred to as the well-known cubic law.

The linear relation between and in (4) can only be anticipated for parallel-plate models under a limited range of flow rates. Besides, natural rock fractures are characterized by fracture intersections and complex fracture surface geometries, and the inertia items are often not zero. With increasing hydraulic head differences, the flow rate is nonlinearly related to the pressure drop and (4) is no longer applicable.

Some empirical expressions were proposed to describe the nonlinear flow in fractures, and Forchheimer’s law is the most extensively used approach where the pressure gradient is a quadratic function of the flow rate, written as [3234]where and are model coefficients, representing pressure drop components caused by linear and nonlinear effects, respectively.

Since the hydraulic gradient is proportional to as , (5) can then be written as follows:where and .

The hydraulic gradient is calculated by dividing the hydraulic head difference by the flow lengthwhere is the flow length, which denotes the horizontal distance between the left and right specimen boundary in this study, and is the gravitational acceleration.

To quantitatively estimate the nonlinearity of fluids flowing through the fracture networks, a nonlinear effect factor was introduced to determine the fluid flow regime [1, 35, 36]. It can be presented bywhere and are energy losses due to viscous and inertial dissipation mechanisms in the fractures. denotes the ratio of nonlinear pressure drop to the overall pressure drop.

For engineering purposes, a critical value of = 10% has been defined as the critical condition for onset of flow nonlinearity in fracture networks, where the nonlinear effect can be appreciable and cannot be neglected [23, 24, 29].

The Forchheimer number is another widely accepted parameter to characterize the onset of flow transition to nonlinearity, which is defined as the ratio of nonlinear to linear pressure drops [1, 13, 17],

Combinations of (8) and (9) yield the following equation:

Transmissivity () is an important hydraulic property to characterize the macroscopically observed flow resistances in fractures, which can be defined as [37]

When the flow rate is sufficiently low and the inertial forces are negligible, the intrinsic transmissivity () is commonly regarded as a constant value. As the hydraulic head difference increases, the transmissivity can be used to evaluate the nonlinear flow through the fractures, and the normalized transmissivity () was then determined [23],

Therefore, as the nonlinear term () contributes to 10% of the total pressure drop or else = 10%, the critical value of is equal to 0.9. The corresponding is defined as the critical hydraulic gradient .

3. Experimental Methodology

3.1. Rock Fracture Specimen Preparation

Granite specimens (495 × 495 × 17 mm in size) containing artificial fracture networks were established by using a high-pressure water jet cutting system [25]. The granites were taken from Linyi city in Shandong Province of China. The granite is a porphyritic monzonite with an average unit weight of about 2.69 g/cm3. The uniaxial compressive strength of the rock is about 97.54 MPa, and its permeability is in the order of magnitude of  m2.

The parallelism between the upper and lower surfaces of the plate specimens was controlled within an error of 0.02 mm. The water pressure of the water jet cutting system was held constant for all rock fractures, and the fractures penetrated the plate specimens thoroughly. Detailed descriptions of the rock fracture specimens are illustrated in Figure 2. Each specimen consists of two sets of parallel fractures with a constant included angle of 60°. By varying the spacing of the fractures, fracture networks with different numbers of intersections ( = 1, 4, 7, and 12) were, respectively, created.

Using a high-resolution three-dimensional laser scanning profilometer system, the fracture surface topography was measured before the hydromechanical tests. We considered the average of joint roughness coefficient (JRC) values for a series of 2D profiles along a fracture surface in the length direction, which is a suggested method by the International Society for Rock Mechanics and Rock Engineering (ISRM) to calculate the JRC of a 3D single fracture [38]. The spacing of 2D profiles is 3 mm. In total, the JRC values of 5 2D profiles are calculated and their average value is regarded as the JRC of a 3D single fracture. The JRC values of all 2D profiles were estimated based on the values of using (13) proposed by Tse and Cruden [39]. Here, JRC is a dimensionless index ranging from 1 to 20 [40], which has been widely accepted in the field of rock mechanics and rock engineering [22, 41].

The results indicate that JRC values of the fractures fluctuate within a very small range around 3.47. The fractures closely approximate parallel-plate models with a small JRC value.where and represent the coordinates of a fracture surface profile and is the number of sampling points along the length of a fracture. is the root mean square slope of the 2D profile.

3.2. Experimental Setup and Procedure

Before the hydromechanical tests, sealing of the fractured specimens was first conducted. Based on the model specimen size, a 3 mm thick rubber jacket was produced using the ethylene propylene diene monomer (EPDM) waterproof rubber [25]. Then, a rock fracture specimen was placed in the jacket, with the extra rubber seamlessly glued to the rock surface. For a better sealing of the fractures, a layer of glass glue was evenly coated on the specimen surface to leave the fractures open. Then, a piece of transparent crystal plate with a suitable size was attached to the glass glue. When the glass glue was dried, a hand-operated electric drill with a diameter of 10 mm was used to drill circular holes at the water inlet and outlet positions of the fractures through the rubber jacket. After a high-strength plexiglass cover plate and the horizontal load devices with uniform flow chambers were installed, the specimen was transferred to the test area.

The stress-dependent flow tests of the fracture networks were conducted by using a self-developed flow test apparatus [25], as shown in Figure 3. It mainly consists of three units: (i) a platform for fluid flow in fracture networks; (ii) a water supply and measurement unit; (iii) a pneumatic-hydraulic cylinder unit. During the flow tests, a constant water pressure was maintained at the flow inlet manifold using an air compressor connected to a water tank, with the pressure range of 0–2 MPa. Both water inlet directions (- and -directions) evenly featured 12 flow distribution chambers to achieve a variable but uniform flow field. The hydraulic sources of the water inlet or outlet manifold can be individually switched on or off using the shut-off valves (Figure 4(a)). Each of the 12 water chambers provided equal water pressure at the rock specimen boundary. The effluents flowing out of the fractures were measured in real time using four glass tube rotameters with a range of 0.0004–11.0 L/min. Both the horizontal (- and -directions) and vertical (-direction) loads on the rock specimens were created using a pneumatic-hydraulic cylinder connected to an air compressor that has a maximum air pressure of 3 MPa. Before the test, a vertical load of 20 kN was applied on the specimen surface to balance the vertical water pressure in the fractures.

During the flow test, water was fed through the inflow manifold connected to a water tank that can supply a constant water head at all times by using an air compressor. Both horizontal boundary loads in the -direction () and -direction () were applied and increased from 7 to 35 kN in a 7 kN interval, with a fixed lateral pressure coefficient of to of 1.0. Applications of the boundary load conditions are displayed in Figure 4(b). In this study, we just discuss the -directional nonlinear flow behaviors through the fracture networks over a range of inlet water pressures from 0 to 0.6 MPa with increasing boundary loads. To achieve this investigation, the side valves (-direction) were opened while the top and bottom valves (-direction) were closed, which forced water to flow horizontally from the right to left specimen boundary (-direction) through the fracture networks. Other boundaries were considered impermeable.

For a certain boundary load and inlet hydraulic pressure, the total volume flow rate of the fracture networks at the water outlet boundary can be obtained when the glass rotor was relatively stable with no fluctuations. The effluent was then collected in another storage tank and recirculated. The entire hydraulic experiments were performed under isothermal conditions at room temperature of approximately 20°C. In addition, the fluid was assumed to be viscous with of 1.0 × 10−3 Pa·s and incompressible with of 1000 kg/m3.

4. Test Results and Discussion

Due to the low permeability of the intact granite matrix, fluid was assumed to flow through the fractures only during the hydromechanical tests. For the rock fracture specimen with various subjected to different (here, refers to the boundary load due to = ), a series of hydraulic tests were conducted with different inlet fluid pressures ranged from 0 to 0.6 MPa. In this study, the hydraulic head at the water outlet specimen boundary was assumed to be equal to zero. Thus, according to (7), in the range of 0–0.6 MPa, ranged from 0 to 123.69.

The experimental data obtained in the hydraulic tests on the rock fracture specimens in the form of hydraulic gradient () versus discharge () curves are displayed in Figure 5. The quadratic polynomial regression of Forchheimer (6) is used to best fit the raw experimental data, with the residual squared values all larger than 0.99 (Table 1). From Figure 5, with an increase in , shows an increasing trend. For a certain , as increases, the slope of - curves becomes steeper, indicating a higher flow resistance mainly as a result of fracture closure. Hence, more flow energy is required to achieve the same flow rate for fracture networks subjected to a larger boundary load. However, at a given , as increases, flowing out of the fracture networks produced by an identical hydraulic head shows an increasing trend (Figure 6). Taking = 14 kN as an example, produced by = 123.69 is increased from 5.04 × 10−6 m3/s () to 10.86 × 10−6 m3/s (), increasing by approximately a factor of 1.15. With an increase in , the overall discharge capacity of the fracture networks enhances.

The regression linear and nonlinear coefficients and in Forchheimer (6) for all test cases were then calculated and tabulated in Table 1. Figures 7(a)7(d) show the variations of the coefficients and in terms of and . For a certain , as increases, both and show an increasing trend. In the range of from 7 to 35 kN, increases by 7.10–13.99 times and increases by 3.97–7.73 times. The increase of coefficient value is due to fracture closure induced by the boundary loads. However, for a certain , as increases, both and decrease. The variation process can be divided into two stages. For = 1–7, the coefficients decrease significantly as varies. However, in the range of 7–12, the coefficients vary gradually and approach essentially constant values. Using coefficient as an example, as increases from 1 to 12, decreases by 71.91% ( = 7 kN), 63.67% ( = 14 kN), 57.51% ( = 21 kN), 74.47% ( = 28 kN), and 72.90% ( = 35 kN), respectively.

In view of the fact that the variations of linear and nonlinear coefficients and against and are very similar, as a function of was analyzed [29], as plotted in Figure 7(e). The experimental data can be fitted very well using the following empirical function:

Although (14) fits very well the relationship between the coefficients and , with the generating of 0.9927, the applicability of this equation needs to be further verified.

To quantify the degree of the nonlinear effect of fluid flow through the fracture networks, variations of for all test cases were calculated and plotted as a function of based on (8), as shown in Figure 8. As increases, displays an increasing trend, while the rate of its increase steadily diminishes. The variation process can be well described using a power functionwhere the fitting coefficients ( and ) are related to both the boundary load conditions and number of intersections. For a certain , as increases, the coefficient shows a decreasing trend, but the coefficient increases. However, at a given , with an increase in , plots an increasing trend while decreases.

By defining a critical value of = 0.1, of the fracture networks with various numbers of intersections was calculated, as listed in Table 1. Figure 9 displays the changes in due to the changes in and . As increases, shows an increasing trend. In the range of from 7 to 35 kN, increases by 6.35 times (), 7.73 times (), 3.97 times (), and 6.09 times (), respectively. The main reason for these variations may be as follows. For = , it is expected that the effect of normal loads plays a dominant role on the overall fluid flow field [42]. An increase in the normal load will undoubtedly lead to corresponding fracture closure. For each single fracture, a small change in fracture aperture can result in a large variation in the flow rate due to proportional relationships between and in (4), which would then impact for the onset of flow transition. However, for a certain , as increases, shows a decreasing trend. Taking = 35 kN as an example, in the range of from 1 to 12, decreases from 118.57 to 13.01, or by 89.03%.

For fluid flow through fractured and porous media, the transmissivity has also been applied to estimate the nonlinear flow regime. Using (11), of the fracture networks was calculated. Figure 10 shows the relationships between and . Notably, is not a constant value but exhibits a decreasing trend with , which further validates the existence of nonlinear flow behaviors in the fracture networks. With an increase in , decreases due to fracture closure. However, for a certain , shows an increasing trend with . Using the least square method, as a function of can be well fitted with the following exponential function:where the regression coefficients , , and were all presented in Figure 10. The units of coefficients and are m4 and is a dimensionless constant.

When assessing the fluid flow through a single rock fracture, the Reynolds number (Re), which is defined as the ratio of inertial forces to viscous forces, is typically used to quantify the onset of nonlinear flow. As noted above, the variations of normalized transmissivity () versus Re can be determined by extending the linear Darcy law with [23, 25, 34, 43],where is a dimensionless coefficient determined by the coefficients and in Forchheimer (6).

However, in engineering practices, there are hundreds to thousands of fractures, and the Re of each fracture generally cannot be ascertained. Instead, in a model typically has known values and can be easily obtained. Therefore, in this study, the relationships between and were used to evaluate the flow regimes in the fracture networks.

By fitting the experimental data sets, the relationships between and for rock fracture networks in this study were analyzed, as shown in Figure 11. Here, denotes the transmissivity corresponding to = 0 in Figure 10, where the hydraulic head difference is sufficiently low and the inertial forces are negligible. The results show that variations in against can be expressed as follows [14]:where is a dimensionless coefficient. Note that the exponent of is −0.45, which does not change with the number of intersections or boundary load conditions.

From Figure 11, as increases, shows a decreasing trend. For a certain , as increases, the transmissivity relationships generally shift upward. However, for a certain , the curves shift downward as increases. In addition, the variations in with can be divided into three stages. When is small (i.e., less than 0.2), with the increase of , approximately holds a constant value of 1.0; thus, the fluid flow is linear. Then, with continuously increasing , decreases and the reduction rate of first increases and then decreases. Based on (12) and (18), when = 0.9, is calculated, which is in the ranges 8.98–103.16 (), 3.46–61.40 (), 2.69–28.37 (), and 1.10–15.50 (), respectively. Generally, the range of is larger than that calculated based on (6) and (8) shown in Figure 9.

The variations of the coefficient with varying and are shown in Figure 12. For a certain , as increases, increases. In the range of 7–35 kN, generally increases by a factor of 2. However, at a given , the coefficient decreases with . As increases from 1 to 12, shows a decrease of 61.14% ( = 7 kN), 66.37% ( = 14 kN), 42.98% ( = 21 kN), 55.29% ( = 28 kN), and 57.38% ( = 35 kN), respectively.

From the quadratic relationships between and of the fracture networks displayed in Figure 5, at a given , the increase in varies with both and , which leads to the changes in the equivalent permeability of the fracture networks. To assess the discharge capacity of the fracture networks, the -directional equivalent permeability of the fracture networks was calculated at a specified inlet hydraulic pressure of 0.05 MPa, with the corresponding of 10.31. The fluid flow took place in fractures only and was calculated by the cubic law.

Calculations were made with the following equation, assuming no gravity term:where is the cross-sectional area.

The changes in the equivalent permeability of the fracture networks in terms of and are displayed in Figure 13. For a certain , as increases, of the model region decreases. In the range of from 7 to 35 kN, decreases by 90.31% (), 90.87% (), 82.28% (), and 84.82% (), respectively. In addition, the - variation curves can be divided into two stages. As = 7–21 kN, the permeability changes are sensitive. However, when > 21 kN, changes gradually and approaches constant values. The main reason for these variations may be as follows. As increases, the effective stress in the fractures increases, resulting in corresponding fracture closure. The overall volume flow rate decreases. As the fracture aperture continues to decrease until the residual fracture aperture, the equivalent permeability of the fracture networks generally retains a constant value. However, at a given F, the equivalent permeability increases with mainly due to richer connectivity in the fracture networks.

5. Conclusions

This paper experimentally examines the impacts of number of intersections ( = 1, 4, 7, and 12) on stress-dependent nonlinearity of fluid flow through real rock fracture networks. For each intersection number, a series of hydromechanical tests with respect to different inlet hydraulic pressures (0–0.6 MPa) and increasing boundary loads from 7 to 35 kN were conducted. The degree of nonlinearity, critical hydraulic gradient, normalized transmissivity, and the equivalent permeability of the fracture networks were all evaluated. The main conclusions are as follows:(1)Forchheimer’s law offers a good description for the relationships between volume flow rate and hydraulic gradient of water flow through the rock fracture networks. Both the linear coefficient and nonlinear coefficient in Forchheimer’s law decrease with the number of intersections but increase with the boundary load. The changes of nonlinear coefficient as a function of linear coefficient can be well fitted using a power function based on the experimental data.(2)The critical hydraulic gradient of the fracture networks is calculated by taking a critical nonlinear effect factor value of 10%. With an increase in the number of intersections, the critical hydraulic gradient decreases mainly due to richer flowing pathways in the fracture networks. However, as the boundary load increases, the critical hydraulic gradient increases due to fracture closure caused by increasing effective stress in the fractures. For all cases in this study, the critical hydraulic gradient is ranged between 0.62 and 118.57.(3)With an increase in the hydraulic gradient, the transmissivity of the fracture networks displays an exponential decrease trend. In addition, the transmissivity increases with the number of intersections but decreases with the applied boundary load. The variations in normalized transmissivity as a function of hydraulic gradient were estimated with a mathematical expression . The coefficient increases with the boundary load but decreases with the number of intersections. As the boundary load increases, the equivalent permeability of the fracture networks shows a decreasing trend.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

The financial support from the Fundamental Research Funds for the Central Universities, China (no. 2018QNA32), is gratefully acknowledged.