# A Review of Critical Conditions for the Onset of Nonlinear Fluid Flow in Rock Fractures

Authors: Liyuan Yu; Richeng Liu; Yujing Jiang
Journal: Geofluids
 (2017)

---

## Abstract

Selecting appropriate governing equations for fluid flow in fractured rock masses is of special importance for estimating the permeability of rock fracture networks. When the flow velocity is small, the flow is in the linear regime and obeys the cubic law, whereas when the flow velocity is large, the flow is in the nonlinear regime and should be simulated by solving the complex Navier-Stokes equations. The critical conditions such as critical Reynolds number and critical hydraulic gradient are commonly defined in the previous works to quantify the onset of nonlinear fluid flow. This study reviews the simplifications of governing equations from the Navier-Stokes equations, Stokes equation, and Reynold equation to the cubic law and reviews the evolutions of critical Reynolds number and critical hydraulic gradient for fluid flow in rock fractures and fracture networks, considering the influences of shear displacement, normal stress and/or confining pressure, fracture surface roughness, aperture, and number of intersections. This review provides a reference for the engineers and hydrogeologists especially the beginners to thoroughly understand the nonlinear flow regimes/mechanisms within complex fractured rock masses.


---

## Body

## 1. Introduction

Rock fracture network controls the main paths of fluid flow and contaminant migration in deep underground, and the estimation of permeability of fractured rock masses has been extensively studied during the past several decades in many geoengineering and geosciences such as CO2 sequestration, enhanced oil recovery, and geothermal energy development [1–8]. The fluid flow in rock fractures and/or fracture networks is commonly assumed to obey the cubic law, in which the flow rate is linearly proportional to the pressure drop [9–13]. However, in the karst systems and/or in the vicinity of wells during pump tests, when the flow rate/velocity is large, fluid flow enters the nonlinear flow regime and flow rate is nonlinearly correlated with pressure drop [11, 14–16]. In such case, using the cubic law to calculate fluid flow will overestimate the conductivity of rock fractures and/or fracture networks [17–20]. Therefore, a thorough understanding of the nonlinear flow properties of fluid within fractures contributes to accurately predicting permeability of fractured rock masses [21].Previous studies have reported that there are three representative types of nonlinear flow behaviors in rock fractures induced by inertial effect, fracture dilation, and solid-water interaction [16, 22, 23], and the present study focuses on the nonlinear flow behaviors induced by inertial effect. Many factors can affect the magnitude of permeability of fractured rock masses, including fracture length [24–27], aperture [28–30], surface roughness [31, 32], dead-end [33], number of intersections [34, 35], hydraulic gradient [36], boundary stress [37, 38], anisotropy [39–42], scale [43–46], stiffness [47], coupled thermo-hydro-mechanical-chemical (HTMC) processes [48–51], and precipitation-dissolution and biogeochemistry [52]. The discrete fracture network (DFN) model, which can consider most of the above parameters, has been increasingly utilized to simulate fluid flow in the complex fractured rock masses [53–56], although it cannot model the aperture heterogeneity of each fracture [57–59]. In the numerical simulations and/or analytical analysis, the linear governing equation such as the cubic law is solved to simulate fluid flow in fractures by applying constant hydraulic gradients (J) on the two opposing boundaries, such as J=1 [53, 60–64], J=0.1 [37], J=0.001 [65, 66], and J = unknown constants [10, 30, 42, 67–69]. This assumption that fluid flow obeys the cubic law is suitable for characterizing hydraulic behaviors of deep underground engineering, in which the flow rate is sufficiently small. For the fractured rock masses such as the karst system and/or the in situ hydraulic tests, the flow rate is relatively high and nonlinearly correlated with the pressure drop [14, 70]. With increasing pressure drop, the fracture network permeability decreases. To accurately simulate the nonlinear flow in fractures, the Navier-Stokes (NS) equations should be solved, which are composed of a set of coupled nonlinear partial derivatives of varying orders [17, 71, 72] and are more complex than solving the cubic law. When fluid flow is in the linear regime, each fracture in the two-dimensional DFNs is represented using a line segment [73], whereas when fluid flow enters the nonlinear regime and the NS equations are solved, the real geometry (void space) of each fracture that is formed with two walls should be incorporated, which to some extent increases the difficulty of establishing the models [36, 74]. As a result, both the yearly published papers and yearly cited times with the keywords “nonlinear flow” and “rock mass” are much smaller than those with the keywords “linear flow” and “rock mass,” as shown in Figure 1. With the development of computing power, more researches are contributing to the nonlinear flow characteristics of fractures, which needs solving the NS equations and is often time-consuming. Therefore, it is a vital issue about how to determine the critical condition (i.e., critical Reynolds number or critical hydraulic gradient) for the onset of nonlinear flow.Figure 1
Yearly published papers and cited times with the keywords “linear flow and rock mass” and “nonlinear flow and rock mass.” These data are collected from the Web of Science™ Core Collection (SCI-Expanded) over the past twenty years from 1997 to 2016.This study firstly reviews the governing equations of fluid flow in fractures and then reviews the effects of shear displacement, normal stress and/or confining pressure, fracture surface roughness, fracture aperture, and fracture intersection on nonlinear flow characteristics of fractured rock masses. This work aims at providing a reference for engineers and researchers to quickly assess the magnitudes of critical Reynolds number or critical hydraulic gradient and to clearly understand the nonlinear flow mechanisms within complex fractured rock masses.


## 2. Governing Equations of Fluid Flow in Fractures

### 2.1. Navier-Stokes Equations

The flow of incompressible Newtonian fluid is governed by the NS equations, written in a tensor form as [75–77](1)ρ∂ui∂t+uiui,j=ρfi-P,i+μui,jj,where ui is the velocity of fluid in the i-direction in the Cartesian coordinate system with i=x,y,z and j=x,y,z, respectively, the body force is fi=gi=0,0,-g, ρ is the fluid density, t is the time, P is the pressure, and μ is the dynamic viscosity.For steady-state flow, the parameters related to timet can be ignored and (1) gives(2)ρuiui,j=ρfi-P,i+μui,jj.Hydraulic head(h) is defined as the summation of elevation height (ze) and hydraulic pressure head (P/(ρg)), as follows:(3)h=ze+Pρg.Hydraulic gradient (J) is defined as the ratio of hydraulic head difference to flow length, written as(4)Ji=h,i=1ρg-ρgi+P,i.Substitution of (4) into (2) leads to a more abbreviated form of NS equations, as below [78]:(5)ρuiui,j=-ρgh,i+μui,jj.


### 2.2. Stokes Equations

When the flow rate/velocity is very small, which corresponds to a small Reynolds number (Re), the inertial term (the left term of (5)) can be negligible, resulting in the expression of the Stokes equations [79]:(6)h,i=μρgui,jj.Equation (6) can be unfolded to the following forms:(7)h,x=μρgux,xx+ux,yy+ux,zz,h,y=μρguy,xx+uy,yy+uy,zz,h,z=μρguz,xx+uz,yy+uz,zz.


### 2.3. Reynolds Equation

Dimensional analysis is performed and some definitions are made:Λ is the dimension of fracture in xy-plane, U is the dimension of flow velocity in both x- and y-directions, V is the dimension of flow velocity in z-direction, and a is the dimension of fracture aperture.(8)mag∂2ux∂x2≈mag∂2ux∂y2≈UΛ2,mag∂2ux∂z2≈Ua2,mag∂2uy∂x2≈mag∂2uy∂y2≈UΛ2,mag∂2uy∂z2≈Ua2,mag∂2uz∂x2≈mag∂2uz∂y2≈UΛ2,mag∂2uz∂z2≈Ua2.When the dimension of fracture aperture is much less than the dimension of fracture in thexy-plane (a≪Λ), in which case the tortuosity of fracture void space is very small, the first two terms of (8) is much less than the last term and can be negligible. As a result, (7) can be simplified to the following equations:(9)h,x=μρgux,zz,(10)h,y=μρguy,zz,(11)h,z=μρguz,zz.When the fracture aperture and the change in hydraulic head in thez-direction are small, (11) yields(12)h,z≈0.Substitution of (12) into (9)~(11) gives(13)uxx,y,z=ρg2μh,xz2+b1z+c1,uyx,y,z=ρg2μh,yz2+b2z+c2,h=hx,y,where b1,b2,c1, and c2 are the integration constants.Importing non-slip boundary conditions results in(14)ux=0,uy=0,whenz=±E2,where E is the fracture aperture, and z=±E/2 are the lower and upper bounds of fracture void space in the z-direction.Substitution of (14) into (13) leads to the values of integration constants (b1,b2,c1, and c2), and (13) can be rewritten as(15)uxx,y,z=ρg2μh,xz+E2z-E2,(16)uyx,y,z=ρg2μh,yz+E2z-E2,(17)h=hx,y.Integrations of (15) and (16) along the z-direction leads to(18)qx=∫-E/2E/2uxdz=-ρgE312μ∂h∂x,qy=∫-E/2E/2uydz=-ρgE312μ∂h∂y,where qx is the flow rate in the x-direction, and qy is the flow rate in the y-direction.The law of conservation of mass of flow rate is expressed as [80](19)qx,x+qy,y=0.Thus, substituting (18) into (19) results in the Reynolds equation, given as [17](20)∂∂xρgE312μ∂h∂x+∂∂yρgE312μ∂h∂y=0.


### 2.4. Cubic Law

An assumption is made that fluid flows through two infinite smooth parallel plates with a constant apertureE, as shown in Figure 2. Because fluid flows along the x-coordinate, the flow velocities in the y- and z-directions equal zero.(21)uy=uz=0.Figure 2
Flow velocity distribution in a parallel plate model.The continuity equation of fluid is(22)∂ρ∂t+ρui,i=0.For incompressible fluid, (22) can be simplified to(23)ui,i=0.Substitution of (21) into (5) gives the governing equation of fluid flow in thex-direction, written as(24)ρuxux,x=-ρgh,x+μux,xx.Substituting (23) into (24) gives rise to the fact that the left term equals zero.(25)0=-ρgh,x+μux,xx.Importing the non-slip boundary conditions (see (14)) and integrating (25) lead to(26)Q=∫-E/2E/2uxdz=-E312μ∂P∂x,where Q is the flow rate through a fracture.Equation (26) is the expression for the cubic law, in which fracture width (w) is assigned a constant value of 1. When w≠1,w should be introduced in (26), as follows:(27)Q=-wE312μ∂P∂x.


### 2.5. Forchheimer Equation

With increasing flow velocity orRe, the flow rate is nonlinearly correlated with pressure drop and (27) is no longer applicable. Some empirical expressions were proposed to describe this nonlinear relationship, among which the Forchheimer equation is widely used where the pressure drop is a quadratic function of flow rate, expressed as [82–87](28)-∇P=AQ+BQ2,where A and B are two coefficients that represent the linear and nonlinear terms of fluid flow, respectively. When Q is very small, the nonlinear term (BQ2) can be negligible and drops out. In this case, (28) reduces to(29)-∇P=AQ,where A is a constant that equals 12μ/wE3 if the fracture is not deformable. Equation (28) can probably be applied over the entire range of flow rate/velocity, including that it reduces to linear Darcy’s law (29) at a sufficiently low flow rate/velocity. Another macroscopic nonlinear flow equation to characterize the nonlinear flow in rock fractures is Izbash’s law (Izbash 1931), but it cannot well quantify the linear flow properties at a sufficiently low Re. Therefore, the two equations can both be utilized to describe the nonlinear relationships between flow rate and pressure drop, especially for strong inertial regime; however, only the Forchheimer equation is used to depict the critical Reynolds number by directly plotting the normalized transmissivity–Reynolds number curves [88–90].To describe the nonlinearity of fluid flow, a nonlinear factor (Nf) is defined as(30)Nf=BQ2AQ+BQ2.N f denotes the ratio of pressure drop caused by the nonlinear term (BQ2). In the previous works, it is assumed that the critical Reynolds number (Rec) is the Re that corresponds to Nf = 0.1 [88, 91–93].Besides the dimensionlessRe that is used for characterizing the transition from linear to nonlinear flow regimes, the dimensionless Forchheimer number (F0) is another widely accepted parameter, which is defined as the ratio of nonlinear to linear pressure losses, written as [16, 90, 94](31)F0=BQ2AQ=BQA.Thus, (30) can be rewritten as(32)Nf=F01+F0.The transmissivity (T), which is a commonly used parameter to describe the conductivity of a fracture, is defined as(33)T=E312=-μQw∇P.When fluid flow is in the linear regime,T is a constant and is represented using T0. Thus the ratio of T to T0 can be calculated according to(34)TT0=-μQ/wAQ+BQ2-μQ/wAQ=AQAQ+BQ2=1-Nf.Therefore,T/T0 = 0.9 and Nf = 0.1 have the same physical meaning that the nonlinear term (BQ2) contributes to 10% of the pressure drop, in which the current Re is Rec.


## 2.1. Navier-Stokes Equations

The flow of incompressible Newtonian fluid is governed by the NS equations, written in a tensor form as [75–77](1)ρ∂ui∂t+uiui,j=ρfi-P,i+μui,jj,where ui is the velocity of fluid in the i-direction in the Cartesian coordinate system with i=x,y,z and j=x,y,z, respectively, the body force is fi=gi=0,0,-g, ρ is the fluid density, t is the time, P is the pressure, and μ is the dynamic viscosity.For steady-state flow, the parameters related to timet can be ignored and (1) gives(2)ρuiui,j=ρfi-P,i+μui,jj.Hydraulic head(h) is defined as the summation of elevation height (ze) and hydraulic pressure head (P/(ρg)), as follows:(3)h=ze+Pρg.Hydraulic gradient (J) is defined as the ratio of hydraulic head difference to flow length, written as(4)Ji=h,i=1ρg-ρgi+P,i.Substitution of (4) into (2) leads to a more abbreviated form of NS equations, as below [78]:(5)ρuiui,j=-ρgh,i+μui,jj.


## 2.2. Stokes Equations

When the flow rate/velocity is very small, which corresponds to a small Reynolds number (Re), the inertial term (the left term of (5)) can be negligible, resulting in the expression of the Stokes equations [79]:(6)h,i=μρgui,jj.Equation (6) can be unfolded to the following forms:(7)h,x=μρgux,xx+ux,yy+ux,zz,h,y=μρguy,xx+uy,yy+uy,zz,h,z=μρguz,xx+uz,yy+uz,zz.


## 2.3. Reynolds Equation

Dimensional analysis is performed and some definitions are made:Λ is the dimension of fracture in xy-plane, U is the dimension of flow velocity in both x- and y-directions, V is the dimension of flow velocity in z-direction, and a is the dimension of fracture aperture.(8)mag∂2ux∂x2≈mag∂2ux∂y2≈UΛ2,mag∂2ux∂z2≈Ua2,mag∂2uy∂x2≈mag∂2uy∂y2≈UΛ2,mag∂2uy∂z2≈Ua2,mag∂2uz∂x2≈mag∂2uz∂y2≈UΛ2,mag∂2uz∂z2≈Ua2.When the dimension of fracture aperture is much less than the dimension of fracture in thexy-plane (a≪Λ), in which case the tortuosity of fracture void space is very small, the first two terms of (8) is much less than the last term and can be negligible. As a result, (7) can be simplified to the following equations:(9)h,x=μρgux,zz,(10)h,y=μρguy,zz,(11)h,z=μρguz,zz.When the fracture aperture and the change in hydraulic head in thez-direction are small, (11) yields(12)h,z≈0.Substitution of (12) into (9)~(11) gives(13)uxx,y,z=ρg2μh,xz2+b1z+c1,uyx,y,z=ρg2μh,yz2+b2z+c2,h=hx,y,where b1,b2,c1, and c2 are the integration constants.Importing non-slip boundary conditions results in(14)ux=0,uy=0,whenz=±E2,where E is the fracture aperture, and z=±E/2 are the lower and upper bounds of fracture void space in the z-direction.Substitution of (14) into (13) leads to the values of integration constants (b1,b2,c1, and c2), and (13) can be rewritten as(15)uxx,y,z=ρg2μh,xz+E2z-E2,(16)uyx,y,z=ρg2μh,yz+E2z-E2,(17)h=hx,y.Integrations of (15) and (16) along the z-direction leads to(18)qx=∫-E/2E/2uxdz=-ρgE312μ∂h∂x,qy=∫-E/2E/2uydz=-ρgE312μ∂h∂y,where qx is the flow rate in the x-direction, and qy is the flow rate in the y-direction.The law of conservation of mass of flow rate is expressed as [80](19)qx,x+qy,y=0.Thus, substituting (18) into (19) results in the Reynolds equation, given as [17](20)∂∂xρgE312μ∂h∂x+∂∂yρgE312μ∂h∂y=0.


## 2.4. Cubic Law

An assumption is made that fluid flows through two infinite smooth parallel plates with a constant apertureE, as shown in Figure 2. Because fluid flows along the x-coordinate, the flow velocities in the y- and z-directions equal zero.(21)uy=uz=0.Figure 2
Flow velocity distribution in a parallel plate model.The continuity equation of fluid is(22)∂ρ∂t+ρui,i=0.For incompressible fluid, (22) can be simplified to(23)ui,i=0.Substitution of (21) into (5) gives the governing equation of fluid flow in thex-direction, written as(24)ρuxux,x=-ρgh,x+μux,xx.Substituting (23) into (24) gives rise to the fact that the left term equals zero.(25)0=-ρgh,x+μux,xx.Importing the non-slip boundary conditions (see (14)) and integrating (25) lead to(26)Q=∫-E/2E/2uxdz=-E312μ∂P∂x,where Q is the flow rate through a fracture.Equation (26) is the expression for the cubic law, in which fracture width (w) is assigned a constant value of 1. When w≠1,w should be introduced in (26), as follows:(27)Q=-wE312μ∂P∂x.


## 2.5. Forchheimer Equation

With increasing flow velocity orRe, the flow rate is nonlinearly correlated with pressure drop and (27) is no longer applicable. Some empirical expressions were proposed to describe this nonlinear relationship, among which the Forchheimer equation is widely used where the pressure drop is a quadratic function of flow rate, expressed as [82–87](28)-∇P=AQ+BQ2,where A and B are two coefficients that represent the linear and nonlinear terms of fluid flow, respectively. When Q is very small, the nonlinear term (BQ2) can be negligible and drops out. In this case, (28) reduces to(29)-∇P=AQ,where A is a constant that equals 12μ/wE3 if the fracture is not deformable. Equation (28) can probably be applied over the entire range of flow rate/velocity, including that it reduces to linear Darcy’s law (29) at a sufficiently low flow rate/velocity. Another macroscopic nonlinear flow equation to characterize the nonlinear flow in rock fractures is Izbash’s law (Izbash 1931), but it cannot well quantify the linear flow properties at a sufficiently low Re. Therefore, the two equations can both be utilized to describe the nonlinear relationships between flow rate and pressure drop, especially for strong inertial regime; however, only the Forchheimer equation is used to depict the critical Reynolds number by directly plotting the normalized transmissivity–Reynolds number curves [88–90].To describe the nonlinearity of fluid flow, a nonlinear factor (Nf) is defined as(30)Nf=BQ2AQ+BQ2.N f denotes the ratio of pressure drop caused by the nonlinear term (BQ2). In the previous works, it is assumed that the critical Reynolds number (Rec) is the Re that corresponds to Nf = 0.1 [88, 91–93].Besides the dimensionlessRe that is used for characterizing the transition from linear to nonlinear flow regimes, the dimensionless Forchheimer number (F0) is another widely accepted parameter, which is defined as the ratio of nonlinear to linear pressure losses, written as [16, 90, 94](31)F0=BQ2AQ=BQA.Thus, (30) can be rewritten as(32)Nf=F01+F0.The transmissivity (T), which is a commonly used parameter to describe the conductivity of a fracture, is defined as(33)T=E312=-μQw∇P.When fluid flow is in the linear regime,T is a constant and is represented using T0. Thus the ratio of T to T0 can be calculated according to(34)TT0=-μQ/wAQ+BQ2-μQ/wAQ=AQAQ+BQ2=1-Nf.Therefore,T/T0 = 0.9 and Nf = 0.1 have the same physical meaning that the nonlinear term (BQ2) contributes to 10% of the pressure drop, in which the current Re is Rec.


## 3. Nonlinear Flow Characteristics of Rock Fractures

### 3.1. Effect of Shear Displacement

Under a constant normal stress, increasing shear displacement (us) can change the void space of a fracture. As a result, both the average aperture and permeability/transmissivity increase by up to 3 orders of magnitude with a maximum of us of 16 mm [95–99]. This variation of fracture void space during shear can also change the nonlinear flow properties of fluid such as the critical Reynolds number Rec, which is used for characterizing the onset of nonlinear flow. Xia et al. [81] conducted laboratory experiments on three artificial rock joints having natural joint surface characteristics to investigate the nonlinear flow behaviors under different shear displacements. The cylindrical specimens that were cored from marble blocks taken from the construction site of Jin-Ping-II Hydropower Station in China have a diameter of 50 mm and a length of 100 mm. The average values of joint roughness coefficient (JRC) of the three samples are 13.28, 15.09, and 11.23, respectively. By fitting their results, it is found that when Re is small (i.e., less than 10), with the increment of Re, T/T0 holds approximately constant value of 1 (see Figure 3). With continuously increasing Re, T/T0 decreases. When T/T0=0.9, Rec is calculated, which is in the ranges 71.33~340.01 (us = 1~5 mm), 75.93~242.67 (us = 1~6 mm), and 97.80~326.07 (us = 1~5 mm), respectively. The range of Rec (from 71.33 to 340.01) is much less than that (from 1408.2 to 5674.4) calculated in the original works, because they used T/T0 = 0.1 to quantify the onset of nonlinear flow, in which the nonlinear term contributes to 90% of the pressure drop. They also reported that, with increasing us, both the linear and nonlinear coefficients (A and B) in Forchheimer equation decrease. Xiong et al. [95] simulated the coupled shear-flow behaviors by directly solving the NS equations using the geometries of two natural rock fractures labelled with J3 and J7, whose JRC equals 17~18 and 3~4, respectively. The relationships betweenT/T0 and Re as shown in Figure 4 exhibit the same variation trend with those in Figure 3. With increasing us from 4 mm to 10 mm, Rec increases from 9.55 to 18.07, because the larger us gives rise to the larger fracture aperture and a smaller number of contacts, which results in the increment of Rec. Javadi et al. [90] experimentally investigated the role of shear processes on Rec for nonlinear flow through rough-walled fractures. The normal stress (σn) ranges from 1.0 to 5.0 MPa and us = 0~20 mm. The results show that when us = 0~1 mm, Rec has very small values approaching 0 for all cases (see Figure 5). When us = 1~5 mm, Rec dramatically increases for σn = 1.0 MPa and 5.0 MPa but gently increases for σn = 3.0 MPa. For us = 5~20 mm, the Rec holds approximately constant values for σn = 5.0 MPa and slightly decreases forσn = 1.0 MPa; however the Rec continuously increases for σn = 3.0 MPa. The different variations of Rec~us curves may be induced by the statistical distributions of void space of fractures during shear. The range of Rec is from 0.001 to 25. Rong et al. [100] prepared 6 granite fractures with JRC = 6.67~8.18 and conducted coupled shear-flow tests with σn = 0.5~3.0 MPa and us = 0.0~10.9 mm. The results shown in Figure 6 depict that when us = 0~0.5 mm, with the increment of us, Rec decreases, which is very different from those presented in Figure 5. With increasing us from 0.5 to 3.0 mm, Rec increases dramatically and then varies negligibly small when us exceeds 3.0 mm. This variation trend is similar to that (with σn = 1.0 or 5.0 MPa) reported in Javadi et al. [90]. As an example, when σn = 1.0 MPa, Rec that corresponds to a small us in the work of Javadi et al. [90] is much smaller than that in the work of Rong et al. [100]; however when us is relatively large, Rec in the work of Javadi et al. [90] is approximately 2.0~2.5 times that in the work of Rong et al. [100]. The differences of Rec evolution are dependent on fracture’s morphology such as mean JRC. Zimmerman et al. [88] presented a high-resolution NS simulation and laboratory measurements of fluid flow in a natural sandstone fracture. The two opposing fracture surfaces were made using epoxy casts [101, 102], which has a size of 2 cm × 2 cm. By fitting the tested and computed results, Figure 7 depicts that Rec equals 12.56 for the experiment and 21.92 for the simulation, which is near the range of Rec reported in literature [90, 95].Figure 3
Relationships betweenT/T0 and Re based on the data replotted from Xia et al. ([81], Figure 10).
(a)
Specimen #1 (b)
Specimen #2 (c)
Specimen #3Figure 4
Relationships betweenT/T0 and Re based on the data replotted from Xiong et al. ([95], Figure 12(a)).Figure 5
Relationships betweenRec and us based on the data replotted from Javadi et al. ([90], Figure 8).Figure 6
Relationships betweenRec and us based on the data replotted from Rong et al. ([100], Figure 9).Figure 7
Relationships betweenT/T0 and Re based on the data replotted from Zimmerman et al. ([88], Figure 2).


### 3.2. Effects of Normal Stress and/or Confining Pressure

The normal stress applied on a cuboid sample and/or confining pressure applied on a cylindrical sample can decrease the permeability of rock mass, following exponential functions [103–108]. With increasing the normal stress σn, both the fracture closure un and the increasing rate increase as shown in Figure 8, where unmax is the maximum closure and kn0 is the initial stiffness [109, 110]. This aperture variation/closure caused by normal stress contributes to the magnitudes of Rec and/or Jc. Zhang and Nemcik [111] studied the fluid flow regimes and nonlinear flow characteristics in deformable rock fractures through flow tests on four cylindrical grain sandstone samples under confining stresses from 1.0 MPa to 3.5 MPa. They analyzed the variations of nonlinear factor B describing flow in deformable rock fractures and reported that the confining stress does not change both the linear and nonlinear flow patterns but has a significant influence on flow characteristics. With increasing JRC from 5.5 to 15.4, the range of Rec under different confining stresses is 21.96~40.34, 46.73~66.84, 40.35~51.70, and 26.13~32.38, respectively, as shown in Figure 9. The calculated ranges of Rec are several times larger than those calculated by Zhang and Nemcik [111], in which the corresponding Rec = 3.5~4.5, 13.1~17.6, 19.3~24.8, and 6.3~8.6, respectively. The potential reason may be that, in the study of Zhang and Nemcik [111], the flow has entered the strong nonlinear regimes, where the factor A calculated using the fitted Forchheimer equation cannot be accurately estimated. That means, in such a situation, changing the magnitude of A will have negligible effect on the correlation coefficient between the tested results and fitted curves of Forchheimer equation, but it can robustly change the variations of Nf, resulting in the smaller ranges of Rec. This is the reason why we would like to use the variations ofT/T0 rather than Nf in the present study. Ranjith and Darlington [89] conducted nonlinear single-phase flow tests on real cylindrical rock joints with 110 mm height and 55 mm diameter under different confining pressures. The results depicted that, with increasing σn from 1.0 MPa to 5.0 MPa, Rec decreases from 31.29 to 28.77 (see Figure 10), which is in a smaller range compared with those in Figure 9. Chen et al. [112] presented the high-pressure packer test (HPPT) observations performed in fractured rock masses, which is located at about 450 m in depth in Qiongzhong County, China. Four test boreholes numbered ZK124, ZK126, ZK129-1, and ZK129-2, respectively, were obtained from in situ cites in the depths ranging from 33.4 m to 142.1 m with numbers of test intervals from 6 to 16, and the flow rate-pressure drop relationships as well as the magnitude of Rec were calculated. By collecting the tested results and replotting them in Figure 11, the variations of Rec with depth show that, for all boreholes, Rec fluctuates significantly in the range 25~66 and does not exhibit strong relationships with depth. The reasons may be that the boreholes at different depths contain fractures having different geometric properties such as aperture, persistence, orientation, and surface roughness. To simplify the model and clearly investigate the influence of normal stress and/or confining pressure on Rec, Zhou et al. [94] prepared four granite rock samples (G1~G4) taken from the depth 450~550 m in the Beishan area in Gansu Province, China, and two sandstone rock samples (S1~S2) taken from the depth of about 20 m in the Three Gorges Reservoir area, Hubei Province, China. The confining pressure varies from 1.0 to 30.0 MPa and Rec under each confining pressure is calculated. The results depicted in Figure 12 show that, with increasing confining pressure, Rec firstly increases and then decreases for the samples G2~G4 and S1~S2, because when the confining pressure is small, the decreased aperture due to the elastic compression will increase Rec [36]. When continuously increasing confining pressure, the contacts between two fracture walls are destroyed and the microfractures may propagate within rocks and the split rock fragments will be clamped by the two walls, which gives rise to complex flow behaviors and decreases Rec. For the sample G1, with increasing confining pressure, Rec steadily decreases and shows a different variation trend with other samples. They reported that Rec for the samples G1~G4 and S1~S2 varies in the ranges 0.075~9.243, 0.120~4.467, 0.160~5.108, 0.039~4.509, 0.191~4.090, and 0.026~2.976, respectively. The results of Rong et al. [100] as reported in Section 3.1 show that when us = 0 mm, Rec steadily decreases, presenting the same trend with the sample G1 in Figure 12. However, for samples under shear with us > 0 mm, with the increment of confining pressure, Rec firstly increases and then decreases as shown in Figure 13, presenting the same trends with most of the samples in Figure 12. For σn = 1.0~3.0 MPa, the range of Rec (Rec = 1.50~13.03) in the work of Rong et al. [100] is a little larger than that (Rec = 0.16~9.24) of Zhou et al. [94].Figure 8
Relationships betweenσn and un with different kn0.Figure 9
Relationships betweenT/T0 and Re based on the data replotted from Zhang and Nemcik ([111], Figure 8).
(a)
Sample number 1, JRC = 5.5 (b)
Sample number 2, JRC = 7.0 (c)
Sample number 3, JRC = 9.2 (d)
Sample number 2, JRC = 15.4Figure 10
Relationships betweenT/T0 and Re based on the data replotted from Ranjith and Darlington ([89], Figure 8).Figure 11
Relationships between depth and Re based on the data replotted from Chen et al. ([112], Tables 2–5).Figure 12
Relationships betweenRec and σn based on the data replotted from Zhou et al. ([94], Figure 8).Figure 13
Relationships betweenRec and σn based on the data replotted from Rong et al. ([100], Figure 9).


### 3.3. Effect of Fracture Surface Roughness

Fracture surface roughness can alter the distributions of void spaces of fractures, especially during shear [97, 114–117]. Besides, the rougher fracture surface contributes to the longer flow paths that a particle moves within fractures, resulting in a weaker conductivity/permeability if the same pressure difference is applied on the opposing boundaries [32, 73, 118, 119]. The previous works have shown that the fracture surface roughness, especially the secondary roughness, plays a significant role on the nonlinear flow properties of rock fractures, because the eddy flow occurs due to the surface roughness [77]. Wang et al. [120] established a series of 3D self-affine rock fractures using the successive random additions (SPA) method, in which the geometry of fracture surface is fractal [121, 122]. The fractal dimension (D3) is incorporated to characterize fracture surface roughness, which is correlated with the Hurst exponent (H) asD3 = 3-H [123–127]. Then a 3D Lattice Boltzmann method (LBM) that has been proven to be capable of solving the NS equations [128–131] is adopted to characterize the nonlinear flow regimes in rock fractures. The results show that, with increasing fracture surface roughness represented by D3, Rec decreases from 47.29 to 3.78 following an exponential function (see Figure 14), because the fluid flow within the rougher fractures will enter the nonlinear flow regime from the linear regime at a lower Re. The range of Rec (3.78~47.29) corresponding to Nf = 0.1 orT/T0 = 0.9 is larger than that (1.8~22.0) shown in the work of Wang et al. [120], since they utilized Nf = 0.05 orT/T0 = 0.95 as the threshold. Liu et al. [36] studied the transition of fluid flow from the linear to the nonlinear regime and quantified the critical hydraulic gradient (Jc) for flow in simple crossed fractures and complex fracture networks. The results show that Jc decreases significantly when JRC < 5 and then decreases slightly when JRC ≥ 5, following power law functions as presented in Figure 15. Note that, in their study, the onset of nonlinear flow is categorized using Nf = 0.01 orT/T0 = 0.99 rather than the commonly used Nf = 0.10 orT/T0 = 0.90 due to their high-precision flow testing system that allows them to adopt a much lower critical condition. By fitting the tested and simulated results, a multivariable regression function was proposed to calculate Jc, written as(35)Jc=λE-0.5·exp⁡300.8λE-0.01-0.03JRC-0.3Ni0.5-303.3,where Ni is the number of intersections and λ is the coefficient that quantifies the reduction of E due to the variation of fracture apertures in a DFN.Figure 14
Relationships between (a)T/T0 and Re and (b) Rec andD3 based on the data replotted from Wang et al. ([47], Figure 16).
(a)
Relationship betweenT/T0 and Re (b)
Relationship betweenRec andD3Figure 15
Relationships betweenJ and Re and JRC based on the data replotted from Liu et al. ([36], Figure 6(c)).The predictedJc using (35) fits well with the simulation results of DFNs with well-known geometric characteristics of fractures, which verifies the validity of (35).


### 3.4. Effect of Fracture Aperture

In the cubic law as described in (27), flow rate is linearly proportional to the cubic of fracture aperture, indicating that fracture aperture plays a vital role in controlling fluid flow in fractures. The fracture aperture has been observed to follow lognormal distributions [132–134] and incorporated into the theoretical models [28, 135–138]. Some previous works show that fracture aperture exhibits a power law relationship with fracture length with the exponent varying from 0.5 to 2.0 [62, 139–143]. The nonlinear flow properties of fractures also depend on the magnitude of fracture aperture. Liu et al. [113] took four rock samples of granite (S1 and S2), marble (S3), and limestone (S4) from the underground caverns of the Huangdao State Oil Reserves, China, and scanned their surface geometries using a three-dimensional laser-scanning rock surface instrument [144]. A total of twenty 2D rock fractures were prepared and fracture aperture varied from 1.0 mm to 10.0 mm. Fluid flow was simulated by directly solving the NS equations, and the nonlinear flow regimes were investigated. The results show that, with the increasing J, the reduction rate ofT/T0 increases and then decreases, existing an inflection point that has been utilized to quantify the transition from linear to nonlinear flow regime (see Figure 16(a)). With the increment of fracture aperture (E) from 1.0 mm to 10.0 mm, Jc decreases approximately four orders of magnitude following a power law function as shown in  (36) and Figure 16(b).(36)Jc=0.0627×E-3.87.Figure 16
Relationships between (a)T/T0 and J and (b) Jc and E based on the data replotted from Liu et al. ([113], Figure 4).
(a)
Relationship betweenT/T0 and J (b)
Relationship betweenJc and EJ c calculated corresponding toT/T0 = 0.9 is smaller than that calculated corresponding to the inflection point. Based on the simple DFNs established by Liu et al. [36], the relationships between Jc and E were established as shown in Figure 17. With the increment of E from 0.5 mm to 10.0 mm, Jc decreases and spans about five orders of magnitude following power law relationships, despite the magnitude of Ni. The coefficient and exponent of the fitted equations are close to those in  (36). The results show that, for both single fractures and fracture networks, Jc and E have power law correlations with a negative exponent.Figure 17
Relationships betweenJc and E based on the data replotted from Liu et al. ([36], Figure 6(a)).


### 3.5. Effect of Fracture Intersection

Fluid is redistributed at the intersections depending on the number and geometry of fracture segments connected to the intersection [74]. Wilson and Witherspoon [145] experimentally investigated the magnitude of laminar flow interference effects at sing fracture intersections. They found that when flow is in the laminar flow regime, the inertial effects at intersections can be negligible and when Re = 100, the head loss equals about five times the pipe diameters. Kosakowski and Berkowitz [146] established numerical models with a wide variability of fracture intersection geometry and studied the streamline distributions as well as the validity of Darcy’s law, by solving the NS equations using the available FEATFLOW package [147]. Their results indicate that the nonlinear inertial effects are significant for Re = 1~100, which is the range that often exists in the karst systems [14] and/or in the vicinity of wells during pump tests [15]. As introduced in the previous Sections 3.3 and 3.4, Liu et al. [36] estimated the influence of number of fracture intersections on the basis of a series of simple and complex DFNs. The results show that, with increasing Ni from 1 to 12, Jc decreases significantly and then slowly showing power law functions as shown in Figure 18. The reduction rate of Jc for JRC = 20 is about five times faster than that for JRC = 0. Comparisons among Figures 15, 17, and 18 indicate that Jc is the most sensitive to fracture aperture, followed by the number of intersections and fracture surface roughness.Figure 18
Relationships betweenJc and Ni based on the data replotted from Liu et al. ([36], Figure 6(b)).


## 3.1. Effect of Shear Displacement

Under a constant normal stress, increasing shear displacement (us) can change the void space of a fracture. As a result, both the average aperture and permeability/transmissivity increase by up to 3 orders of magnitude with a maximum of us of 16 mm [95–99]. This variation of fracture void space during shear can also change the nonlinear flow properties of fluid such as the critical Reynolds number Rec, which is used for characterizing the onset of nonlinear flow. Xia et al. [81] conducted laboratory experiments on three artificial rock joints having natural joint surface characteristics to investigate the nonlinear flow behaviors under different shear displacements. The cylindrical specimens that were cored from marble blocks taken from the construction site of Jin-Ping-II Hydropower Station in China have a diameter of 50 mm and a length of 100 mm. The average values of joint roughness coefficient (JRC) of the three samples are 13.28, 15.09, and 11.23, respectively. By fitting their results, it is found that when Re is small (i.e., less than 10), with the increment of Re, T/T0 holds approximately constant value of 1 (see Figure 3). With continuously increasing Re, T/T0 decreases. When T/T0=0.9, Rec is calculated, which is in the ranges 71.33~340.01 (us = 1~5 mm), 75.93~242.67 (us = 1~6 mm), and 97.80~326.07 (us = 1~5 mm), respectively. The range of Rec (from 71.33 to 340.01) is much less than that (from 1408.2 to 5674.4) calculated in the original works, because they used T/T0 = 0.1 to quantify the onset of nonlinear flow, in which the nonlinear term contributes to 90% of the pressure drop. They also reported that, with increasing us, both the linear and nonlinear coefficients (A and B) in Forchheimer equation decrease. Xiong et al. [95] simulated the coupled shear-flow behaviors by directly solving the NS equations using the geometries of two natural rock fractures labelled with J3 and J7, whose JRC equals 17~18 and 3~4, respectively. The relationships betweenT/T0 and Re as shown in Figure 4 exhibit the same variation trend with those in Figure 3. With increasing us from 4 mm to 10 mm, Rec increases from 9.55 to 18.07, because the larger us gives rise to the larger fracture aperture and a smaller number of contacts, which results in the increment of Rec. Javadi et al. [90] experimentally investigated the role of shear processes on Rec for nonlinear flow through rough-walled fractures. The normal stress (σn) ranges from 1.0 to 5.0 MPa and us = 0~20 mm. The results show that when us = 0~1 mm, Rec has very small values approaching 0 for all cases (see Figure 5). When us = 1~5 mm, Rec dramatically increases for σn = 1.0 MPa and 5.0 MPa but gently increases for σn = 3.0 MPa. For us = 5~20 mm, the Rec holds approximately constant values for σn = 5.0 MPa and slightly decreases forσn = 1.0 MPa; however the Rec continuously increases for σn = 3.0 MPa. The different variations of Rec~us curves may be induced by the statistical distributions of void space of fractures during shear. The range of Rec is from 0.001 to 25. Rong et al. [100] prepared 6 granite fractures with JRC = 6.67~8.18 and conducted coupled shear-flow tests with σn = 0.5~3.0 MPa and us = 0.0~10.9 mm. The results shown in Figure 6 depict that when us = 0~0.5 mm, with the increment of us, Rec decreases, which is very different from those presented in Figure 5. With increasing us from 0.5 to 3.0 mm, Rec increases dramatically and then varies negligibly small when us exceeds 3.0 mm. This variation trend is similar to that (with σn = 1.0 or 5.0 MPa) reported in Javadi et al. [90]. As an example, when σn = 1.0 MPa, Rec that corresponds to a small us in the work of Javadi et al. [90] is much smaller than that in the work of Rong et al. [100]; however when us is relatively large, Rec in the work of Javadi et al. [90] is approximately 2.0~2.5 times that in the work of Rong et al. [100]. The differences of Rec evolution are dependent on fracture’s morphology such as mean JRC. Zimmerman et al. [88] presented a high-resolution NS simulation and laboratory measurements of fluid flow in a natural sandstone fracture. The two opposing fracture surfaces were made using epoxy casts [101, 102], which has a size of 2 cm × 2 cm. By fitting the tested and computed results, Figure 7 depicts that Rec equals 12.56 for the experiment and 21.92 for the simulation, which is near the range of Rec reported in literature [90, 95].Figure 3
Relationships betweenT/T0 and Re based on the data replotted from Xia et al. ([81], Figure 10).
(a)
Specimen #1 (b)
Specimen #2 (c)
Specimen #3Figure 4
Relationships betweenT/T0 and Re based on the data replotted from Xiong et al. ([95], Figure 12(a)).Figure 5
Relationships betweenRec and us based on the data replotted from Javadi et al. ([90], Figure 8).Figure 6
Relationships betweenRec and us based on the data replotted from Rong et al. ([100], Figure 9).Figure 7
Relationships betweenT/T0 and Re based on the data replotted from Zimmerman et al. ([88], Figure 2).


## 3.2. Effects of Normal Stress and/or Confining Pressure

The normal stress applied on a cuboid sample and/or confining pressure applied on a cylindrical sample can decrease the permeability of rock mass, following exponential functions [103–108]. With increasing the normal stress σn, both the fracture closure un and the increasing rate increase as shown in Figure 8, where unmax is the maximum closure and kn0 is the initial stiffness [109, 110]. This aperture variation/closure caused by normal stress contributes to the magnitudes of Rec and/or Jc. Zhang and Nemcik [111] studied the fluid flow regimes and nonlinear flow characteristics in deformable rock fractures through flow tests on four cylindrical grain sandstone samples under confining stresses from 1.0 MPa to 3.5 MPa. They analyzed the variations of nonlinear factor B describing flow in deformable rock fractures and reported that the confining stress does not change both the linear and nonlinear flow patterns but has a significant influence on flow characteristics. With increasing JRC from 5.5 to 15.4, the range of Rec under different confining stresses is 21.96~40.34, 46.73~66.84, 40.35~51.70, and 26.13~32.38, respectively, as shown in Figure 9. The calculated ranges of Rec are several times larger than those calculated by Zhang and Nemcik [111], in which the corresponding Rec = 3.5~4.5, 13.1~17.6, 19.3~24.8, and 6.3~8.6, respectively. The potential reason may be that, in the study of Zhang and Nemcik [111], the flow has entered the strong nonlinear regimes, where the factor A calculated using the fitted Forchheimer equation cannot be accurately estimated. That means, in such a situation, changing the magnitude of A will have negligible effect on the correlation coefficient between the tested results and fitted curves of Forchheimer equation, but it can robustly change the variations of Nf, resulting in the smaller ranges of Rec. This is the reason why we would like to use the variations ofT/T0 rather than Nf in the present study. Ranjith and Darlington [89] conducted nonlinear single-phase flow tests on real cylindrical rock joints with 110 mm height and 55 mm diameter under different confining pressures. The results depicted that, with increasing σn from 1.0 MPa to 5.0 MPa, Rec decreases from 31.29 to 28.77 (see Figure 10), which is in a smaller range compared with those in Figure 9. Chen et al. [112] presented the high-pressure packer test (HPPT) observations performed in fractured rock masses, which is located at about 450 m in depth in Qiongzhong County, China. Four test boreholes numbered ZK124, ZK126, ZK129-1, and ZK129-2, respectively, were obtained from in situ cites in the depths ranging from 33.4 m to 142.1 m with numbers of test intervals from 6 to 16, and the flow rate-pressure drop relationships as well as the magnitude of Rec were calculated. By collecting the tested results and replotting them in Figure 11, the variations of Rec with depth show that, for all boreholes, Rec fluctuates significantly in the range 25~66 and does not exhibit strong relationships with depth. The reasons may be that the boreholes at different depths contain fractures having different geometric properties such as aperture, persistence, orientation, and surface roughness. To simplify the model and clearly investigate the influence of normal stress and/or confining pressure on Rec, Zhou et al. [94] prepared four granite rock samples (G1~G4) taken from the depth 450~550 m in the Beishan area in Gansu Province, China, and two sandstone rock samples (S1~S2) taken from the depth of about 20 m in the Three Gorges Reservoir area, Hubei Province, China. The confining pressure varies from 1.0 to 30.0 MPa and Rec under each confining pressure is calculated. The results depicted in Figure 12 show that, with increasing confining pressure, Rec firstly increases and then decreases for the samples G2~G4 and S1~S2, because when the confining pressure is small, the decreased aperture due to the elastic compression will increase Rec [36]. When continuously increasing confining pressure, the contacts between two fracture walls are destroyed and the microfractures may propagate within rocks and the split rock fragments will be clamped by the two walls, which gives rise to complex flow behaviors and decreases Rec. For the sample G1, with increasing confining pressure, Rec steadily decreases and shows a different variation trend with other samples. They reported that Rec for the samples G1~G4 and S1~S2 varies in the ranges 0.075~9.243, 0.120~4.467, 0.160~5.108, 0.039~4.509, 0.191~4.090, and 0.026~2.976, respectively. The results of Rong et al. [100] as reported in Section 3.1 show that when us = 0 mm, Rec steadily decreases, presenting the same trend with the sample G1 in Figure 12. However, for samples under shear with us > 0 mm, with the increment of confining pressure, Rec firstly increases and then decreases as shown in Figure 13, presenting the same trends with most of the samples in Figure 12. For σn = 1.0~3.0 MPa, the range of Rec (Rec = 1.50~13.03) in the work of Rong et al. [100] is a little larger than that (Rec = 0.16~9.24) of Zhou et al. [94].Figure 8
Relationships betweenσn and un with different kn0.Figure 9
Relationships betweenT/T0 and Re based on the data replotted from Zhang and Nemcik ([111], Figure 8).
(a)
Sample number 1, JRC = 5.5 (b)
Sample number 2, JRC = 7.0 (c)
Sample number 3, JRC = 9.2 (d)
Sample number 2, JRC = 15.4Figure 10
Relationships betweenT/T0 and Re based on the data replotted from Ranjith and Darlington ([89], Figure 8).Figure 11
Relationships between depth and Re based on the data replotted from Chen et al. ([112], Tables 2–5).Figure 12
Relationships betweenRec and σn based on the data replotted from Zhou et al. ([94], Figure 8).Figure 13
Relationships betweenRec and σn based on the data replotted from Rong et al. ([100], Figure 9).


## 3.3. Effect of Fracture Surface Roughness

Fracture surface roughness can alter the distributions of void spaces of fractures, especially during shear [97, 114–117]. Besides, the rougher fracture surface contributes to the longer flow paths that a particle moves within fractures, resulting in a weaker conductivity/permeability if the same pressure difference is applied on the opposing boundaries [32, 73, 118, 119]. The previous works have shown that the fracture surface roughness, especially the secondary roughness, plays a significant role on the nonlinear flow properties of rock fractures, because the eddy flow occurs due to the surface roughness [77]. Wang et al. [120] established a series of 3D self-affine rock fractures using the successive random additions (SPA) method, in which the geometry of fracture surface is fractal [121, 122]. The fractal dimension (D3) is incorporated to characterize fracture surface roughness, which is correlated with the Hurst exponent (H) asD3 = 3-H [123–127]. Then a 3D Lattice Boltzmann method (LBM) that has been proven to be capable of solving the NS equations [128–131] is adopted to characterize the nonlinear flow regimes in rock fractures. The results show that, with increasing fracture surface roughness represented by D3, Rec decreases from 47.29 to 3.78 following an exponential function (see Figure 14), because the fluid flow within the rougher fractures will enter the nonlinear flow regime from the linear regime at a lower Re. The range of Rec (3.78~47.29) corresponding to Nf = 0.1 orT/T0 = 0.9 is larger than that (1.8~22.0) shown in the work of Wang et al. [120], since they utilized Nf = 0.05 orT/T0 = 0.95 as the threshold. Liu et al. [36] studied the transition of fluid flow from the linear to the nonlinear regime and quantified the critical hydraulic gradient (Jc) for flow in simple crossed fractures and complex fracture networks. The results show that Jc decreases significantly when JRC < 5 and then decreases slightly when JRC ≥ 5, following power law functions as presented in Figure 15. Note that, in their study, the onset of nonlinear flow is categorized using Nf = 0.01 orT/T0 = 0.99 rather than the commonly used Nf = 0.10 orT/T0 = 0.90 due to their high-precision flow testing system that allows them to adopt a much lower critical condition. By fitting the tested and simulated results, a multivariable regression function was proposed to calculate Jc, written as(35)Jc=λE-0.5·exp⁡300.8λE-0.01-0.03JRC-0.3Ni0.5-303.3,where Ni is the number of intersections and λ is the coefficient that quantifies the reduction of E due to the variation of fracture apertures in a DFN.Figure 14
Relationships between (a)T/T0 and Re and (b) Rec andD3 based on the data replotted from Wang et al. ([47], Figure 16).
(a)
Relationship betweenT/T0 and Re (b)
Relationship betweenRec andD3Figure 15
Relationships betweenJ and Re and JRC based on the data replotted from Liu et al. ([36], Figure 6(c)).The predictedJc using (35) fits well with the simulation results of DFNs with well-known geometric characteristics of fractures, which verifies the validity of (35).


## 3.4. Effect of Fracture Aperture

In the cubic law as described in (27), flow rate is linearly proportional to the cubic of fracture aperture, indicating that fracture aperture plays a vital role in controlling fluid flow in fractures. The fracture aperture has been observed to follow lognormal distributions [132–134] and incorporated into the theoretical models [28, 135–138]. Some previous works show that fracture aperture exhibits a power law relationship with fracture length with the exponent varying from 0.5 to 2.0 [62, 139–143]. The nonlinear flow properties of fractures also depend on the magnitude of fracture aperture. Liu et al. [113] took four rock samples of granite (S1 and S2), marble (S3), and limestone (S4) from the underground caverns of the Huangdao State Oil Reserves, China, and scanned their surface geometries using a three-dimensional laser-scanning rock surface instrument [144]. A total of twenty 2D rock fractures were prepared and fracture aperture varied from 1.0 mm to 10.0 mm. Fluid flow was simulated by directly solving the NS equations, and the nonlinear flow regimes were investigated. The results show that, with the increasing J, the reduction rate ofT/T0 increases and then decreases, existing an inflection point that has been utilized to quantify the transition from linear to nonlinear flow regime (see Figure 16(a)). With the increment of fracture aperture (E) from 1.0 mm to 10.0 mm, Jc decreases approximately four orders of magnitude following a power law function as shown in  (36) and Figure 16(b).(36)Jc=0.0627×E-3.87.Figure 16
Relationships between (a)T/T0 and J and (b) Jc and E based on the data replotted from Liu et al. ([113], Figure 4).
(a)
Relationship betweenT/T0 and J (b)
Relationship betweenJc and EJ c calculated corresponding toT/T0 = 0.9 is smaller than that calculated corresponding to the inflection point. Based on the simple DFNs established by Liu et al. [36], the relationships between Jc and E were established as shown in Figure 17. With the increment of E from 0.5 mm to 10.0 mm, Jc decreases and spans about five orders of magnitude following power law relationships, despite the magnitude of Ni. The coefficient and exponent of the fitted equations are close to those in  (36). The results show that, for both single fractures and fracture networks, Jc and E have power law correlations with a negative exponent.Figure 17
Relationships betweenJc and E based on the data replotted from Liu et al. ([36], Figure 6(a)).


## 3.5. Effect of Fracture Intersection

Fluid is redistributed at the intersections depending on the number and geometry of fracture segments connected to the intersection [74]. Wilson and Witherspoon [145] experimentally investigated the magnitude of laminar flow interference effects at sing fracture intersections. They found that when flow is in the laminar flow regime, the inertial effects at intersections can be negligible and when Re = 100, the head loss equals about five times the pipe diameters. Kosakowski and Berkowitz [146] established numerical models with a wide variability of fracture intersection geometry and studied the streamline distributions as well as the validity of Darcy’s law, by solving the NS equations using the available FEATFLOW package [147]. Their results indicate that the nonlinear inertial effects are significant for Re = 1~100, which is the range that often exists in the karst systems [14] and/or in the vicinity of wells during pump tests [15]. As introduced in the previous Sections 3.3 and 3.4, Liu et al. [36] estimated the influence of number of fracture intersections on the basis of a series of simple and complex DFNs. The results show that, with increasing Ni from 1 to 12, Jc decreases significantly and then slowly showing power law functions as shown in Figure 18. The reduction rate of Jc for JRC = 20 is about five times faster than that for JRC = 0. Comparisons among Figures 15, 17, and 18 indicate that Jc is the most sensitive to fracture aperture, followed by the number of intersections and fracture surface roughness.Figure 18
Relationships betweenJc and Ni based on the data replotted from Liu et al. ([36], Figure 6(b)).


## 4. Conclusions

Before permeability estimation of fractured rock masses, the critical conditions for the onset of nonlinear flow such as critical Reynolds number and critical hydraulic gradient should be firstly calculated, below which the cubic law in conjunction with modifications is applicable for giving reasonable solutions and beyond which the complex Navier-Stokes equations rather than cubic law should be solved. However there exist many difficulties of solving the Navier-Stokes equations in three-dimensional single fractures and/or fracture networks, because a set of coupled nonlinear partial derivatives of varying orders should be solved and discrete fracture network models having complex geometries should be established. As a result, there are still limited works focusing on nonlinear flow properties of fluid within fractures. This is the motivation of this work that summarizes the relative works on the evolution of critical conditions for the onset of nonlinear flow and provides a reference for those who are interested in this topic.The results show that, with the increment of shear displacement, critical Reynolds number increases significantly due to the robust effect of dilation and then increases slowly due to the weak effect of dilation. The variation trend of critical Reynolds number is close to that of fracture aperture during shear. When shear displacement varies from 0 mm to 20 mm, the critical Reynolds number increases from 0.001 to 25. The normal stress and/or confining pressure that is perpendicularly applied on the fracture would give rise to fracture aperture closure and affect the magnitude of critical Reynolds number. For most of the cases, with increasing normal stress and/or confining pressure from 0 to 30 MPa, critical Reynolds number increases firstly and then decreases, in the range 0.026~9.243. Fractures that have rougher surfaces will lose more energies due to eddy flow, resulting in the stronger nonlinear inertial effects and the smaller values of critical Reynolds number. As the fractal dimension of single three-dimensional fractures that is utilized to describe fracture surface roughness increases from 2.2 to 2.5, the critical Reynolds number decreases from 47.29 to 3.78 following an exponential function. Fracture aperture, which is a vital parameter to influence the transmissivity of fractures, plays an important role on critical hydraulic gradient. With increasing fracture aperture from 1.0 to 10.0 mm, the critical hydraulic gradient decreases for about 4~5 orders of magnitude following power law functions for both single fracture profiles and complex fracture networks. The flow paths of particles will be redistributed at fracture intersections; therefore, the fracture intersection would enhance the nonlinearity of fluid flow, resulting in a smaller critical hydraulic gradient with a larger number of intersections. The critical hydraulic gradient and number of intersections show exponential relationships. Among fracture aperture, surface roughness, and number of intersections, the critical hydraulic gradient is most sensitive to fracture aperture, followed by number of intersections and surface roughness.Future works should be devoted on the influences of coupled (thermo-hydro-mechanical-chemical) THMC processes, precipitation-dissolution, and biogeochemistry on the evolutions of critical conditions such as critical Reynolds number and critical hydraulic gradient.



---

*Source: 2176932-2017-07-06.xml*