# Experimental Study on Stress-Dependent Nonlinear Flow Behavior and Normalized Transmissivity of Real Rock Fracture Networks

Authors: Qian Yin; Hongwen Jing; Richeng Liu; Guowei Ma; Liyuan Yu; Haijian Su
Journal: Geofluids
 (2018)

---

## 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 coefficienta and nonlinear coefficient b in Forchheimer’s law J=aQ+bQ2 generally decrease with the number of intersections but increase with the boundary load. The relationships between a and b can be well fitted with a power function. A nonlinear effect factor E=bQ2/(aQ+bQ2) was used to quantitatively characterize the nonlinear behaviors of fluid flow through fracture networks. By defining a critical value of E = 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 T/T0=1-exp⁡(-αJ-0.45) for decreased normalized transmissivity T/T0 against the hydraulic gradient J was established. When the hydraulic gradient is small, T/T0 holds a constant value of 1.0. With increasing hydraulic gradient, the reduction rate of T/T0 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.


---

## Body

## 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 [1–3], CO2 sequestration [4, 5], geothermal energy extraction [6–8], and hazardous wastes isolation [9–11]. 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 [12–14]. 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 [15–17]. 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 [18–22]. Therefore, a thorough understanding of nonlinear flow properties within fracture networks is of great significance to ensure performance and safety of engineering activities.The transmissivityT, 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 (T/T0) remains constant when Re<1 but shows a decreasing trend when Re=1–10 (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 J (see Figure 1(b)). Liu et al. [14] calculated the T/T0-J 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 J increases, T/T0 decreases. When J=10-7–10-5, T/T0 remains constant approaching 1.0; the fluid flow is linear. For 10-4<J<100, T/T0 decreases dramatically; fluid flow enters a transitional regime. When J>102, T/T0 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 T 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, T exhibits a decreasing trend with Re and the rate of its decrease diminishes. T 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.Figure 1
Relationships between hydraulic gradientJ (Reynolds number Re) and the (normalized) transmissivity in previous studies. (a) Relationships between T/T0 and Re (after [23]). (b) Relationships between T and J (after [24]). (c) Relationships between T/T0 and J (after [14]). (d) Relationships between T and Re (after [25]).
(a) (b) (c) (d)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 [27–29]. 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]:(1)ρ∂u∂t+u·∇u=-∇p+∇·T+ρf,where u is the flow velocity vector, p is the hydraulic pressure, T is the shear stress tensor, ρ is the fluid density, t is the time, f is the body force acting on the fluid, and in gravity environment f denotes the gravitational acceleration term.For the case of incompressible and steady-state Newtonian flow, the terms involving timet drop out and the NS equations can be expressed as [1, 17](2)ρu·∇u=-∇p+μ∇2u+f,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](3)∇·u=0.In (2), the convective acceleration term (u·∇)u 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(u·∇)u in (2) vanishes, and the Navier-Stokes equations can be reduced to solvable forms as [12](4)Q=-wbh312μ∇p,where Q is the total volumetric flow rate, ∇p is the pressure gradient, w is the fracture width, and bh is the hydraulic aperture. The flow rate is proportional to the cube of bh. This equation is commonly referred to as the well-known cubic law.The linear relation betweenQ and ∇p 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 [32–34](5)-∇P=a′Q+b′Q2,where a′ and b′ are model coefficients, representing pressure drop components caused by linear and nonlinear effects, respectively.Since the hydraulic gradientJ is proportional to ∇p as J=∇p/(ρg), (5) can then be written as follows:(6)J=aQ+bQ2,where a=-ρga′ and b=-ρgb′.The hydraulic gradientJ is calculated by dividing the hydraulic head difference by the flow length(7)J=PρgL,where L is the flow length, which denotes the horizontal distance between the left and right specimen boundary in this study, and g is the gravitational acceleration.To quantitatively estimate the nonlinearity of fluids flowing through the fracture networks, a nonlinear effect factorE was introduced to determine the fluid flow regime [1, 35, 36]. It can be presented by(8)E=bQ2aQ+bQ2,where aQ and bQ2 are energy losses due to viscous and inertial dissipation mechanisms in the fractures. E denotes the ratio of nonlinear pressure drop to the overall pressure drop.For engineering purposes, a critical value ofE = 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 numberF0 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],(9)F0=bQ2aQ=bQa.Combinations of (8) and (9) yield the following equation:(10)E=F01+F0.Transmissivity (T) is an important hydraulic property to characterize the macroscopically observed flow resistances in fractures, which can be defined as [37](11)-∇P=μQTw.When the flow rate is sufficiently low and the inertial forces are negligible, the intrinsic transmissivity (T0) 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 (T/T0) was then determined [23],(12)TT0=-μQ/waQ+bQ2-μQ/waQ=1-E.Therefore, as the nonlinear term (bQ2) contributes to 10% of the total pressure drop or else E = 10%, the critical value of T/T0 is equal to 0.9. The corresponding J is defined as the critical hydraulic gradient Jc.


## 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 10-19 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 Figure2. 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 (N = 1, 4, 7, and 12) were, respectively, created.Figure 2
Plate granite specimens of fracture networks with different numbers of intersections (N = 1, 4, 7, and 12). The black lines denote the geometric locations the fractures.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 Z2 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.(13)Z2=1M∑zi-1-zixi-1-xi21/2JRC=32.2+32.47log⁡Z2,where xi and zi represent the coordinates of a fracture surface profile and M is the number of sampling points along the length of a fracture. Z2 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 (x- and y-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 (x- and y-directions) and vertical (z-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.Figure 3
Schematic view of the stress-dependent fluid flow test apparatus.Figure 4
(a) A top view of hydraulic setup of the test apparatus. (b) Loading diagram of the rock fracture specimen.
(a) (b)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 thex-direction (Fx) and y-direction (Fy) were applied and increased from 7 to 35 kN in a 7 kN interval, with a fixed lateral pressure coefficient λ of Fy to Fx of 1.0. Applications of the boundary load conditions are displayed in Figure 4(b). In this study, we just discuss the x-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 (x-direction) were opened while the top and bottom valves (y-direction) were closed, which forced water to flow horizontally from the right to left specimen boundary (x-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−3Pa·s and incompressible with ρ of 1000 kg/m3.


## 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 10-19 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 Figure2. 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 (N = 1, 4, 7, and 12) were, respectively, created.Figure 2
Plate granite specimens of fracture networks with different numbers of intersections (N = 1, 4, 7, and 12). The black lines denote the geometric locations the fractures.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 Z2 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.(13)Z2=1M∑zi-1-zixi-1-xi21/2JRC=32.2+32.47log⁡Z2,where xi and zi represent the coordinates of a fracture surface profile and M is the number of sampling points along the length of a fracture. Z2 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 (x- and y-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 (x- and y-directions) and vertical (z-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.Figure 3
Schematic view of the stress-dependent fluid flow test apparatus.Figure 4
(a) A top view of hydraulic setup of the test apparatus. (b) Loading diagram of the rock fracture specimen.
(a) (b)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 thex-direction (Fx) and y-direction (Fy) were applied and increased from 7 to 35 kN in a 7 kN interval, with a fixed lateral pressure coefficient λ of Fy to Fx of 1.0. Applications of the boundary load conditions are displayed in Figure 4(b). In this study, we just discuss the x-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 (x-direction) were opened while the top and bottom valves (y-direction) were closed, which forced water to flow horizontally from the right to left specimen boundary (x-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−3Pa·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 variousN subjected to different F (here, F refers to the boundary load due to Fx = Fy), 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 P of 0–0.6 MPa, J 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 (J) versus discharge (Q) 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 R2 values all larger than 0.99 (Table 1). From Figure 5, with an increase in J, Q shows an increasing trend. For a certain N, as F increases, the slope of J-Q 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 F, as N increases, Q flowing out of the fracture networks produced by an identical hydraulic head shows an increasing trend (Figure 6). Taking F = 14 kN as an example, Q produced by J = 123.69 is increased from 5.04 × 10−6m3/s (N=1) to 10.86 × 10−6m3/s (N=12), increasing by approximately a factor of 1.15. With an increase in N, the overall discharge capacity of the fracture networks enhances.Table 1
Values ofa, b, Jc, and R2 of fracture networks with different numbers of intersections N subjected to various boundary loads F.
N F (kN) a(kg⋅Pa−1⋅s−1⋅m−4) b(kg⋅Pa−1⋅m−7) R 2 J c 1 7 8.14 E + 06 1.36 E + 12 0.9993 6.01 14 1.37 E + 07 2.15 E + 12 0.9994 10.78 21 2.84 E + 07 3.93 E + 12 0.9993 25.34 28 5.74 E + 07 7.13 E + 12 0.9994 57.05 35 9.80 E + 07 1.00 E + 13 0.9996 118.57 4 7 4.13 E + 06 9.17 E + 11 0.9993 2.29 14 7.38 E + 06 1.40 E + 12 0.9990 4.80 21 1.90 E + 07 2.94 E + 12 0.9995 15.16 28 3.25 E + 07 4.80 E + 12 0.9995 27.17 35 6.19 E + 07 8.01 E + 12 0.9994 59.06 7 7 3.16 E + 06 7.57 E + 11 0.9991 1.63 14 3.57 E + 06 8.03 E + 11 0.9995 1.96 21 1.25 E + 07 2.01 E + 12 0.9994 9.58 28 2.26 E + 07 3.42 E + 12 0.9996 18.44 35 2.56 E + 07 3.76 E + 12 0.9975 21.50 12 7 1.39 E + 06 3.81 E + 11 0.9994 0.62 14 3.03 E + 06 7.80 E + 11 0.9989 1.45 21 9.70 E + 06 1.67 E + 12 0.9995 6.96 28 1.11 E + 07 1.82 E + 12 0.9994 8.36 35 1.69 E + 07 2.71 E + 12 0.9996 13.01Figure 5
Regression analysis of hydraulic gradient as a function of measured flow rate using Forchheimer’s law of fracture networks subjected to increased boundary loads.
(a)
N = 1 (b)
N = 4 (c)
N = 7 (d)
N = 12Figure 6
Q as a function of N for different F at a fixed J of 123.69.The regression linear and nonlinear coefficientsa and b 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 a and b in terms of N and F. For a certain N, as F increases, both a and b show an increasing trend. In the range of F from 7 to 35 kN, a increases by 7.10–13.99 times and b increases by 3.97–7.73 times. The increase of coefficient a value is due to fracture closure induced by the boundary loads. However, for a certain F, as N increases, both a and b decrease. The variation process can be divided into two stages. For N = 1–7, the coefficients decrease significantly as N varies. However, in the range of 7–12, the coefficients vary gradually and approach essentially constant values. Using coefficient b as an example, as N increases from 1 to 12, b decreases by 71.91% (F = 7 kN), 63.67% (F = 14 kN), 57.51% (F = 21 kN), 74.47% (F = 28 kN), and 72.90% (F = 35 kN), respectively.Figure 7
Variations in linear and nonlinear coefficientsa and b, and the nonlinear relation between a and b.
(a)
N-a (b)
F-a (c)
N-b (d)
F-b (e)
a-bIn view of the fact that the variations of linear and nonlinear coefficientsa and b against N and F are very similar, b as a function of a was analyzed [29], as plotted in Figure 7(e). The experimental data can be fitted very well using the following empirical function:(14)b=7.899×106a0.766.Although (14) fits very well the relationship between the coefficients a and b, with the generating R2 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 ofE for all test cases were calculated and plotted as a function of J based on (8), as shown in Figure 8. As J increases, E displays an increasing trend, while the rate of its increase steadily diminishes. The variation process can be well described using a power function(15)E=cJd,where the fitting coefficients (c and d) are related to both the boundary load conditions and number of intersections. For a certain N, as F increases, the coefficient c shows a decreasing trend, but the coefficient d increases. However, at a given F, with an increase in N, c plots an increasing trend while d decreases.Figure 8
Evolution of nonlinear effect factorE with hydraulic gradient J.
(a)
N = 1 (b)
N = 4 (c)
N = 7 (d)
N = 12By defining a critical value ofE = 0.1, Jc of the fracture networks with various numbers of intersections was calculated, as listed in Table 1. Figure 9 displays the changes in Jc due to the changes in N and F. As F increases, Jc shows an increasing trend. In the range of F from 7 to 35 kN, Jc increases by 6.35 times (N=1), 7.73 times (N=4), 3.97 times (N=7), and 6.09 times (N=12), respectively. The main reason for these variations may be as follows. For Fx = Fy, 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 Q and bh3 in (4), which would then impact Jc for the onset of flow transition. However, for a certain F, as N increases, Jc shows a decreasing trend. Taking F = 35 kN as an example, in the range of N from 1 to 12, Jc decreases from 118.57 to 13.01, or by 89.03%.Figure 9
Relationships betweenJc and N.For fluid flow through fractured and porous media, the transmissivity has also been applied to estimate the nonlinear flow regime. Using (11), T of the fracture networks was calculated. Figure 10 shows the relationships between T and J. Notably, T is not a constant value but exhibits a decreasing trend with J, which further validates the existence of nonlinear flow behaviors in the fracture networks. With an increase in F, T decreases due to fracture closure. However, for a certain F, T shows an increasing trend with N. Using the least square method, J as a function of T can be well fitted with the following exponential function:(16)T=A+BCJ,where the regression coefficients A, B, and C were all presented in Figure 10. The units of coefficients A and B are m4 and C is a dimensionless constant.Figure 10
Regression analysis of transmissivity (T) as a function of hydraulic gradient J.
(a)
N = 1 (b)
N = 4 (c)
N = 7 (d)
N = 12When 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 (T/T0) versus Re can be determined by extending the linear Darcy law with T=μQ/(-∇P) [23, 25, 34, 43],(17)TT0=11+βRe,where β is a dimensionless coefficient determined by the coefficients a and b 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,J in a model typically has known values and can be easily obtained. Therefore, in this study, the relationships between J and T/T0 were used to evaluate the flow regimes in the fracture networks.By fitting the experimental data sets, the relationships betweenJ and T/T0 for rock fracture networks in this study were analyzed, as shown in Figure 11. Here, T0 denotes the transmissivity corresponding to J = 0 in Figure 10, where the hydraulic head difference is sufficiently low and the inertial forces are negligible. The results show that variations in T/T0 against J can be expressed as follows [14]:(18)TT0=1-exp⁡-αJ-0.45,where α is a dimensionless coefficient. Note that the exponent of J is −0.45, which does not change with the number of intersections or boundary load conditions.Figure 11
Relationships between hydraulic gradientJ and the normalized transmissivity (T/T0).
(a)
N = 1 (b)
N = 4 (c)
N = 7 (d)
N = 12From Figure11, as J increases, T/T0 shows a decreasing trend. For a certain N, as F increases, the transmissivity relationships generally shift upward. However, for a certain F, the J – T/T0 curves shift downward as N increases. In addition, the variations in T/T0 with J can be divided into three stages. When J is small (i.e., less than 0.2), with the increase of J, T/T0 approximately holds a constant value of 1.0; thus, the fluid flow is linear. Then, with continuously increasing J, T/T0 decreases and the reduction rate of T/T0 first increases and then decreases. Based on (12) and (18), when T/T0 = 0.9, Jc is calculated, which is in the ranges 8.98–103.16 (N=1), 3.46–61.40 (N=4), 2.69–28.37 (N=7), and 1.10–15.50 (N=12), respectively. Generally, the range of Jc is larger than that calculated based on (6) and (8) shown in Figure 9.The variations of the coefficientα with varying N and F are shown in Figure 12. For a certain N, as F increases, α increases. In the range of 7–35 kN, α generally increases by a factor of 2. However, at a given F, the coefficient α decreases with N. As N increases from 1 to 12, α shows a decrease of 61.14% (F = 7 kN), 66.37% (F = 14 kN), 42.98% (F = 21 kN), 55.29% (F = 28 kN), and 57.38% (F = 35 kN), respectively.Figure 12
Variations in the coefficientα with varying N and F.
(a)
F-α (b)
N-αFrom the quadratic relationships betweenQ and J of the fracture networks displayed in Figure 5, at a given J, the increase in Q varies with both F and N, which leads to the changes in the equivalent permeability of the fracture networks. To assess the discharge capacity of the fracture networks, the x-directional equivalent permeability k of the fracture networks was calculated at a specified inlet hydraulic pressure of 0.05 MPa, with the corresponding J 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:(19)Q=Akμ∂P∂x,where A is the cross-sectional area.The changes in the equivalent permeability of the fracture networks in terms ofN and F are displayed in Figure 13. For a certain N, as F increases, k of the model region decreases. In the range of F from 7 to 35 kN, k decreases by 90.31% (N=1), 90.87% (N=4), 82.28% (N=7), and 84.82% (N=12), respectively. In addition, the F-k variation curves can be divided into two stages. As F = 7–21 kN, the permeability changes are sensitive. However, when F > 21 kN, k changes gradually and approaches constant values. The main reason for these variations may be as follows. As F increases, the effective stress σfne 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 givenF, the equivalent permeability k increases with N mainly due to richer connectivity in the fracture networks.Figure 13
Equivalent permeabilityk changes due to the changes in F.


## 5. Conclusions

This paper experimentally examines the impacts of number of intersections (N = 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 coefficienta and nonlinear coefficient b in Forchheimer’s law decrease with the number of intersections but increase with the boundary load. The changes of nonlinear coefficient b as a function of linear coefficient a can be well fitted using a power function b=7.899×106a0.766 based on the experimental data.(2)
The critical hydraulic gradient of the fracture networks is calculated by taking a critical nonlinear effect factorE 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 expressionT/T0=1-exp⁡(-αJ-0.45). 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.



---

*Source: 8217921-2018-04-24.xml*