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.

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 [18]. 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 [913]. 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, 1416]. In such case, using the cubic law to calculate fluid flow will overestimate the conductivity of rock fractures and/or fracture networks [1720]. 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 [2427], aperture [2830], surface roughness [31, 32], dead-end [33], number of intersections [34, 35], hydraulic gradient [36], boundary stress [37, 38], anisotropy [3942], scale [4346], stiffness [47], coupled thermo-hydro-mechanical-chemical (HTMC) processes [4851], 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 [5356], although it cannot model the aperture heterogeneity of each fracture [5759]. 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 [53, 6064], [37], [65, 66], and = unknown constants [10, 30, 42, 6769]. 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.

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 [7577]where is the velocity of fluid in the -direction in the Cartesian coordinate system with and , respectively, the body force is , 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]:

2.2. Stokes Equations

When the flow rate/velocity is very small, which corresponds to a small Reynolds number (), 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:

2.3. 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) giveswhere , and are the integration constants.

Importing non-slip boundary conditions results inwhere is the fracture aperture, and 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 (, and ), and (13) can be rewritten as

Integrations of (15) and (16) along the -direction leads towhere 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]

Thus, substituting (18) into (19) results in the Reynolds equation, given as [17]

2.4. 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 towhere 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 should be introduced in (26), as follows:

2.5. Forchheimer Equation

With increasing flow velocity or , 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 [8287]where and are two coefficients that represent the linear and nonlinear terms of fluid flow, respectively. When is very small, the nonlinear term () can be negligible and drops out. In this case, (28) reduces towhere is a constant that equals 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 . 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 [8890].

To describe the nonlinearity of fluid flow, a nonlinear factor () is defined as

denotes the ratio of pressure drop caused by the nonlinear term (). In the previous works, it is assumed that the critical Reynolds number () is the that corresponds to = 0.1 [88, 9193].

Besides the dimensionless that is used for characterizing the transition from linear to nonlinear flow regimes, the dimensionless Forchheimer number () 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 . Thus the ratio of to can be calculated according to

Therefore, = 0.9 and = 0.1 have the same physical meaning that the nonlinear term () contributes to 10% of the pressure drop, in which the current is .

3. Nonlinear Flow Characteristics of Rock Fractures

3.1. 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 [9599]. This variation of fracture void space during shear can also change the nonlinear flow properties of fluid such as the critical Reynolds number , 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, holds approximately constant value of 1 (see Figure 3). With continuously increasing Re, decreases. When , 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 (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.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/T0 and as shown in Figure 4 exhibit the same variation trend with those in Figure 3. With increasing from 4 mm to 10 mm, 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 . Javadi et al. [90] experimentally investigated the role of shear processes on 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, has very small values approaching 0 for all cases (see Figure 5). When = 1~5 mm, dramatically increases for = 1.0 MPa and 5.0 MPa but gently increases for = 3.0 MPa. For = 5~20 mm, the holds approximately constant values for = 5.0 MPa and slightly decreases for σn = 1.0 MPa; however the continuously increases for = 3.0 MPa. The different variations of ~ curves may be induced by the statistical distributions of void space of fractures during shear. The range of 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 = 0.5~3.0 MPa and = 0.0~10.9 mm. The results shown in Figure 6 depict that when = 0~0.5 mm, with the increment of , decreases, which is very different from those presented in Figure 5. With increasing from 0.5 to 3.0 mm, 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, 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, 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 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 equals 12.56 for the experiment and 21.92 for the simulation, which is near the range of reported in literature [90, 95].

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 [103108]. With increasing the normal stress , both the fracture closure and the increasing rate increase as shown in Figure 8, where is the maximum closure and is the initial stiffness [109, 110]. This aperture variation/closure caused by normal stress contributes to the magnitudes of 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.4, the range of 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 are several times larger than those calculated by Zhang and Nemcik [111], in which the corresponding = 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 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 . This is the reason why we would like to use the variations of T/T0 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, 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 were calculated. By collecting the tested results and replotting them in Figure 11, the variations of with depth show that, for all boreholes, 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 , 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 under each confining pressure is calculated. The results depicted in Figure 12 show that, with increasing confining pressure, 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 [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 . For the sample G1, with increasing confining pressure, steadily decreases and shows a different variation trend with other samples. They reported that 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, 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, firstly increases and then decreases as shown in Figure 13, presenting the same trends with most of the samples in Figure 12. For = 1.0~3.0 MPa, the range of ( = 1.50~13.03) in the work of Rong et al. [100] is a little larger than that ( = 0.16~9.24) of Zhou et al. [94].

3.3. Effect of Fracture Surface Roughness

Fracture surface roughness can alter the distributions of void spaces of fractures, especially during shear [97, 114117]. 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 () is incorporated to characterize fracture surface roughness, which is correlated with the Hurst exponent (H) as D3 = [123127]. Then a 3D Lattice Boltzmann method (LBM) that has been proven to be capable of solving the NS equations [128131] is adopted to characterize the nonlinear flow regimes in rock fractures. The results show that, with increasing fracture surface roughness represented by , 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 (3.78~47.29) corresponding to = 0.1 or T/T0 = 0.9 is larger than that (1.8~22.0) shown in the work of Wang et al. [120], since they utilized = 0.05 or T/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 () 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/T0 = 0.99 rather than the commonly used = 0.10 or T/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 , written aswhere 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).

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 [132134] and incorporated into the theoretical models [28, 135138]. 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, 139143]. 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 , the reduction rate of T/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 () 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).

calculated corresponding to T/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 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.

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 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.

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.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

This study has been partially funded by National Basic Research 973 Program of China, China (Grant no. 2013CB036003), and National Natural Science Foundation of China, China (Grant nos. 51579239, 51323004). These supports are gratefully acknowledged.