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

1State Key Laboratory for Geomechanics and Deep Underground Engineering, China University of Mining and Technology, Xuzhou 221116, China 2State Key Laboratory of Mining Disaster Prevention and Control Co-Founded by Shandong Province and the Ministry of Science and Technology, Shandong University of Science and Technology, Qingdao 266590, China 3School of Engineering, Nagasaki University, Nagasaki 8528521, Japan


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 CO 2 sequestration, enhanced oil recovery, and geothermal energy development [1][2][3][4][5][6][7][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][10][11][12][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][15][16].In such case, using the cubic law to calculate fluid flow will overestimate the conductivity of rock fractures and/or fracture networks [17][18][19][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][25][26][27], aperture [28][29][30], surface roughness [31,32], dead-end [33], number of intersections [34,35], hydraulic gradient [36], boundary stress [37,38], anisotropy [39][40][41][42], scale [43][44][45][46], stiffness [47], coupled thermo-hydro-mechanical-chemical (HTMC) processes [48][49][50][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][54][55][56], although it cannot model the aperture heterogeneity of each fracture [57][58][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 () on the two opposing boundaries, such as  = 1 [53,[60][61][62][63][64],  = 0.1 [37],  = 0.001 [65,66], and  = unknown constants [10,30,42,[67][68][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.
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.

Governing Equations of Fluid Flow in Fractures
. .Navier-Stokes Equations.The flow of incompressible Newtonian fluid is governed by the NS equations, written in a tensor form as [75][76][77]  ( where   is the velocity of fluid in the -direction in the Cartesian coordinate system with  = , ,  and  = , , , respectively, the body force is   =   = [0, 0, −],  is the fluid density,  is the time,  is the pressure, and  is the dynamic viscosity. For steady-state flow, the parameters related to time  can be ignored and (1) gives Hydraulic head (ℎ) is defined as the summation of elevation height (  ) and hydraulic pressure head (/()), as follows: Hydraulic gradient () is defined as the ratio of hydraulic head difference to flow length, written as Substitution of (4) into (2) leads to a more abbreviated form of NS equations, as below [78]: . .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]: Equation ( 6) can be unfolded to the following forms: . .Reynolds Equation.Dimensional analysis is performed and some definitions are made: Λ is the dimension of fracture in -plane,  is the dimension of flow velocity in both and -directions,  is the dimension of flow velocity in direction, and  is the dimension of fracture aperture.
When the dimension of fracture aperture is much less than the dimension of fracture in the -plane ( ≪ Λ), 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: When the fracture aperture and the change in hydraulic head in the -direction are small, (11) yields Substitution of ( 12) into (9)∼ (11) gives where  1 ,  2 ,  1 , and  2 are the integration constants.Importing non-slip boundary conditions results in where  is the fracture aperture, and  = ±/2 are the lower and upper bounds of fracture void space in the -direction.Substitution of ( 14) into (13) leads to the values of integration constants ( 1 ,  2 ,  1 , and  2 ), and (13) can be rewritten as Integrations of ( 15) and ( 16) along the -direction leads to where   is the flow rate in the -direction, and   is the flow rate in the -direction.
The law of conservation of mass of flow rate is expressed as [80]  , +  , = 0.
Thus, substituting (18) into (19) results in the Reynolds equation, given as [17] . .Cubic Law.An assumption is made that fluid flows through two infinite smooth parallel plates with a constant aperture , as shown in Figure 2. Because fluid flows along the -coordinate, the flow velocities in the and -directions equal zero.The continuity equation of fluid is For incompressible fluid, (22) can be simplified to Substitution of ( 21) into (5) gives the governing equation of fluid flow in the x-direction, written as Substituting ( 23) into (24) gives rise to the fact that the left term equals zero.
Importing the non-slip boundary conditions (see ( 14)) and integrating (25) lead to where  is the flow rate through a fracture.Equation ( 26) is the expression for the cubic law, in which fracture width () is assigned a constant value of 1.When  ̸ = 1,  should be introduced in ( 26), as follows: . .Forchheimer Equation.With increasing flow velocity or Re, 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][83][84][85][86][87] where  and  are two coefficients that represent the linear and nonlinear terms of fluid flow, respectively.When  is very small, the nonlinear term ( 2 ) can be negligible and drops out.In this case, (28) reduces to where  is a constant that equals 12/ 3 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][89][90].
To describe the nonlinearity of fluid flow, a nonlinear factor (  ) is defined as denotes the ratio of pressure drop caused by the nonlinear term ( 2 ).In the previous works, it is assumed that the critical Reynolds number (Re  ) is the Re that corresponds to   = 0.1 [88,[91][92][93].
Besides the dimensionless Re that is used for characterizing the transition from linear to nonlinear flow regimes, the dimensionless Forchheimer number ( 0 ) is another widely accepted parameter, which is defined as the ratio of nonlinear to linear pressure losses, written as [16,90,94] Thus, (30) can be rewritten as The transmissivity (), which is a commonly used parameter to describe the conductivity of a fracture, is defined as When fluid flow is in the linear regime,  is a constant and is represented using  0 .Thus the ratio of  to  0 can be calculated according to Therefore, / 0 = 0.9 and   = 0.1 have the same physical meaning that the nonlinear term ( 2 ) contributes to 10% of the pressure drop, in which the current Re is Re  .

Nonlinear Flow Characteristics of Rock Fractures
. .Effect of Shear Displacement.Under a constant normal stress, increasing shear displacement (  ) 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   of 16 mm [95][96][97][98][99].This variation of fracture void space during shear can also change the nonlinear flow properties of fluid such as the critical Reynolds number Re  , 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, / 0 holds approximately constant value of 1 (see Figure 3).With continuously increasing Re, / 0 decreases.When / 0 = 0.9, Re  is calculated,  which is in the ranges 71.33∼340.01(  = 1∼5 mm), 75.93∼ 242.67 (  = 1∼6 mm), and 97.80∼326.07(  = 1∼5 mm), respectively.The range of Re  (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 / 0 = 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   , both the linear and nonlinear coefficients ( and ) 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 between T/T 0 and Re as shown in Figure 4 exhibit the same variation trend with those in Figure 3.With increasing   from 4 mm to 10 mm, Re  increases from 9.55 to 18.07, because the larger   gives rise to the larger fracture aperture and a smaller number of contacts, which results in the increment of Re  .Javadi et al. [90] experimentally investigated the role of shear processes on Re  for nonlinear flow through rough-walled fractures.The normal stress (  ) ranges from 1.0 to 5.0 MPa and   = 0∼ 20 mm.The results show that when   = 0∼1 mm, Re  has very small values approaching 0 for all cases (see Figure 5).When   = 1∼5 mm, Re  dramatically increases for   = 1.0 MPa and 5.0 MPa but gently increases for   = 3.0 MPa.For   = 5∼20 mm, the Re  holds approximately constant values for   = 5.0 MPa and slightly decreases for   = 1.0 MPa; however the Re  continuously increases for   = 3.0 MPa.The different variations of Re  ∼   curves may be induced by the statistical distributions of void space of fractures during shear.The range of Re  is from 0.001 to 25. Rong et al. [100] prepared 6 granite fractures with JRC = 6.67∼8.18and conducted coupled shear-flow tests with   = 0.5∼3.0MPa and   = 0.0∼ 10.9 mm.The results shown in Figure 6 depict that when   = 0∼0.5 mm, with the increment of   , Re  decreases, which is very different from those presented in Figure 5.With increasing   from 0.5 to 3.0 mm, Re  increases dramatically and then varies negligibly small when   exceeds 3.0 mm.This variation trend is similar to that (with   = 1.0 or 5.0 MPa) reported in Javadi et al. [90].As an example, when   = 1.0 MPa, Re  that corresponds to a small   in the work of Javadi et al. [90] is much smaller than that in the work of Rong et al. [100]; however when   is relatively large, Re  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 Re  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 Re  equals 12.56 for the experiment and 21.92 for the simulation, which is near the range of Re  reported in literature [90,95].
. .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][104][105][106][107][108].With increasing the normal stress   , both the fracture closure   and the increasing rate increase as shown in Figure 8, where  max is the maximum closure and  0 is the initial stiffness [109,110].This aperture variation/closure caused by normal stress contributes to the magnitudes of Re  and/or   .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  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.The results depicted in Figure 12 show that, with increasing confining pressure, Re  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 Re  [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 Re  .For the sample G1, with increasing confining pressure, Re  steadily decreases and shows a different variation trend with other samples.They reported that Re  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   = 0 mm, Re  steadily decreases, presenting the same trend with the sample G1 in Figure 12.However, for samples under shear with   > 0 mm, with the increment of confining pressure, Re  firstly increases and then decreases as shown in Figure 13, presenting the same trends with most of the samples in Figure 12.
. .Effect of Fracture Surface Roughness.Fracture surface roughness can alter the distributions of void spaces of fractures, especially during shear [97,[114][115][116][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 ( 3 ) is incorporated to characterize fracture surface roughness, which is correlated with the Hurst exponent (H) as D 3 = 3 −  [123][124][125][126][127]. Then a 3D Lattice Boltzmann method (LBM) that has been proven to be capable of solving the NS equations [128][129][130][131] is adopted to characterize the nonlinear flow regimes in rock fractures.The results show that, with increasing fracture surface roughness represented by  3 , Re  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 Re  (3.78∼47.29)corresponding to   = 0.1 or T/T 0 = 0.9 is larger than that (1.8∼22.0)shown in the work of Wang et al. [120],   16).
y = 0.040x −0.12 , R 2 = 0.9768 y = 0.036x −0.12 , R 2 = 0.7189 y = 0.037x −0.17 since they utilized   = 0.05 or T/T 0 = 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 (  ) for flow in simple crossed fractures and complex fracture networks.The results show that   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   = 0.01 or T/T 0 = 0.99 rather than the commonly used   = 0.10 or T/T 0 = 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   , written as where   is the number of intersections and  is the coefficient that quantifies the reduction of  due to the variation of fracture apertures in a DFN.The predicted   using (35) fits well with the simulation results of DFNs with well-known geometric characteristics of fractures, which verifies the validity of (35).
. .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][133][134] and incorporated into the theoretical models [28,[135][136][137][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][140][141][142][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 threedimensional 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 , the reduction rate of T/T 0 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 () from 1.0 mm to 10.0 mm,   decreases approximately four orders of magnitude following a power law function as shown in (36) and Figure 16(b).
= 0.0627 ×  −3.87 . ( calculated corresponding to T/T 0 = 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   and  were established as shown in Figure 17.With the increment of  from 0.5 mm to 10.0 mm,   decreases and spans about five orders of magnitude following power law relationships, despite the magnitude of   .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,   and  have power law correlations with a negative exponent.
. .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   from 1 to 12,   decreases significantly and then slowly showing power law functions as shown in Figure 18.The reduction rate of   for JRC = 20 is about five times faster than that for JRC = 0. Comparisons among Figures 15,17,and 18 indicate that   is the most sensitive to fracture aperture, followed by the number of intersections and fracture surface roughness.

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.

Figure 1 :
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 Science6 Core Collection (SCI-Expanded) over the past twenty years from 1997 to 2016.

Figure 2 :
Figure 2: Flow velocity distribution in a parallel plate model.
Re(a) Specimen #1 u s = 1 mm, fitting curve u s = 1 mm, experiment u s = 3 mm, fitting curve u s = 3 mm, experiment u s = 4 mm, fitting curve u s = 4 mm, experiment u s = 6 mm, fitting curve u s = 6 mm, experiment Re c Re (b) Specimen #2 u s = 1 mm, fitting curve u s = 1 mm, experiment u s = 2 mm, fitting curve u s = 2 mm, experiment u s = 3 mm, fitting curve u s = 3 mm, experiment u s = 4 mm, fitting curve u s = 4 mm, experiment u s = 5 mm, fitting curve u s = 5 mm, experiment Re c
Relationship between Re  and D 3

Figure 14 :
Figure 14: Relationships between (a) T/T 0 and Re and (b) Re  and D 3 based on the data replotted from Wang et al. ([47], Figure16).

Figure 15 :
Figure 15: Relationships between  and Re and JRC based on the data replotted from Liu et al. ([36], Figure 6(c)).

Figure 16 :
Figure 16: Relationships between (a) T/T 0 and  and (b)   and  based on the data replotted from Liu et al. ([113], Figure 4).
Figure 8: Relationships between   and   with different  0 .
[89]strong nonlinear regimes, where the factor  calculated using the fitted Forchheimer equation cannot be accurately estimated.That means, in such a situation, changing the magnitude of  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   , resulting in the smaller ranges of Re  .This is the reason why we would like to use the variations of T/T 0 rather than   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   from 1.0 MPa to 5.0 MPa, Re  decreases from 31.29 to 28.77 (see Figure10), which is in a smaller range compared with those in Figure9.Chen et al.  under each confining pressure is calculated.