Experimental and numerical investigation on nonlinear flow behaviour through three dimensional fracture intersections and fracture networks

Computers and Geotechnics 121 (2020) 103446

Contents lists available at ScienceDirect

Computers and Geotechnics journal homepage: www.elsevier.com/locate/compgeo

Research Paper

Experimental and numerical investigation on nonlinear flow behaviour through three dimensional fracture intersections and fracture networks Feng Xionga,b,c, Wei Weia,b, Chaoshui Xuc, , Qinghui Jianga,b, ⁎

T

⁎

a

School of Civil Engineering, Wuhan University, Wuhan 430072, PR China Key Laboratory Safety Geotechnical & Structural Engineering Hubei Province, Wuhan 430072, PR China c School of Civil, Environmental and Mining Engineering, University of Adelaide, Adelaide 5005, Australia b

ARTICLE INFO

ABSTRACT

Keywords: Nonlinear flow through rock fracture networks Fracture intersection Fracture intersection angle Critical hydraulic gradient Fluid mechanics

Fluid flow tests were conducted on four kinds of fracture intersections, namely straight fracture intersection, buckling fracture intersection, crossing fracture intersection and furcating fracture intersection. Corresponding numerical simulations were performed by solving the Navier-Stokes equations. The resultant flow force perpendicular to fracture walls were derived. The results show that at high hydraulic gradients the increased normal resultant force perpendicular to fracture walls will change the original flow direction and consequently cause the formation of eddies which triggers flow nonlinearity. Greater buckling angle and smaller crossing angle would result in the onset of nonlinear flow at a higher critical hydraulic gradient for the cases of buckling and crossing fracture intersections. The critical hydraulic gradient for the case of furcating fracture intersection depends on the absolute angle between furcating flow and the original flow direction. Based on the analyses, three nonlinear models for coefficient B in the Forchheimer law corresponding to three cases of fracture intersections are proposed, which are related to the fracture intersection angle and hydraulic aperture. These models are further validated using flow experimental data from a real fracture network created in a rock specimen with the corresponding fracture model constructed using the computed tomography (CT) method.

1. Introduction Due to geological processes, weathering, water invasion and tectonic actions [1,16,22,27], a complex discrete fracture network (DFN) system is commonly formed within rock masses. These fractures are generally characterized as randomly distributed, inhomogeneous and morphologically diverse. The permeability of most intact rocks is usually very small, and therefore fractures within the rock mass are the main conduit channels for fluid flows. Many published studies (e.g., [31,10,11,7,33,32,30,28]) used DFN to investigate fluid flow behaviour in fractured rock masses under different application contexts. However, issues related to fracture intersections are generally ignored or oversimplified. As fracture intersections are the main components in a DFN, understanding the flow behaviours around these intersections is critical to the generation of more realistic flow models for DFNs, particularly for the nonlinear flow regime. Many flow tests in fractures [5,39,13,3,23,18] show that the relationship between the flow rate and pressure gradient becomes nonlinear at high Reynolds number and the cubic law is no longer applicable. The mathematical description of nonlinear flow through fractures

⁎

widely adopted is the Forchheimer’s law [38,17], which consists of a linear term and a nonlinear term. Even though the Forchheimer’s law was empirically established based on experimental results, its theoretical support was established in Wang et al. [26], where the ratio of aperture variation to length was used as the perturbation parameter to derive the perturbation solution identical to the Forchheimer’s law. The linear term in the Forchheimer’s law represents the intrinsic fracture permeability [4,2], while the nonlinear term represents the component triggered by inertial effect due to complicated rock fracture surfaces [2,24] and contact areas [29]. Sævik and Nixon [21] used node classification scheme to characterize the hydraulic connectivity of DFNs for the estimation of effective permeability. They showed that there were four basic types of fracture intersections in DFNs, namely, I-node, Y-node, V-node and Xnode. These fracture intersections are referred to as straight fracture intersection (SFI), buckling fracture intersection (BFI), crossing fracture intersection (CFI) and furcating fracture intersection (FFI) in this paper, as shown in Fig. 1. When fluid is injected into a fracture intersection, the flow pressure and flowrate distribution are affected by many factors such as fracture intersection angle, fracture surface roughness and

Corresponding authors at: Wuhan University, University of Adelaide, Australia. E-mail addresses: [email protected] (C. Xu), [email protected] (Q. Jiang).

https://doi.org/10.1016/j.compgeo.2020.103446 Received 26 June 2019; Received in revised form 7 January 2020; Accepted 7 January 2020 0266-352X/ © 2020 Elsevier Ltd. All rights reserved.

Computers and Geotechnics 121 (2020) 103446

F. Xiong, et al.

Nomenclature

θ crossing angle γ furcating angle g gravity acceleration critical hydraulic gradient Jc e hydraulic aperture head loss hLocal φ exterior angle λ1, λ2 coefficients of head loss equation v fluid velocity m1,m2, m3,m4, m5, m6 coefficients of proposed B model L flow path length N number of fracture intersection v1, v2 mean velocity vector F external force of fracture drag force of fracture Fτ Fn normal force of fracture Fx, Fy x- and y-direction external force Fnx, Fny x- and y-direction normal force Fτx, Fτy x- and y-direction drag force Abs () absolute value function

List of symbols Ρ P μ u b w Q ▽P Re J A B Ta T0 β E α

fluid density fluid pressure fluid viscosity coefficient velocity vector mechanical aperture fracture width flow rate pressure gradient Reynolds number hydraulic gradient linear coefficient of Forchheimer law nonlinear coefficient of Forchheimer law apparent transmissivity apparent transmissivity value in Darcy’s flow regime Forchheimer coefficient non-Darcy effect factor buckling angle

Fig. 1. . A schematic diagram of a discrete fracture network containing different types of fracture intersections: (b) Straight fracture intersection (SFI), (c) Bucking fracture intersection (BFI), (d) Crossing fracture intersection (CFI) and (e) Furcating fracture intersection (FFI).

pressure gradient. The fluid flow behaviour derived for a single fracture is not applicable to the case of fracture intersections, especially when the nonlinear flow regime is considered. Numerical simulations were performed by Liu et al. [14,15] to estimate the effects of fracture intersections on the nonlinear flow behaviours in crossing fractures. It was found that when the hydraulic gradient J is greater than 10–4, a nonlinear flow regime occurs in which the inertial effects caused by the fracture intersections should not be neglected. Li et al. [12] further conducted flow tests and numerical simulations on crossed fractures and the results showed that larger apertures, rougher fracture surfaces and a greater number of intersections in a DFN would result in the onset of nonlinear flow at a lower critical hydraulic gradient. Although many studies described were conducted experimentally and numerically, the influence mechanisms of fracture intersection angle on fluid flow behaviour were rarely discussed. In addition, the corresponding nonlinear flow model was also rarely reported. In this paper, four groups of three dimensional fracture intersections, SFI, BFI, CFI and FFI, were studied with particular focus on the nonlinear flow

behaviour. Flow tests were carried out on these fracture models and the corresponding numerical simulations were carried out by solving the Naiver-Stokes equations to characterize the nonlinear relationship and flow distribution behaviour. Semi-empirical expression for the nonlinear term parameter in the Forchheimer law were proposed for different cases of fracture intersections. The proposed models were validated using flow experimental data from a real rock fracture network created in a cylindrical rock specimen with the corresponding fracture model constructed using the computed tomography (CT) technique. 2. Theoretical background The incompressible steady state fluid flow is governed by the Naiver-Stokes equations, written as [40]: (1)

u =0

u u +

P =µ

2u

(2)

where ρ, P, μ, u denoting the density, pressure, viscosity coefficient and 2

Computers and Geotechnics 121 (2020) 103446

F. Xiong, et al.

velocity of fluid, respectively. For low velocity fluid flow within a fracture, the Naiver-Stokes equation can be simplified to the cubic law:

P=

12µQ wb3

E=

Q µw

Jc =

(4)

3.1. Specimen preparation Flow tests were conducted using granite specimens to help understand the behaviour of flow at rock fracture intersections. The granite specimens were sourced from the site of Dagangshan Hydropower Project, China. Cylindrical specimens were prepared with the diameter of 50 mm and lengths of 50 mm and 100 mm respectively. Fracture intersections were prepared carefully within the specimen using a diamond wire saw (Fig. 2). The maximum cut depth and length of the machine are 110 mm and 200 mm, respectively. The radius of diamond wire is 0.1 mm and the travelling speed of the wire is 5.0 m/s. Four different groups of fracture intersections were prepared as shown in Fig. 3. The SFI group includes fractures (S-1 and S-2) with 80° and 90° dip angle. The BFI group includes fractures (B-1, B-2 and B-3) with 90°, 120° and 150° buckling angles (α). The CFI group includes fractures (C1, C-2 and C-3) with 5°, 30° and 60° crossing angles (θ). The FFI group includes fractures (F-1, F-2 and F-3) with 20°, 40° and 75° furcating angles (γ). In these fractures, γ1 is 20° and γ2 is 0° for F-1. γ1 is 20° and γ2 is 20° for F-2. γ1 is 20° and γ2 is 55° for F-3 (shown in Fig. 3d). The length of the specimen for C-3 and F-3 is 50 mm, and that for others is 100 mm.

(5)

1 1+

Re

(8)

3. Flow tests

where A and B represent the linear and nonlinear coefficients. Coefficient A can be determined by cubic law, i.e., A = 12μ/we3ρg. g is the acceleration of gravity. However, no model is currently available to the evaluation of B for intersecting fractures. In Zimmerman et al. [38], a normalized transmissivity (Ta/T0) was determined as:

Ta/ T0 =

A2 E A2 E 2 + B (1 E ) B (1 E )2

where E = 0.1.

Re is commonly used to predict the onset of nonlinear flow in a single fracture or a fracture intersection. However, the prediction is not directly applicable for DFNs where many fractures and fracture intersections are involved. For a DFN, the hydraulic gradient J of interest for the onset of the nonlinear flow is that defined over the region, not over an individual fracture. The nonlinear relationships were revealed in Zimmerman et al. [38] when the value of Reynolds number (Re) is greater than 10 for roughwalled fractures. In such conditions, the Forchheimer law [38,2] can be applied to describe this relationship. In the present study, the hydraulic gradient linearly proportional to ▽P/ρg is represented by J, the Forchheimer law can be rewritten as:

J = AQ + BQ2

(7)

This factor E represents the proportional contribution of the nonlinear term to the total pressure drop. In general, E = 0.1 was used as the threshold to quantify the transition of the flow regime from linear to nonlinear, see for examples, Zhang et al. [41], Zou et al. [40] and Zhou et al. [37]. Substituting Eq. (7) into Eq. (5) yields the estimation for the critical hydraulic gradient Jc:

(3)

where b is the fracture aperture, w is the fracture width, Q is the flow rate. ▽P is the pressure gradient along flow direction. For flow in fracture, the Reynolds number, Re, is defined as [38,9]:

Re =

BQ2 AQ + BQ 2

(6)

where Ta is the apparent transmissivity, T0 is a special apparent transmissivity in Darcy’s flow regime. β is the Forchheimer coefficient used to describe the nonlinear behavior of flow. There are significant research efforts to estimate the onset point of flow rate where the deviation from linear flow occurs. A Non-Darcy effect factor E [35] was introduced to quantify the nonlinearity of fluid flow, which is defined as:

3.2. Experimental setup Fig. 4 shows a schematic view of the experimental system used in our tests, which is specifically manufactured. The system consists of a

Fig. 2. . The diamond wire saw used for specimen preparation. 3

Computers and Geotechnics 121 (2020) 103446

F. Xiong, et al.

Fig. 3. . Specimens of four types of fracture intersections. Specimens were labelled according to the intersection: (a) SFI (S-1, S-2). (b) BFI (B-1, B-2 and B-3). (c) CFI (C-1, C-2 and C-3). (d) FFI (F-1, F-2 and F-3).

4

Computers and Geotechnics 121 (2020) 103446

F. Xiong, et al.

Fig. 3. (continued)

pressure pump with the pressure range of 0–3.0 MPa, a pressure gauge with a measurement accuracy of ± 1.0 Pa, plus a flow cell and a flow rate measurement apparatus. The flow rate is measured by recording the height changes of water column in the flowmeter.

4. Numerical simulation To understand better the flow behaviour at fracture intersections, numerical models were created based on the geometry of tested specimens. Note fracture surface roughness could have significant impact on the non-linear behaviour of flow within rock fractures [25]. However, the study of this impact is not the focus of this work. Fractures in our specimens were created using a diamond wire saw and therefore the resulted fracture surfaces were relatively smooth. Therefore, the influence of fracture surface roughness was not considered in our experimental analyses and numerical models. Note fractures in these samples were created using a diamond saw. The cut surfaces are not perfectly smooth and therefore the two walls of the fracture will not match perfectly when they come into contact under confining pressure. Void spaces exist between the walls that provide conduits for the fluid to flow through. As the influence of fracture surface roughness on flow characteristics is not the focus of this work, the fracture aperture discussed hereafter is the equivalent hydraulic aperture which is backcalculated from experimental measurements (J-Q curves). The aperture of each fracture branch for fracture intersection model could not be measured separately in the present flow test. Therefore, the fracture

3.3. Test procedure The specimen enclosed in a rubber jacket was put in the flow cell and a small confining pressure (1.0 MPa) was applied to the specimen to ensure there is no gap between the specimen and the rubber jacket to prevent fluid leakage. The small confining pressure also ensure fracture walls are always in contact during the flow tests. A flow run was conducted to ensure that no air enters the pipes and the cell. Water was then introduced into the fracture intersection through the inlet and collected at the outlet (Fig. 4). A series of flow pressures at the inlet in the range of 0–1 MPa were tested and their corresponding flowrates were recorded. The pressure was transformed into hydraulic gradient using the equation J = P/ρgL and therefore the relationships between flowrate Q and hydraulic gradient J can be obtained, where L is flow path length. In the tests, the temperature was maintained at 20 °C so the density and viscosity of water are 997.298 kg/m3 and 0.00108 Pa·s, respectively.

Fig. 4. . Schematic diagram of the fracture flow testing system. 5

Computers and Geotechnics 121 (2020) 103446

F. Xiong, et al.

Fig. 5. . Flow diagram in (a) BFI, (b) CFI and (c) FFI.

Q2 =

wb3 P3 P2 12µ L 2

(10)

Considering the flow equilibrium, Q1 equals to Q2. Hence, a function of flow rate and pressure can be obtained:

Q = Q1 = Q2 =

wb3 P1 P2 12µ (L1 + L 2)

(11)

Likewise, two equations can be induced for CFI (Fig. 5b) and FFI (Fig. 5c), respectively:

branches were assumed to have the same apertures. The same assumption was also done by Ishibashi et al. [8]. For example, according to cubic law, there are two equations existing in BFI (shown in Fig. 5a):

wb3 P1 P3 12µ L1

wb3 2(P1 P2) 12µ (L1 + L2 )

Q =Q2 + Q3 =

wb3 12µ

P1 P2 L2 L3 + L1 L3 + L1 L2 L2 + L3

(12)

(13)

Hence, flow rate and pressure data, original from flow experiments, could be taken into Eqs. (11)–(13) to back calculate the hydraulic aperture value, which is designated as the aperture of numerical models. Even though this setup simplifies the fracture surface topography, the established numerical model is same with experimental model generally. The intersecting fractures were meshed using ANSYS code. A structured hexahedral mesh was used for fractures and fracture intersections. The mesh at the intersection and the boundary layers were refined until stable accurate results were obtained. This in general

Fig. 6. . An example of meshing used in numerical simulations.

Q1 =

Q =Q1 + Q2 =

(9)

6

Computers and Geotechnics 121 (2020) 103446

F. Xiong, et al.

Fig. 7. . Comparison of experimental and numerical results for four groups of fractures (E: experimental results, N: Numerical simulation results).

resulted in the total number of elements in each model greater than 104. A mesh example is showed in Fig. 6 for the case of a crossing fracture intersection. The Viscous-Laminar module was selected for water flow simulation. The constant velocity condition was used for the inlet boundary of the fracture intersection, and the pressure boundary condition was used for the outlet. Non-slip boundary conditions were assigned to the fracture walls. In numerical simulations, the momentum and pressure were solved using a 2nd-order upwind scheme and a 2nd-order scheme, respectively. The pressure-velocity coupling was calculated using the SIMPLE algorithm. The solution gradients were evaluated by direct interpolation with a least square method at the center of each cell. Steady state solutions were considered to be converged when the normalized velocity and continuity residual ratios were less than 1.0e−4.

5. Results from the flow experiments and extended numerical simulations 5.1. Nonlinear flow behaviour at fracture intersections The experimental results of Q-J curves are shown in Fig. 7 for the cases of SFI, BFI, CFI and FFI, respectively. These results suggest that hydraulic gradient J increases linearly initially, and then nonlinearly as flowrate Q increases, except for the case of SFI, where the relationship remains linear. The Q-J data in linear stage is used to calculate hydraulic apertures, which are listed in Table 1. The aperture of the corresponding numerical model uses the same value of hydraulic aperture obtained by experimental data. Simulation results are also plotted in Fig. 7, which show a good agreement with experimental results (5–10% discrepancy) and therefore numerical simulations can be

Table 1 The values of A, B and Jc for the Forchheimer law derived based on experimental results. SFI

A(s/m3)

B(s2/m6)

e/mm

β

Jc

CFI

A(s/m3)

B(s2/m6)

e/mm

β

Jc

S_1 S_2

2.14e4 1.94e4

/ /

1.04 1.07

/ /

/ /

C_1 C_2 C_3

4.00e4 9.64e4 1.42e5

4.20e7 2.30e9 1.80e10

0.84 0.66 0.55

5.3e−5 1.2e−3 6.3e−3

4.700 0.500 0.138

BFI B_1 B_2 B_3

A(s/m3) 3.48e4 3.82e4 3.36e4

B(s2/m6) 4.37e9 1.27e9 1.80e8

e/mm 0.88 0.86 0.67

β 0.0062 0.0017 0.00027

Jc 0.034 0.142 0.774

FFI F_1 F_2 F_3

A(s/m3) 9.30e4 9.45e4 8.20e4

B(s2/m6) 1.41e10 3.1e9 7.30e10

e/mm 0.66 0.63 0.57

β 0.008 0.002 0.045

Jc 0.075 0.356 0.011

7

Computers and Geotechnics 121 (2020) 103446

F. Xiong, et al.

intersections. When flowrate is small (Fig. 7), the inertial force is much less than the viscous force and J linearly increases with increasing Q. In this condition, the nonlinear term BQ2 could be ignored and the cubic law could describe adequately the relationship between Q and J. However, when the flowrate increases to the critical value (i.e. Q = 1e−6 m3/s as shown in Fig. 7), the impact of inertial effects increases and the cubic law is no longer applicable and therefore the Forchheimer law could be employed. The nonlinearity and the deviations from the cubic law are enhanced as flowrate increases. As described in [36], there are three flow stages for fracture intersections: viscous regime, weak inertia regime and strong inertia regime. According to the definition of critical hydraulic gradient Jc (Eq. (8)) and Forchheimer coefficient β, those values are listed in Table 1. The Jc is range from 0.011 to 4.7, and β from 5.3e−5 to 0.045. 5.2. Flow patterns Fig. 7 and Table 1 does not give a clear indication about the flow characteristics at fracture intersections and their influences on the nonlinear flow behavior. For these analyses, numerical simulations are used. The velocity profile for the case of SFIs is shown in Fig. 8. Clearly all flow velocity vectors are parallel to fracture walls in this case. No nonlinear flow features can be observed for the hydraulic gradients examined. Therefore, the flow behaviour in SFI can be described by the cubic law for the hydraulic gradients tested. Fig. 9 shows the pressure and flow profiles for the case of BFI at

Fig. 8. . Flow velocity profiles for SFI (S_1 and S_2).

used to further analyze the nonlinear flow behaviour at fracture intersections. The quadratic polynomial regression in the form of Eq. (5) is used to fit the Q-J curves in the cases of BFI, CFI and FFI. The high coefficients of determination of R2 (> 0.95) demonstrate that the Forchheimer law adequately describes the nonlinear flow behaviour at fracture

Fig. 9. . Flow velocity profiles for BFI (B_1, B_2 and B_3) at different hydraulic gradients.

8

Computers and Geotechnics 121 (2020) 103446

F. Xiong, et al.

Fig. 10. . Flow velocity profiles for CFI (C_1, C_2 and C_3) at different hydraulic gradients.

different hydraulic gradients. For the buckling angles α of 90° and 120°, as the hydraulic gradient increases, eddies appear at the buckling position of top surface and lower part of the bottom surface, causing the effective flow channel to be narrowed. However, for the buckling angle of 150°, no eddy can be observed in the fracture, indicating a very low degree of nonlinearity. Fig. 10 shows the pressure and flow profiles for the case of CFI at J = 0.01 and J = 0.3. For crossing angle θ of 5°, the flow streamlines follow closely the curvatures of the fracture walls. For crossing angles θ of 30° and 60°, two eddies emerge in lower parts of the fracture which cause the reduction in the effective area available for flow when hydraulic gradient J is 0.3. This result is consistent with that of Liu et al. [15]. Fig. 11 shows the pressure and flow distribution profiles for the case of FFI at J = 0.01 and J = 0.3. A different nonlinear flow behavior can be observed, even at a high hydraulic gradient where the flow is in the Forchheimer regime. Two eddies appear in the lower parts of the fracture with the furcating angle γ of 20° and 75°. But for the furcating angles of 40°, no eddy can be observed within the fracture even when the flow is in the Forchheimer regime

Therefore, the formation of eddies is related to the fracture intersecting angle, which demonstrate that the intersecting angle α, θ and γ play a significant role in the nonlinear behavior of flow in intersecting fractures. 5.3. Onset of nonlinear flow To systematically investigate the nonlinear flow behavior in single fracture intersections, many extended numerical models beyond the flow experimental cases were established. Five BFI models have been established with buckling angle α 115°, 130°, 140°, 160° and 170°. Six CFI models have been established with crossing angle β 10°, 20°, 40°, 50°, 70° and 80°. Six FFI models have been established with furcating angle γ1 = 20° and γ2 = 10°, γ1 = 20°, γ2 = 50°; γ1 = 10°, γ2 = 55°; γ1 = 5°, γ2 = 60°; γ1 = 5°, γ2 = 55° and γ1 = 10°, γ2 = 50°. The numerical simulations by solving the Naiver-Stokes equations were also performed by ANSYS FLUENT. The simulation results are listed in Table 2. In order to understand better the influence of fracture intersections on the nonlinear flow behaviour, the Forchheimer coefficient β was

9

Computers and Geotechnics 121 (2020) 103446

F. Xiong, et al.

Fig. 11. . Flow velocity profiles for FFI (F_1, F_2 and F_3) at different hydraulic gradients. Table 2 The values of A, B and Jc for the Forchheimer law of extended simulations. BFI

A(s/m3)

B(s2/m6)

e/mm

β

Jc

CFI

α = 115° α = 130° α = 140° α = 160° α = 170° FFI γ1 = 20° γ2 = 10° γ1 = 20° γ2 = 50° γ1 = 10° γ2 = 55° γ1 = 5° γ2 = 60° γ1 = 5° γ2 = 55° γ1 = 10° γ2 = 50°

4.69e4 4.69e4 4.69e4 4.69e4 4.69e4 A(s/m3) 1.11e5 1.11e5 1.11e5 1.11e5 1.11e5 1.11e5

2.10e9 8.50e8 3.90e8 3.90e7 5.80e6 B(s2/m6) 6.40e9 5.70e10 1.20e11 1.80e11 1.50e11 9.90e10

0.80 0.80 0.80 0.80 0.80 e/mm 0.60 0.60 0.60 0.60 0.60 0.60

2.2e−3 9.0e−4 4e−4 4.2e−5 6.2e−6 β 2.88e−3 2.57e−2 5.40e−2 8.10e−2 6.75e−2 4.46e−2

1.29e−1 3.19e−1 6.96e−1 6.96 4.68e1 Jc 2.38e−1 2.67e−2 1.27e−2 8.47e−3 1.02e−2 1.54e−2

θ θ θ θ θ θ

calculated and the values for different cases are listed in Table 2. The relationships between β and buckling angle α in BFIs, crossing angle θ in CFIs and Furcating angle γ in FFIs are plotted in Fig. 12 . As expected, β decreases as α increases for BFIs, but β increases as θ increases for CFIs

= = = = = =

10° 20° 40° 50° 70° 80°

A(s/m3)

B(s2/m6)

e/mm

β

Jc

4.69e4 4.69e4 4.69e4 4.69e4 4.69e4 4.69e4

1.80e8 7.40e8 2.80e9 4.30e9 1.60e10 2.0e10

0.80 0.80 0.80 0.80 0.80 0.80

1.92e−4 7.89e−4 2.99e−3 4.59e−3 1.71e−2 2.13e−2

1.51 3.67e−1 9.69e−2 6.31e−2 1.70e−2 1.36e−2

or as |γ1-γ2| increases for FFIs. All these relationships are non-linear. The values of Jc for all the cases can be calculated using Eq. (8), which are also listed in Table 2. The relationships between Jc and α for BFIs, θ for CFIs and γ for FFIs are displayed in Fig. 13. As can be seen,

10

Computers and Geotechnics 121 (2020) 103446

F. Xiong, et al.

Fig. 12. . The relationships of (a) β − α, (b) β − θ and (c) β − γ based on experiments and simulations (blue points: experimental data, black points: numerical results). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

the relationships are also non-linear. As α increases, Jc increases for BFIs. But as θ (CFIs) or |γ1-γ2| (FFIs) increases, Jc decreases, as expected (see Fig. 13). Hence, it is also proved that the intersecting angle has significant influences on nonlinear flow behavior of single fracture intersections.

flow phenomenon would occur at high flow rate. In addition, Fnx and Fny are functions of the buckling angle α. As the buckling angle α increases, Fnx decreases while Fny increases. Therefore, the forces that can potentially change the flow direction become lower and the flow nonlinearity thus becomes weaker. This trend is consistent with experimental and numerical results. From Eq. (a.15) and Eq. (a.16), Fnx for the case of CFIs increases with flow velocity. Higher flow velocity results in larger eddies being formed in the lower parts of the crossing fracture. In addition, Fnx also increases with the crossing angle θ. Greater crossing angle θ leads to stronger nonlinear flow, which is consistent with experimental and numerical results. Similarly, from Eq. (a.23) and Eq. (a.24), Fny for the case of FFIs increases with flow velocity. Higher velocity causes the formation of larger eddies in lower parts of the fracture, indicating increasing nonlinear flow behaviour. The flow nonlinearity in this case also depends on the value of furcating angle γ (γ1 and γ2 in Fig. a3). Greater abs(γ1γ2) will lead to greater Fny hence stronger nonlinear flow.

5.4. Flow mechanics analysis in fracture intersections There is no nonlinear flow behaviour observed for the case of straight fracture for the hydraulic gradients examined. But strong flow nonlinearity is observed for the cases of BFI, CFI and FFI. Therefore, three different models of flow physics were established to analyze the nonlinear flow mechanics for the latter three cases. The normal resultant forces Fn (x-axis component Fnx and y-axis component Fny) perpendicular to fracture walls can affect direction of fluid flow. Fn can be obtained from the flow momentum equation and the Forchheimer law, see detailed explanations in Appendix A. The expressions of Fn for three different cases are given in Eq. (a.7) and Eq. (a.8), Eq. (a.15) and Eq. (a.16), and Eq. (a.23) and Eq. (a.24) respectively. From the expressions of Fn in Eq. (a.7) and Eq. (a.8), as the flow velocity increases, Fn for the case of BFIs increases and therefore a lateral flow would occur in the direction perpendicular to fracture walls. When Fn is large enough, this perpendicular flow will cause the formation of eddies, which in turn reduces the effective area for flow hence the effective aperture of the fracture. Consequently, nonlinear

5.5. Nonlinear coefficient B for BFI, CFI and FFI Significant research efforts have been made to examine the nonlinear coefficient B or Forchheimer coefficient β for fluid flow in a single fracture (e.g., [2,19]). However, the study of B is rarely reported for fluid flow in rock intersecting fracture.

11

Computers and Geotechnics 121 (2020) 103446

F. Xiong, et al.

Fig. 13. . The relationships of (a) Jc − α, (b) Jc − θ and (c) Jc − γ based on experiments and simulations (blue points: experimental data, black points: numerical results). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

According to the theory of flow loss [20] in angular pipes, the relationship between of angle φ, velocity v and head loss hLocal can be written as:

hLocal = ( 1 sin2

+

2

4 2 sin

2

)

v2 2g

where w is fracture width; m1, m2, m3, m4, m5 and m6 are regression coefficients; e is the fracture hydraulic aperture as defined above. Also plotted in Fig. 14 are the best-fitted curves of experimental data for BFI, CFI and FFI. The results show that Eq. (15), Eq. (16) and Eq. (17) fit well the experimental data with a high coefficient of determination (R2 > 0.80). m1, m2, m3, m4, m5 and m6 in this case are 0.103, 4.138, 3.93, 7.37e−7, 38.4 and 7.38e−7, respectively. The inevitable experimental error leads to the point of θ = 60° in CFI deviated from the predicted curve, shown in Fig. 14(b). Note these equations include the effects of fracture intersecting angle and fracture aperture. The effect of fracture surface roughness is not considered as discussed above.

(14)

where φ is exterior angle of the angular pipe and λ1 = 0.945 and λ2 = 2.047 for angular pipe. Based on those experimental data and extended numerical simulation, those relationships of B verse pipe intersecting angle are plotted in Fig. 14. Motivated by the characteristics of the curves plotted in Fig. 14 and knowing Q = evw, the semi-empirical relationship between B and α, θ or γ is proposed below using the similar form of Eq. (14):

B=

B=

B=

m1 cos2

( ) + m cos ( ) 4

2

2

2

m3

( )+m 2

4

sin4

() 2

(

1

2

2

) ) + m sin (abs ( 6

Le 2w 2

A sandstone block also from the site of Dagangshan Hydropower Project, China was used to produce a fractured rock specimen, which is a cylindrical specimen with the length and diameter of 100 mm and 50 mm. An increasing axial load was applied to the specimen until it just reaches the peak load when the load will be relived. Reasonably intact specimens with visible developed fractures were selected to do CT scans. On the CT scan image, fracture pixels were separated from

(16)

Le 2w 2

m5sin2 (abs

6.1. The fracture network creation in a rock specimen

(15)

Le 2w 2 sin2

6. Nonlinear flows in a fractured rock specimen

4

1

2

2

))

(17)

12

Computers and Geotechnics 121 (2020) 103446

F. Xiong, et al.

Fig. 14. . Curve-fitting of nonlinear coefficient B as a function of intersecting angle α, θ or γ based on experimental data and extended simulation results: (a) BFI (R2 = 0.99), (b) CFI (R2 = 0.85), (c) FFI (R2 = 0.90) (blue points: experimental data, black points: numerical results). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

those of rock matrices and pores. A serial of consecutive scan images was then combined to obtain the three dimensional geometry of fractures, which seen in Fig. 15(a). This method was also used by Yu et al. [34]. Pores and isolated fractures were then removed and only the clusters of fractures connected to both the top and bottom surfaces of the specimen were retained in the final DFN model, see the red surfaces shown in Fig. 15(b). There are 18 fractures with 28 intersections, and the dip angles of these fractures range from 20° to 80° in final fracture geometry model. ANSYS FLUENT was used to model the flow in the fracture network model. A total of 495,600 nodes and 1,690,000 elements were generated.

inside the fracture network was plotted. Fig. 17 shows the streamlines of fluid flow through the fracture network at the cross-section of y = 0.40 when the hydraulic gradient J = 1e−5 and 1e−3. It can be seen that at J = 1e−5, there is no eddy within the fracture and all the flow stream lines are parallel to the fracture walls. When J increases to 1e−3, eddies emerge for BFIs, CFIs and FFIs. These eddies cause the reduction of the effective flow area within the fracture. As a simplified approach, the proposed B model for fracture intersections is used to calculate the B value for the fracture network model:

6.2. Analysis of the experimental and numerical results

where N is number of fracture intersections, Bi is value of B calculated by Eqs. (15)–(17). In this example, the B value calculated by Eq. (13) is 9.0e10 s2/m6, which differs by < 20.3% compared with the result obtained from the experiment. Noted that the effect of fracture surface roughness and contacts are not considered in Eqs. (15)–(17), which are the main reasons for the deviation. However, this comparison demonstrated that the proposed B model can be used to estimate at a reasonable accuracy the nonlinear term in the Forchheimer law for a

N

B=

Bi i=1

The simulation results showing the relationship between hydraulic gradient J and flow rate Q are plotted in Fig. 16 where the experimental results also are displayed. Clearly J in this case exhibits a quadratic relationship with Q, which fits well the Forchheimer law. The values of A and B of the fitted curve are 22453.28 s/m3 and 1.13e11 s2/m6. To help understand better the nonlinear flow behaviour, the flow pattern

13

(18)

Computers and Geotechnics 121 (2020) 103446

F. Xiong, et al.

Fig. 15. . Three dimensional reconstruction of the fracture network model using CT: (a) volume modeling technique (b) three dimensional fracture network model used in flow simulations.

fracture network, even though potential interactions between different fracture intersections are not considered in Eqs. (15)–(17).

Corresponding extended numerical simulations were also conducted to supplement the analyses. A fractured rock specimen containing different types of fracture intersections was scanned by CT to construct its fracture network model. Both experimental tests and numerical simulations were performed on the specimen to investigate the nonlinear flow behaviour at the network scale. The experimental and numerical results show that no nonlinear flow was observed for SFI for the hydraulic gradients examined. For other

7. Conclusions In this work, flow tests were carried out on four kinds of fracture intersections: SFI, BFI, CFI and FFI to study the influences of fracture intersections on the nonlinear behaviour of the fluid flow.

14

Computers and Geotechnics 121 (2020) 103446

F. Xiong, et al.

Fig. 16. . The relationship between hydraulic gradient J and flow rate Q for the rock specimen.

three cases of fracture intersections: BFI, CFI and FFI, the flow nonlinearity is strongly affected by the fracture intersecting angle. The Forchheimer law still provides a good relationship to quantify the nonlinear flow regime. The nonlinearity of fluid flow was mainly caused by the formation of eddies at fracture intersections, which are clearly observable in numerical models. The occurrence of eddies is caused by the increasing normal resultant forces Fn perpendicular to fracture walls which changes the flow direction. The expressions of Fn for different types of fracture intersections are derived in this work based on the momentum balance. Based on experimental and extended numerical results, as the buckling angle increases, the Forchheimer coefficient β decreases for the case of BFI and the critical hydraulic gradient Jc increases. In contrary, the β for the case of CFI increases and Jc decreases with increasing crossing angle. The flow nonlinearity for the case of FFI depends on the values of furcating angle γ1 and γ2. The greater the abs(γ1γ2), the greater the value of β and the smaller the value of Jc. Based on laboratory observations and extended numerical results, three semi-empirical relationships (Eqs. (15)–(17)) that relate the nonlinear coefficient B in Forchheimer law to intersecting angle α, θ and γ were proposed for the fluid flow through BFI, CFI and FFI, respectively. At the fracture network scale, high hydraulic gradient causes the formation of more eddies at the fracture intersection of the cases of BFI, CFI and FFI. Consequently, the effective flow paths are narrowed, resulting in the nonlinear flow behaviour. The proposed B model for different cases of single fracture intersections can be used to predict the B value at the fracture network scale. This work focuses on the influences of three types of fracture intersections and the intersecting angles on the nonlinearity of the flow. It should be noted that other factors such as fracture surface roughness and fracture surface contacts can also strongly affect the flow nonlinearity in DFNs. More extensive experimental and numerical studies are required to study flow behavior for fracture intersections considering these factors.

Fig. 17. . The flow velocity profile and streamlines on y = 0.40 cross-section at different hydraulic gradients: (a) J = 1.0e − 5, (b) J = 1.0e − 3.

Qinghui Jiang: Conceptualization, Methodology, Software. Declaration of Competing Interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgement This research was financially supported by the National Science Foundation of China (No. 51679173) and the Natural Science Foundation of Hubei Province (No. 2016CFA083). These supports are gratefully acknowledged. The first author also acknowledges the support from the Chinese Scholarship Council for awarding a scholarship (CSC No. 201806270149) to carry out the study presented in this paper as a visiting research student at the University of Adelaide.

CRediT authorship contribution statement Feng Xiong: Data curation, Writing - original draft. Wei Wei: Software, Validation. Chaoshui Xu: Writing - review & editing.

15

Computers and Geotechnics 121 (2020) 103446

F. Xiong, et al.

Appendix A. Normal force perpendicular to fracture walls based on the Forchheimer law The positive directions of forces and momentums are defined as the positive direction of the X-axis and Y-axis. (A) BFI A representative unit thickness element of a BFI is shown in Fig. a1. The following equations can be obtained based on the conservation of momentums in X and Y directions:

Fx =

b1 v12 cos

1 2

b2 v22 cos

1 2

(a.1)

Fy =

b1 v12 sin

1 2

b2 v22 sin

1 2

(a.2)

where v1 and v2 are mean velocity. F (Fx and Fy) is the external force applied to the element, mainly due to pressure difference, drag force Fτ along fracture walls and normal force Fn perpendicular to fracture walls. Considering the force balance of the element leads to the following equations: (a.3)

Fx = F x + Fnx Fy =

P1 b1 sin

1 2

+ P2 b2 sin

1 2

+ F y + Fny

(a.4)

The drag force Fτ is the viscous resistance of fracture surfaces to fluid flow, which acts in the direction opposite to flow [6]. Hence, Fτ is:

F x = gb1 L1 J1 cos

1 2

+ gb2 L 2 J2 cos

1 2

(a.5)

F y = gb1 L1 J1 sin

1 2

+ gb2 L 2 J2 sin

1 2

(a.6)

where J1 and J2 are hydraulic gradient for segment 1 and segment 2, shown in Fig. a1, respectively. Considering L1 = L2 = L, P2 = 0 as well as the Forchheimer law, Fn can be obtained by combining Eq. (a.1)–(a.6) and setting the flow rate Q = vb:

Fnx = Fny = (

[ (b1 v12 + b2 v22) + gL ((Av1 b1 + Bv12 b12 ) b1 + (Av2 b2 + Bv22 b22 ) b2)] cos (b1 v12 + b2 v22) + P1 b1

1 2

gL ((Av1 b1 + Bv12 b12 ) b1 + (Av2 b2 + Bv22 b22 ) b2)) sin

(a.7)

1 2

(b) CFI

Fig. a1. . Schematic diagram of flow for a BFI. 16

(a.8)

Computers and Geotechnics 121 (2020) 103446

F. Xiong, et al.

A representative unit thickness element of a CFI is shown in Fig. a2. The following equations can be obtained based on the conservation of momentums in X and Y directions:

Fx = b1 v12 sin

1 2

b2 v22 sin

1 2

b3 v32 sin

1 2

+ b4 v42 sin

1 2

(a.9)

Fy = b1 v12 sin

1 2

b2 v22 sin

1 2

b3 v32 sin

1 2

+ b4 v42 sin

1 2

(a.10)

Considering the force balance of the element leads to the following equations:

Fx = P1 b1 sin

Fy =

1 2

P1 b1 cos

P1 b2 sin

1 2

1 2

1 2

P0 b3 sin

1 2

P1 b2 cos

P0 b3 cos

+ P0 b4 sin

1 2

1 2

P0 b4 cos

+ F x + Fnx

1 2

(a.11)

+ F y + Fny

(a.12)

Similarly, the drag force Fτ is:

Fx =

gb1 L1 J1 sin

gb3 L3 J3 sin

( ) 1 2

F y = gb1 L1 J1 cos gb3 L3 J3 cos

( )+ 1 2

gb2 L2 J2 sin

gb4 L4 J4 sin

( )+ 1 2

( )+ 1 2

( )

1 2

1 2

gb2 L 2 J2 cos

gb4 L4 J4 cos

( )+

( )

(a.13)

( )+ 1 2

1 2

(a.14)

where J1, J2, J3 and J4 are hydraulic gradient for segment 1, segment 2, segment 3 and segment 4 shown in Fig. a2, respectively. Given L1 = L2 = L3 = L4 = L, P2 = 0 as well as the Forchheimer law, Fn can be obtained by combining Eqs. (a.9)–(a.14):

Fnx = (b1 v12

b2 v22

(b1 (Av1 b1 + Bv12 b12) Fny =

b3 v32 + b4 v42) sin

( ) + (P b 1 2

1 2

b2 (Av2 b2 + Bv22 b22)

(b1 v12 + b2 v22 + b3 v32 + b4 v42) cos

P1 b1) sin

( )+ 1 2

gL sin

( ) 1 2

b3 (Av3 b3 + Bv32 b32 ) + b4 (Av3 b3 + Bv32 b32 ))

( ) + (P b 1 2

1 1

+ P1 b2) cos

( ) 1 2

gL cos

(a.15)

( ) 1 2

(b1 (Av1 b1 + Bv12 b12) + b2 (Av2 b2 + Bv22 b22) + b3 (Av3 b3 + Bv32 b32 ) + b4 (Av3 b3 + Bv32 b32 ))

(a.16)

(c) FFI A representative unit thickness element of a FFI is shown in Fig. a3. The following equations can be obtained based on the conservation of momentums in X and Y directions:

Fig. a2. . Schematic diagram of flow for a CFI. 17

Computers and Geotechnics 121 (2020) 103446

F. Xiong, et al.

Fig. a3. . Schematic diagram of flow for a FFI.

Fx = b1 v12 + b2 v22 cos( 2) + b2 v22 cos( 1)

(a.17)

b2 v22 sin( 2 ) + b3 v32 sin( 1)

(a.18)

Fy =

where γ1 is angle between segment 1 and segment 3, and γ2 is angle between segment 1 and segment 2. Considering the force balance of the element leads to the following equations:

Fx = P1 b1

P2 b2 cos( 2)

Fy = P2 b2 sin( 2 )

(a.19)

P3 b3 cos( 1) + F x + Fnx

(a.20)

P3 b3 sin( 1) + F y + Fny

The drag force Fτ is:

Fx =

gb1 L1 J1

gb2 L2 J2 cos( 2 )

Fy =

gb2 L2 J2 sin( 2 )

(a.21)

gb3 L3 J3 cos( 1)

(a.22)

gb3 L3 J3 sin( 1)

where J1, J2, J3 are hydraulic gradient for segment 1, segment 2, segment 3 shown in Fig. a3, respectively. Given L1 = L2 = L3 = L, P2 = 0 as well as the Forchheimer law, Fn can be obtained by combining Eqs. (a.17)–(a.22):

Fnx = (b1 v12 + b2 v22 cos( 2) + b3 v32 cos( 1))

P1 b1+

+ gL (b1 (Av1 b1 + Bv12 b12) + b2 (Av2 b2 + Bv22 b22 ) cos( 2) + b3 (Av3 b3 + Bv32 b32 ) cos( 1))

Fny = (

+ gL ( b2 (Av2 b2 +

+ b3 v32 sin( 1)) Bv22 b22) sin( 2 ) +

(a.23)

b2 v22 sin( 2 )

b3 (Av3 b3 + Bv32 b32) sin( 1))

(a.24)

References [15]

[1] Bandis SC, Lumsden AC, Barton NR. Fundamentals of rock joint deformation. Int J Rock Mech Min Sci Geomech Abstr 1983;20(6):249–68. [2] Chen Y, Zhou J, Hu S, Hu R, Zhou C. Evaluation of Forchheimer equation coefficients for non-Darcy flow in deformation rough-walled fractures. J Hydrol 2015;529:993–1006. [3] Cherubini C, Giasi C, Pastore N. Bench scale laboratory tests to analyze non-linear flow in fractured media. Hydrol Earth Syst Sci 2012:2511–22. [4] Ghane E, Fausey N, Brown L. Non-Darcy flow of water through woodchip media. J Hydrol 2014;519:3400–9. [5] Holditch SA, Morse RA. The effects of non-Darcy flow on the behavior of hydraulically fractured gas wells. J Petrol Technol 1976;28(10):1169–79. [6] Huilgol R, Phan-Thien N. Fluid Mechanics of Viscoelasticity. New York: Elsevier; 1997. [7] Hyman J, Karra S, Makedonska N, W.Gable C, Painter S, Viswanathan H. DFNWORKS: A discrete fracture network for modeling subsurface flow and transport. Comput Geosci 2015; 84: 10-19. [8] Ishibashi T, Watanabe N, Hirano N, Okamoto A, Tsuchiya N. GeoFlow: A novel model simulator for prediction of the 3-D channeling flow in a rock fracture network. Water Resour Res 2012;48:W07601. [9] Javadi M, Sharifzadeh M, Shahriar K, Mitani Y. Critical Reynolds number for nonlinear flow through rough-walled fractures: The role of shear processes. Water Resour Res 2014;50(2):1789–804. [10] Jiang Q, Yao C, Ye Z, Zhou C. Seepage flow with free surface in fracture networks. Water Resour Res 2013;49(1):176–86. [11] Karra S, Makedonska N, Viswanathan H, Painter S, Hyman DJ. Effect of advective flow infractures and matrix diffusion on natural gas production. Water Resour Res 2015;51(10):8646–57. [12] Li B, Liu R, Jiang Y. Influences of hydraulic gradient, surface roughness, intersecting angle and scale effect on nonlinear flow behaviour at single fracture. J Hydrol 2016;538:440–53. [13] Li K, Ma M, Wang X. Experimental study of water flow behaviour in narrow fractures of cementitious materials. Cem Concr Compos 2011;33(10):1009–13. [14] Liu R, Jiang Y, Li B. Effects of intersection and dead-end of fractures on nonlinear

[16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28]

18

flow and particle transport in rock fracture networks. Geosci J 2016(a); 20(3): 415–426. Liu R, Li B, Jiang Y. Critical hydraulic gradient for nonlinear flow through rock fracture networks: The roles aperture, surface roughness, and number of intersections. Adv Water Resour 2016(b); 88: 53–65. Qin Z, Fu H, Chen X. A study on altered granite meso-damage mechanisms due to water invasion-water loss cycles. Environ Earth Sci 2019;78:428. Radilla G, Nowamooz A, Fourar M. Modeling non-Darcian single-and two-phase flow in transparent replicas of rough-walled rock fractures. Transp Porous Media 2013;98(2):401–26. Ren S, Gragg S, Zhang YJ, Carr B, Yao G. Borehole characterization of hydraulic properties and groundwater flow in a crystalline fractured aquifer of a headwater mountain wastershed, Laramine Range, Wyoming. J Hydrol 2018; 561: 780–795. Rong G, Yang J, Cheng L, Zhou C. Laboratory investigation of nonlinear flow characteristics in rough fractures during shear process. J Hydrol 2016;541:1385–94. Rennels D, Hudson H. Pipe flow: A practical and comprehensive guide. John Wiley & Sons; 2012. Sævik P, Nixon C. Inclusion of topological measurements into analytic estimates of effective permeability in fractured media. Water Resour Res 2017;53(11):9424–43. Sun W, Du H, Zhou F, Shao J. Experimental study of crack propagation of rock-like specimens containing conjugate fractures. Geomech Eng 2019;17(4):323–31. Tzelepis V, Moutsopoulos NK, Papaspyros J, Tsihrintzis V. Experimental investigation of flow behavior in smooth and rough artificial fractures. J Hydrol 2015; 521: 108–118. Wang M, Chen Y, Ma G, Zhou J, Zhou C. Influence of surface roughness on nonlinear flow behaviours in 3D self-affine rough fractures: Lattice Boltzmann simulations. Adv Water Resour 2016;96:373–88. Wang Z, Xu C, Dowd P. A modified cubic law for single-phase saturated laminar flow in rough rock fractures. Int J Rock Mech Min Sci 2018;103:107–15. Wang Z, Xu C, Dowd P. Perturbation Solutions for Flow in a Slowly Varying Fracture and the Estimation of Its Transmissivity. Transp Porous Media 2019:1–25. Wang J, Li S, Li L, Xu Z, Gao C. Attribute recognition model for risk assessment of water inrush. Bull Eng Geol Environ 2019;78:1057–71. Xiong F, Jiang Q, Xu C, Zhang X. Influences of connectivity and conductivity on nonlinear flow behaviours through three-dimension discrete fracture networks.

Computers and Geotechnics 121 (2020) 103446

F. Xiong, et al. Comput Geotech 2019;107:128–41. [29] Xiong F, Jiang Q, Ye Z, Zhang X. Nonlinear flow behaviour through rough-walled rock fractures: The effect of contact area. Comput Geotech 2018;102:179–95. [30] Xu C, Fidelibus C, Dowd P, Wang Z, Tian Z. An iterative procedure for the simulation of the steady-state fluid flow in rock fracture networks. Eng Geol 2018;242:160–8. [31] Yao C, Jiang Q, Shao J. A numerical analysis of permeability evolution in rocks with multiple fractures. Transp Porous Media 2015;108(2):289–311. [32] Ye Z, Jiang Q, Yao C, Liu Y, Cheng A, Huang S, et al. The Parabolic Variational Inequalities for Variably Saturated Water Flow in Heterogeneous Fracture Networks. Geofluids 2018. [33] Ye Z, Jiang Q, Zhou C, Liu Y. Numerical analysis of unsaturated seepage flow in two-dimensional fracture networks. Int J Geomech 2016;17(5):04016118. [34] Yu Q, Yang S, Ranjith P, Zhu W, Yang T. Numerical modeling of jointed rock under compressive loading using X-ray computerized tomography. Rock Mech Rock Eng 2015;49(3):877–91.

[35] Zeng Z, Grigg R. A criterion for non-Darcy flow in porous media. Transp Porous Media 2006;63(1):57–69. [36] Zhou J, Hu S, Chen Y, Wang M, Zhou C. The friction factor in the Forchheimer equation for rock fractures. Rock Mech Rock Eng 2016;49(8):3055–68. [37] Zhou J, Hu S, Fang S, Chen Y, Zhou C. Nonlinear flow behavior at low Reynolds numbers through rough-walled fractures subjected to normal compressive loading. Int J Rock Mech Min Sci 2015;80:202–18. [38] Zimmerman R, Al-Yaarubi A, Pain C, Grattoni C. Non-linear regimes of fluid flow in rock fractures. Int J Rock Mech Min Sci 2004;41:163–9. [39] Zimmerman R, Bodvarsson G. Hydraulic conductivity of rock fractures. Transp Porous Media 1996;23(1):1–30. [40] Zou L, Jing L, Cvetkovic V. Roughness decomposition and nonlinear fluid flow in a single rock fracture. Int J Rock Mech Min Sci 2015;75:102–18. [41] Zhang Z. Jan Nemcik. Fluid flow regimes and nonlinear flow characteristics in deformable rock fractures. J Hydrol 2013;477:139–51. https://doi.org/10.1016/j. jhydrol.2012.11.024 In this issue.

19