Effects of angle patterns at fracture intersections on fluid flow nonlinearity and outlet flow rate distribution at high Reynolds numbers

International Journal of Rock Mechanics & Mining Sciences 124 (2019) 104136

Contents lists available at ScienceDirect

International Journal of Rock Mechanics and Mining Sciences journal homepage: http://www.elsevier.com/locate/ijrmms

Effects of angle patterns at fracture intersections on fluid flow nonlinearity and outlet flow rate distribution at high Reynolds numbers L.F. Fan a, H.D. Wang a, Z.J. Wu b, *, S.H. Zhao c a

College of Architecture and Civil Engineering, Beijing University of Technology, Beijing, 100124, China School of Civil Engineering, Wuhan University, Wuhan, 430072, China c School of Geosciences and Info-physics, Central South University, Changsha, 410083, China b

A R T I C L E I N F O

A B S T R A C T

Keywords: Nonlinear flow Fracture intersection Reynolds numbers Angle pattern Flow rate distribution

Fluid flow tests on a series of one-inlet-two-outlet specimens with fracture intersections are conducted under conditions of high hydraulic gradient (J ¼ 100–101) and high Reynolds number (Re ¼ 104–105). The fluid flow through the intersections is recorded using a single lens reflex (SLR) camera. It was demonstrated that the intersection shapes can enhance the strong nonlinearity of the fluid flow at high J and high Re. The angle pattern, which is represented by the case of intersecting fractures, affects the ratio of the outlet flow rate distribution (ε). The detailed effects for all cases can be categorized into three models via quantitative analysis. The intersection shapes of model I angle patterns result in larger ε. However, the intersection shapes of model II and model III angle patterns result in smaller ε compared. In addition, the distribution fitting curve of the peak ε values is provided. A modified Gaussian distribution function is formulated to fit the experimental curves. Furthermore, the relationship between ε and the angle patterns is elaborated.

1. Introduction

extensively studied by many researchers.23–34 Scesi and Gattinoni35 developed a relationship that enables estimation of the hydraulic con­ ductivity tensor (an essential parameter for understanding fluid flow in fractured rock masses) by taking the joint roughness into consideration. Ju et al.36 proposed a fractal model to relate the fluid resistance to the fracture roughness, based on which a fractal equivalent permeability coefficient for a single rough fracture was formulated. Meanwhile, the rigorous disturbances in the hydraulic behaviour through a single frac­ ture during shear force also became a research focus.37–47 Apart from fluid flow in a single fracture, fluids flow through various fracture intersections are of equal importance to understand, simulate and analysis fluids flow in fracture networks. It also attracts many re­ searchers’ attentions. Wilson and Witherspoon48 conducted a series of laboratory experiments using circular conduits to determine the magnitude of laminar flow interference effects at fracture intersections. The results indicated that the interference effects at fracture in­ tersections are negligibly small in most fracture systems when flow is in the laminar regime. However, if the flow interference effects are of the same order of magnitude for platelike fractures as for circular conduits, they concluded that the head loss at a fracture intersection would be equivalent to a length of approximately five mechanical apertures (e) for flow rates with Reynolds number (Re) of approximately 100. Tian49

Fractures often are the main flow channels in naturally fractured rock masses.1 Therefore, fluid movement in fracture networks (see Fig. 1) is of critical importance to engineering applications in fractured rock masses. Due to the strongly varying velocity fields in different in­ dustrial and engineering fields, such as geothermal energy extraction,2–6 water resources management,7–11 underground tunneling,12,13 geothermal fluid or hydrocarbon exploitation14,15 and hydraulic splitting,5,10,16,17 the fluid flow in fractured rock masses is highly com­ plex. Reynolds numbers remain a large magnitude for the actual con­ ditions of high-pressure hydraulic splitting or fluid flow under high pressure and high stress in deep rock masses. In these conditions, the angles and geometries at the fracture intersections have a significant effect on the flow state among the fractures.18–20 In addition, the flow rate at the fracture intersections is very important in determining fluid flow through a fracture network.21 Therefore, it is meaningful to sys­ tematically investigate the effects of fracture intersections on the fluid flow characteristics including fluid flow nonlinearity and outlet flow rate distributions. Since Lomize22 first established the cubic law using the parallel plate model, fluid flow through a single fracture without shear has been * Corresponding author. E-mail address: [email protected] (Z.J. Wu).

https://doi.org/10.1016/j.ijrmms.2019.104136 Received 7 October 2018; Received in revised form 12 April 2019; Accepted 23 October 2019 Available online 1 November 2019 1365-1609/© 2019 Elsevier Ltd. All rights reserved.

L.F. Fan et al.

International Journal of Rock Mechanics and Mining Sciences 124 (2019) 104136

fractures as inlets and two fractures as outlets.21,48,51,53,54 As shown in Fig. 2, the form of one fracture as inlet and two fractures as outlets, among which the fluid flow undergoes redistribution at the fracture intersection, is also a fundamental element in the discrete fracture net­ works (DFNs). Li et al..19 used laboratory experiments and numerical simulations to study one-inlet-two-outlet models. The results showed that the nonlinear degree of the hydraulic gradient (J) and Q increases with increasing J and increasing joint roughness coefficient (JRC) while with decreasing fracture length. In addition, they found that the inter­ secting angle affects the effluent fracture flow rate (q) by three inter­ secting angle patterns. Liu et al.56 used numerical simulations of DFNs to conclude that when J drops below a critical hydraulic gradient (Jc), the relationship between J and Q reduces to the widely used cubic law. They also demonstrated that larger apertures, rougher fracture surfaces and a larger number of intersections in a DFN may result in the onset of nonlinear flow at a lower Jc. However, the effects of the intersection shapes on the hydraulic properties of fluid flow at fracture intersections have not been studied systematically. In addition, few efforts have focused on fluid flow through one-inlet-two-outlet intersecting fractures at high Re. The main purpose of this paper is to experimentally investigate the effects of angle patterns at fracture intersections on fluid flow nonline­ arity and outlet flow rate distribution at high Re. First, high-precision fluid flow tests were conducted out on a series of one-inlet-two-outlet specimens with 51 different angle patterns. The syringe pressure gradient was varied in the range of 0.20–1.80 MPa to ensure that the flow nonlinear regime among the intersecting fractures was at a large magnitude Re. For each angle pattern, flow test was performed with different J. Second, three models were used to evaluate the effects of angle patterns on the outlet flow rate distribution. Finally, the redistri­ bution law for all angle patterns in a 2D plane was provided using different q during fluid flow laboratory investigations.

found the deflection flow phenomena at fracture intersections for the case when the fluid at the narrow fracture fully flows into the wide fracture, leading to increased flow at the wide fracture and decreased flow at the narrow fracture which belongs to the water storage struc­ tures. Hull and Koslow21 recognized that the assumption of complete mixing at junctions may be unrealistic and conducted a series of labo­ ratory tests to determine the routing criteria for streamlines through fracture junctions. Based on the test results, they concluded that streamlines cannot cross continuous junctions, while flow along adja­ cent streamlines must be in the same direction for discontinuous junc­ tions. Philip50 applied fluid mechanics to investigate the streamline routing treatment for solute transport in fracture networks and found that the errors of proportional routing may be unacceptably large when the two inlet discharges and the two outlet discharges both differ greatly in magnitude. By using numerical simulations to examine the variability of flow patterns in realistic fracture intersection geometries rather than perfectly orthogonal intersections of two-inlet-two-outlet models, Kosakowski and Berkowitz18 found that the intersection geometries and the boundary conditions can lead to significant nonlinear inertial effects for Re of 1–100 due to a strong interplay among the flow velocity (v). Johnson and Brown51,52 used experimental and numerical methods based on two-inlet-two-outlet models to show that flow channelization through rough-walled intersecting fractures significantly enhances physical mixing compared to intersecting parallel plates, while surface roughness may reduce solute dilution and increase solute dispersion. Based on theoretical formulae and a series of experiments of two intersecting fractures, Hu et al.53 found that local head losses exist at fracture intersections and that the intersection angle affects the flow distribution at fracture intersections. Based on the results of laboratory experiments and mathematical analyses, Ji et al.54 concluded that fracture intersections with differing apertures are a more significant factor in controlling the network-scale phase structure for non-wetting flow than for wetting flow. Zafarani and Detwiler55 conducted a novel probabilistic method that enables accurate and efficient calculation of the Peclet number-dependent flow and transport through fracture in­ tersections. Their model accurately represents the trajectory and travel times of particles passing through a fracture intersection and is also 2–3 orders of magnitude faster than traditional simulations that solve the Stokes equations. The above discussions show that most of the previous studies mainly focused on fluid flow through a fracture intersection in the form of two

2. Methodology 2.1. Theoretical background The general description of fluid flow in a single fracture is provided by the Navier-Stokes (NS) equation, which expresses momentum and mass conservation over the fracture void space.12 For a steady incom­ pressible Newtonian fluid, the NS equation, which is derived from

Fig. 1. Natural fractures and fracture intersections in rock masses. 2

L.F. Fan et al.

International Journal of Rock Mechanics and Mining Sciences 124 (2019) 104136

Fig. 2. A 2D DFN composed of a large number of fracture intersections. Schematic views of fluid flow through fracture intersections are also given for the cases of a one inlet and one outlet, b one inlet and two outlets, c two inlets and one outlet, and d two inlets and two outlets (excerpted from Li et al., 2016).

Newton’s second law,57–59 is written as follows:

ρ½∂u = ∂t þ ðu ⋅ rÞu� ¼

rP þ r ⋅ T þ ρf

where a is linearly proportional to the reciprocal of the fracture permeability.19,56,72 Note that 1D simplification of NS equation is adopted in this study, which is based on the following reasons: (1) the height and width of a fracture is much smaller than its length; (2) the cross section of the fracture is not variable along the length of the fracture; (3) neglecting the variation of fluid velocity in the cross section. For fluid flow among fractures, the ratio of the inertial force to the viscous force is defined as the Reynolds number (Re), which is expressed as32,66:

(1)

where u is the flow velocity, ρ is the fluid density, P is the pressure, T is the shear stress tensor, t is the time, and f is the body force. Considering that the convective acceleration term, ðu ⋅rÞu, is the only source of nonlinearity that is strongly affected by the geometric variation of void spaces due to complex surface roughness of fractures,60 the NS equation can be reduced to the linear Stokes equation by assuming that the in­ ertial forces in the fluid flow are negligible compared with the viscous force and pressure force. Furthermore, for laminar flow in single frac­ tures under a one-dimensional pressure gradient between two smooth parallel plates, which is known as Poiseuille flow, the linear Stokes equation is reduced to: Q¼

wρge3H J; 12μd

Re ¼

(2)

v¼

J¼

h2 1 P1 P2 1 ΔP 1 ¼ rP; ¼ ¼ ΔL ρg ΔL ρg ΔL ρg

(3)

q q2

ε ¼ 1;

(4)

(7)

(8)

where ε is the ratio of the outlet flow rate distribution, q1 is the larger outlet flow rate, and q2 is the smaller outlet flow rate. In this paper, due to the need to know the dispersion degree of a set of data from the tests, the coefficient of variation (also known as the discrete coefficient in statistics) is used to compare the discrete degree of different levels of variables:

where h1 and h2 are the hydraulic heads at each end of a fracture, P1 and P2 are the pressures at each end of a fracture, and ΔL is the effective path under the hydraulic gradient. Eq. (3) can be written in the following Forchheimer’s law-based form as: J ¼ aQ þ bQ2 ;

Q ; S2

where v is the flow velocity and S is the cross-sectional area of the fracture. S can be calculated as S ¼ ew, where e and w are accurately measured using digital vernier caliper. For there are two outlets in the present study, the ratio of q at one outlet to the total Q may change depending on the varying angle pat­ terns. Based on the ratio methods used to study the outlet flow rates for cases of two inlets and two outlets,21,51,52 the ratio of outlet flow rate distribution is defined as:

where rP is the pressure gradient, a’ and b’ are coefficients that are related to the fracture geometry, such as aperture and roughness. Eq. (3) is the famous Forchheimer’s law,67 which adequately describes the nonlinear flow through rock fractures, especially for strong inertia regimes.68–71 J is proportional to rP: h1

(6)

Equation (6) indicates that Re is linearly correlated with Q, while J is quadratically correlated with Re, by substituting Eq. (6) into Eq. (5). J is linearly correlated with Re at small Q. Re is typically utilized to detect the onset of nonlinear flow in single fractures or single-fracture in­ tersections that only involve single or several fracture segments.32,66 Given that the test in this paper is performed in a single intersecting fracture, u (the vector flow velocity) can be reduced to one-dimensional flow velocity:

where Q is the total volumetric flow rate, g is gravitational acceleration, J is the hydraulic gradient, w is the width of the fracture, eH is the hy­ draulic aperture, and μd is the dynamic viscosity. This equation is commonly called the “cubic law”61 and predicts that Q is proportional to J. Eq. (2) also implies that Q is linearly correlated with J, which is valid when Re is sufficiently low. With the increase of Q and the increase of Re, more than the proportional increases in Q that is known as nonlinear flow62–66: rP ¼ a’Q þ b’Q2 ;

ρQ : μd w

(5)

σ μ

CV ¼ � 100%;

where a ¼ a’=ρg and b ¼ b’=ρg. Forchheimer’s law can be applied over the entire range of Q, where it effectively reduces to Darcy’s law at small Q by dropping the nonlinear term, bQ.2 Eq. (5) reduces to J ¼ aQ,

(9)

where CV is the coefficient of variation, σ is the standard deviation of a set of data, and μ is the average of a set of data. 3

L.F. Fan et al.

International Journal of Rock Mechanics and Mining Sciences 124 (2019) 104136

The Gaussian distribution (also called the normal distribution) is an extremely important probability distribution in mathematics, physics, engineering and other fields. Due to the shape of the curve, with low ends, a high middle and axial symmetry, the Gauss distribution is commonly known as the bell curve. The mathematical expectation de­ termines the bell curve’s position, namely, the translation on the x-axis, and the standard deviation determines the bell curve’s amplitude, namely, the single peak width. Because the weighting of the variables in this paper is 1, the mathematical expectation is equal to the average. There is only one dimension for the variables, and the probability den­ sity function of the Gauss distribution based on one-dimensional is given as follows: 1 ðx μÞ2 f ðxÞ ¼ pffiffiffiffiffiffiffiffie 2σ2 ; 2πσ

(10)

where σ determines the horizontal shift of the function, and μ de­ termines the peak width. 2.2. Flow tests

Fig. 4. Final PMMA specimen with a fracture intersection.

2.2.1. Specimen preparation PMMA (polymethyl methacrylate) material was used to machine the specimens of fracture intersections in this present study. In order to validate the PMMA specimens, a low-permeability granite was also used to machine a specimen. It was found that the ε - J curves of the rock and PMMA specimens matched well in all three groups of contrast experi­ ments for l ¼ 283 mm (see Fig. 3). PMMA was adopted to replace the rock to machine specimens for the following reasons: (1) the PMMA material can be machined into accurate specimens using computer nu­ merical control machine tools; (2) the fractures can be sealed well using PMMA adhesive; (3) the rock matrix parameters were not considered in this study. After observing that the intersection shapes have a significant effect on the outlet flow rate distribution when J > 10 2 and Rr (radius) < 50–200 mm19, the tests were conducted with e ¼ 2 mm, l ¼ 4 mm and w ¼ 5.6 mm. The final PMMA specimens are shown in Fig. 4, with e ¼ 2.108–2.126 mm, CVe ¼ 0.432%; l ¼ 3.995–4.007 mm, CVl ¼ 0.153%; and w ¼ 5.648–5.745 mm, CVw ¼ 0.851%. The values of CV are obtained based on Eq. (9), which are all small, and indicates that specimens were machined accurately.

test. Tap water was considered as a sufficient water supply. To avoid contamination from impurities in the tap water, which can deposit among fractures and change the fracture apertures, a water purifier was attached at the start of the experimental system. A high-pressure syringe pump with a designed pressure of 2.0 MPa, effective volume of 40 L, and syringe pressure accuracy of �0.01 MPa was accurately controlled using a digital operation panel. The varying range of the syringe pressure gradient was 0.20–1.80 MPa, which resulted in a syringe flow rate of 2.00–9.00 L/min (3.33 � 10 5–1.50 � 10 4 m3/s). Thin PMMA plates as the same size as the specimens were glued to the fractures specimens to ensure proper sealing of the fractures. The specimens were fixed on a horizontal platform using transparent tape to guarantee that the initial pressure difference between the inlet and outlets was negligibly small. Moreover, the horizontal platform height could be adjusted as needed. Two buckets, each with a side hole, were used as the collectors. The size of the platform was slightly smaller than the specimens so that the outlets could be inserted into the side hole of the buckets to ensure that the fluid flow at the outlets could completely freely jet into the buckets. Experiments indicated that the tubule, which can produce confining pressure, was equivalent to extension of the fracture. Li et al.19 concluded that an excessively long fracture can produce outlet flow rates close to the average distribution, i.e., ε close to a fixed value of 1. In addition, cases with ε close to 1 could result in major challenges when analysing and summarizing the experimental results. Hence, the drainage form was conducted by free jet rather than tubule diversion. The effluent at each outlet was collected by the buckets. Electronic balances with an accuracy of �0.02 g were used to measure the weight of the effluents. Each electronic balance was connected to a computer via an R232 data cable to record the weight of the effluents in real time. A digital pressure gauge with an accuracy of �0.1 kPa was used to measure the pressure between the inlet and outlets. An SLR (single lens reflex) camera, which was installed on the upper part of the experi­ mental specimen, was used to record the experimental phenomena via the transparent polymethyl methacrylate (PMMA) material.

2.2.2. Experimental setup Fig. 5 shows a schematic view of the experimental system of the flow

2.2.3. Test procedure Before the test, the air in the experimental system was removed using a vacuum pump. The entire flow path, including the connecting tubules, inflow butt joint and all fractures, were simultaneously filled with pu­ rified water. After the specimen was fixed by transparent tape, a fracture was chosen as the influent fracture optionally. The inlet of the influent fracture was quickly connected to the inflow butt joint via a transparent pressure-proof PVC tubule with a 6 mm inner diameter and 8 mm outer diameter. The connection between the PVC tubule and the inlet was sealed using a hot glue gun. The two outlets of the specimen were

Fig. 3. Relationship between ε and J for rock and PMMA specimens at l ¼ 283 mm. 4

L.F. Fan et al.

International Journal of Rock Mechanics and Mining Sciences 124 (2019) 104136

Fig. 5. Experimental system for the fluid flow test.

respectively aimed at the hole of the collector to ensure that the fluid flow could jet into the collector. The test procedure followed five steps: (1) open the inlet valve and close the outlet valve of the pump to inject purified water into the pump; (2) after the pump is filled with purified water, set 0.2 MPa as the initial supply pressure through the digital operation panel. Then, open the inlet valve of the influent fracture to supply stable fluid flow for 1 min. At the same time, record the measured data. Finally, close the valve; (3) change the supply pressure to obtain data at different J, referring to Eq. (4); (4) change the influent fracture and repeat the test; (5) use the SLR camera to record the experimental phenomena. According to Eq. (4), where P1 ¼ P_inlet, P2 ¼ P0, because P0 is atmospheric pressure, and ΔL ¼ l, the values of J in both effluent fractures are equal.

3.1. Nonlinearity of the fluid flow among intersecting fractures The specimen whose image is shown in Fig. 7(e) (or Fig. 7(g)) was used to conduct three groups of experiments for three fractures as the influent fracture in turn. The relationship between J and Q (or q) is shown in Fig. 8(a), when fracture_1 & 2, fracture_2 & 3 and fracture_1 & 3 are used as the effluent fractures (also called outlet_1 & 2, outlet_2 & 3 and outlet_1 & 3), in accordance with α ¼ 60� & β ¼ 120� , α ¼ 30� & β ¼ 60� and α ¼ 30� & β ¼ 180� . The trends of Q, measured in the ex­ periments and calculated using the cubic law, with increasing J are shown in Fig. 8(b). Referring to Eq. (4), Eq. (6) and Eq. (7), it can be shown that for α ¼ 60� & β ¼ 120� , with J varying from 1.184 to 30.494, Re increases from 5782.456 to 21981.287 and v increases from 2.690 to 10.224 m/s; for α ¼ 30� & β ¼ 60� , with J varying from 1.270 to 30.753, Re increases from 6107.018 to 23757.894 and v increases from 2.841 to 11.050 m/s; and for α ¼ 30� & β ¼ 180� , with J varying from 1.356 to 30.753, Re increases from 5819.298 to 23769.006 and v increases from 2.707 to 11.055 m/s. As shown in Fig. 8(a), for J (100 - 101) and Re (104 - 105), the rela­ tionship between J and Q (or q) is strongly nonlinear. J is quadratically correlated with Q (or q), which agrees with Forchheimer’s law. Li et al.19 also observed a quadratic relationship between J and Q at J ¼ 10 1 - 100 using the same angle patterns. As shown in Fig. 8(b), for higher J, the deviations in the cubic law are enhanced with increasing J. The inertial force has a larger effect on the nonlinear flow regime than the viscous force for J > 10 1 and Rr < 50–200 mm19. In addition, the pressure difference occurs in the direction perpendicular to v due to the sharp change in e at fracture intersections, which causes the occur­ rence of eddies.73 As a result, for Re ¼ 104 - 105, nonlinear flow becomes turbulent with lots of small eddies,70 which results in strong nonlinearity of the fluid flow through the fracture intersections.

3. Results and discussion As shown in Fig. 6, to analyse the experimental results systemati­ cally, β (0� < β < 180� ) is defined as the intersection angle between both effluent fractures to represent the radiation range of the effluent frac­ tures, while α (0� � α � 180� ) is defined as the included angle between the extended orientation and the angle bisector of β to represent the orientation of the effluent fractures relative to the influent fracture. It is emphasized that α and β are independent of each other, that is, the combinations of α and β represent cases where two random fractures distribute in the 2D plane. The effluent fracture with a larger flow rate is called the advantageous fracture, while the effluent fracture with the smaller flow rate is called the disadvantageous fracture. Fig. 7 exhaustively shows the distribution of the specimen angle patterns in the domain of the α - β plane. In the α - β plane, the blank blocks represent non-existent angle patterns based on the definitions in this paper.

3.2. Mechanism analysis of the effects of angle patterns on the outlet flow rate distribution Liu et al.56,72 concluded that intersection shapes have a serious effect on the fluid flow nonlinear regimes. In fact, the intersection shapes greatly influence ε. The real photos of the specimens in Fig. 7 visually present the various effects when the fluid flow rushes through the fracture intersections under various angle patterns. Based on all these cases, the above effects on fluid flow due to the intersection shapes of the angle patterns can be divided into two types: no effect for α ¼ 0� , 180� and 0� < β < 180� ; and an effect for 0� < α < 180� and 0� < β < 180� . Based on the different mechanisms of the different intersection shapes effects on fluid flow based on the experimental phenomena and the data analysis, the latter type can be further divided into three models, which are itemized in the α - β plane as model I, model II and model III (see Fig. 7). The three models are elaborated below.

Fig. 6. Definition of α and β as well as the effluent fractures. 5

L.F. Fan et al.

International Journal of Rock Mechanics and Mining Sciences 124 (2019) 104136

Fig. 7. Distribution of the angle patterns of specimens with images of special specimens.

(1) Model I

hindered by bulge_A at high Re. A curved surface boundary, which is reduced to a curvilinear boundary in the 2D schematic diagram, is formed by the combination of the hindrance from bulge_A, the inertial force, the viscous force, and the suddenly decreasing confining pressure owing to the sudden expansion of the fracture cross-sectional area. The curvilinear boundary is the envelope of the v direction on the boundary. In addition, the curvilinear boundary is the separatrix of the flow, which means that only a small flow can cross it. With J varying from 2.908 to 29.288, Re increases from 7587.135 to 22100.584, and v increases from 3.162 to 10.391 m/s. Due to the large average radius of the curvilinear boundary that can cross the mechanical aperture of fracture_2 at high v, only a small flow can enter fracture_2, in accordance with the experi­ mental phenomena. As a result, there is non-full flow in fracture_2. The flow in fracture_2 is pulled back by the atmospheric pressure because high v makes the pressure in fracture_2 decrease to less than the atmo­ spheric pressure according to the Bernoulli equation. In addition, v is not stable and it is easy to change the direction at high Re so that a huge eddy_1 (see Fig. 10(a)), which is a “white wandering eddy”, occurs in fracture_273. Eddy_1 makes the flow produce local friction, which leads to local energy loss. As a result, the outflow in fracture_2 moves out more difficultly. Bulge_A and bulge_B increase the value of q in advantageous fracture_3 but reduce the value of q in advantageous fracture_2. There­ fore, bulge_A and bulge_B produce values of ε larger than those in nat­ ural cases. Moreover, owing to the hindrance from bulge_A and the excessively long curvilinear boundary, eddy_2 occurs in fracture_3 (see

Because α ¼ 90� & β ¼ 30� (see Fig. 7(a)) results in the largest value of ε, it is selected as the example to analyse the effect mechanism. For α ¼ 90� & β ¼ 30� , ε increases from 113.75 to 145.89 with J varying from 2.908 to 29.288. As shown in the enlarged view of the experi­ mental phenomena in Fig. 9, the fluid flow almost jets out from outlet_3, while only a little fluid exits from oulet_2. There is evidently a “white wandering eddy” in fracture_3, the length of which is approximately 20 mm, and there is also a white “curved surface boundary” at the intersection. To analyse the mechanism expediently, streamlines to sketch the experimental phenomena for this angle pattern are presented. In Fig. 10(a), fracture_1 is regarded as the influent fracture, fracture_2 & 3 are regarded as the effluent fractures, and bulge_A, bulge_B, and cor­ ner_C are formed by two intersecting fractures. However, it is stressed that bulges and corners at the intersection in natural fractured rock masses are rounded or do not exist due to the washing and transporting action of fluid flow over a long time period. The unordered natural cases of “rounded bulges and corners (or not exist)” may cause the analysis results to not be unique. Therefore, ideal fractures with regular bulges and corners are analysed in this paper, and the results are compared with prediction results based on real situations in natural fractured rock masses. Furthermore, the experimental results based on ideal fractures can be used to generate prediction cases in the nature. The experimental phenomena show that the flow from fracture_1 is 6

L.F. Fan et al.

International Journal of Rock Mechanics and Mining Sciences 124 (2019) 104136

fracture_2, the curvilinear boundary bends towards fracture_2 because the fluid flow from fracture_1 is further hindered by bulge_A after being hindered by bulge_B. Hence, the portion of fluid flow from fracture_1 can enter fracture_2 more easily. The occurrence of eddy_4 is due to reflux, which is caused by the head-on hindrance from bulge_A (see Fig. 10(b)). Eddy_5 (see Fig. 10(b)) occurs in the shadow of bulge_B, where the in­ ertial force is smaller than that in the nearby positions. Similarly, eddy_6 (see Fig. 10(b)) occurs in fracture_1. Moreover, due to the flow running into the fracture wall and producing reflux, eddy_7 (see Fig. 10(b)) oc­ curs at corner_C. Finally, eddies narrow the effective flow path and consequently enhance the nonlinearity of the fluid flow.56 Bulge_A and bulge_B increase the value of q in disadvantageous fracture_2 but decrease the value of q in advantageous fracture_1. Therefore, bulge_A and bulge_B result in smaller values of ε than those in natural cases. The above effect mechanism is called model II, whose distributed area is shown in Fig. 7. Model II applies to the cases where the fluid flow is hindered by bulges at the fracture intersections, making advantageous fracture q smaller than would occur in natural cases and making disadvantageous fracture q larger than would occur in natural cases. Furthermore, the larger the area of the head-on-hindering bulge is, the greater the effect of the intersection shapes on the outlet flow rate distribution. (3) Model III As shown in Fig. 7(c), when fracture_2 is regarded as the influent fracture and fracture_1 & 3 are regarded as the effluent fractures, the case is α ¼ 157.5� & β ¼ 105� . Streamlines to sketch the experimental phenomena at this angle pattern are given in Fig. 10(c). With J varying from 2.908 to 29.632, Re increases from 7645.614 to 22423.860, v in­ creases from 3.556 to 10.430 m/s, and ε increases from 2.410 to 2.489. For high v, the portion of fluid flow from fracture_3 runs into the fracture wall sharply at corner_C and is then rebounded into fracture_1 and fracture_3. The inertial force does not have an effect owing to the angle of 157.5� between the influent orientation and advantageous fracture_1. In addition, v is not stable and it is easy to change the direction at high Re. Therefore, the portion of fluid flow is randomly assigned to both effluent fractures equally after a sharp rebound. Obviously, the sharp rebound of the fluid flow hindered by the fracture wall at corner_C in­ creases the value of q in disadvantageous fracture_3 but decreases the value of q in advantageous fracture_1. Therefore, a sharp fluid flow rebound can make the value of ε smaller than it is in the case of slight hindrance from the fracture wall at fracture intersections. Eddy_8 (see Fig. 10(c)) occurs in the shadow of bulge_B, where the inertial force is smaller than that in nearby positions. Likewise, eddy_9 (see Fig. 10(c)) occurs in fracture_1. Moreover, due to the flow running into the fracture wall to form reflux, eddy_10, eddy_11, eddy_12 and eddy_13 (see Fig. 10(c)) occur at corner_C. Finally, eddies narrow the effective flow path and consequently enhance the nonlinearity of the fluid flow.56 The above effect mechanism is called model III, whose distributed area is shown in Fig. 7. Model III applies to cases where a sharp rebound in the fluid flow hindered by the fracture wall at the fracture in­ tersections makes the advantageous fracture q smaller than that in the case of slight hindrance from the fracture wall at fracture intersections, and makes disadvantageous fracture q larger than that in the case of slight hindrance from the fracture wall at the fracture intersections. Furthermore, the sharper the flow rebound is, the greater the effect of the intersection shapes on the outlet flow rate distribution.

Fig. 8. Comparisons of the total and effluent flow rates for cases of three angle patterns (outlet_1 & 2, outlet_2 & 3, outlet_1 & 3), as well as comparisons of the experimental and cubic law results for the specimen whose image is shown in Fig. 7(e) (or Fig. 7(g)). a J is quadratically correlated with Q (or q) in all three angle patterns; b Outlet_1 & 3.

Fig. 10(a)). Due to the fluid flow running into the fracture wall, resulting in reflux, eddy_3 (see Fig. 10(a)) occurs at corner_C. Finally, eddies narrow the effective flow path and consequently enhance the nonline­ arity of the fluid flow.56 The above effect mechanism is called model I, whose distributed area is shown in Fig. 7. Model I applies to cases where the fluid flow is hin­ dered by bulges at the fracture intersections, making advantageous fracture q larger than it would be in natural cases and making disad­ vantageous fracture q smaller than it would be in natural cases. Furthermore, as α approaches 90� and β approaches 0� , the effects of the intersection shapes on the outlet flow rate distribution become increasingly significant. (2) Model II As shown in Fig. 7(b), when fracture_3 is regarded as the influent fracture and fracture_1 & 2 are regarded as the effluent fractures, the case is α ¼ 112.5� & β ¼ 75� . Streamlines to sketch the experimental phenomena under this angle pattern are shown in Fig. 10(b). With J varying from 2.046 to 29.115, Re increases from 6672.515 to 22519.298, v increases from 3.104 to 10.613 m/s, and ε increases from 2.201 to 2.290. Instead of crossing over the mechanical aperture of

3.3. The deviation of ε due to angle patterns compared with the prediction value in nature The relationship between ε and α with a J gradient for β ¼ 30� , 60� , 90� , 120� , and 150� is shown in Fig. 11. The relationship between ε and α shows few distinct changes in a narrow range of J, and that 7

L.F. Fan et al.

International Journal of Rock Mechanics and Mining Sciences 124 (2019) 104136

Fig. 9. The experimental phenomena for α ¼ 90� & β ¼ 30� .

demonstrates the reliability of this relationship. As shown in Fig. 11(a) (b)(c), for β � 90� , the average increasing rate of ε for α < 90� is slower than the average decreasing rate of ε for α > 90� . However, as shown in Fig. 11(d)(e), for β > 90� , the average increasing rate of ε for α < 90� is faster than the average decreasing rate of ε for α > 90� . As a result, the relationship between ε and α is an irregular symmetric distribution when β is fixed. Therefore, the comparison between the Gaussian distribution curve, which is symmetric, and the experimental curve, which is an irregular symmetric distribution, can be used to analyse the deviation of ε. In fact, the Gaussian distribution curves represent cases of rounded bulges and corners (or cases where bulges and corners not exist) or slight hindrance from the fracture wall at the fracture intersections, which is calculated by Eq. (10). The non-ignorable deviation of ε (Δε) due to α and β is shown in Fig. 11. Moreover, the detailed Δε due to α and β is shown in Fig. 12. For β ¼ 30� and J ¼ 17, the experimental curve and the Gaussian distribution curve of ε with the same symmetry axis are shown in Fig. 12 (a). The ε from the experiments is 7.318 (Δε1) less than the prediction value for α ¼ 15� because these angle patterns belong to model II. With increasing α, the angle patterns gradually transform into model I from model II, and the ε observed experimentally gradually exceeds the prediction values. For example, the experimental ε is 21.541 (Δε2) greater than the prediction value for α ¼ 75� . As α approaches 90� , the intersection shapes have a stronger effect on hindering fluid flow and Δε is larger. Δε reaches the maximum of 34.898 (Δε3) for α ¼ 90� . For α > 90� , as α decreases, the angle patterns gradually transform back into model II from model I, while the experimental ε gradually becomes smaller than the prediction value. For example, the experimental ε is 12.625 (Δε4) less than the prediction value for α ¼ 105� . Finally, as α decreases towards 180� , the angle patterns gradually transform into model III from model II, and the experimental ε remains smaller than the prediction value. For β ¼ 90� and J ¼ 17, the experimental curve and Gaussian distri­ bution curve of ε with the same symmetry axis are shown in Fig. 12(b). The experimental ε is 0.701 (Δε5) less than the prediction value for α ¼ 45� because these angle patterns belong to model II. As α increases, the angle patterns gradually transform into model I from model II, and the experimental ε gradually becomes larger than the prediction value. For example, the experimental ε is 1.139 (Δε6) greater than the pre­ diction value for α ¼ 90� . For α > 90� , as α decreases, the angle patterns gradually transform back into model II from model I, and the experi­ mental ε gradually becomes smaller than the prediction value. For

example, the experimental ε is 5.418 (Δε7) smaller than the prediction value for α ¼ 105� . Finally, as α decreases towards 180� , the angle pat­ terns gradually transform into model III from model II, while the experimental ε remains smaller than the prediction value. For β ¼ 150� and J ¼ 17, the experimental curve and the Gaussian distribution curve of ε with the same symmetry axis are shown in Fig. 12 (c). The experimental ε is 1.820 (Δε8) less than the prediction value for α ¼ 45� because these angle patterns belong to model II. As α increases towards 75� , the deviation degree of ε decreases. Δε reaches the mini­ mum of 0.146 (Δε9) for α ¼ 75� . For α > 75� , as α decreases, the angle patterns gradually transform into model III from model II, and the experimental ε remains smaller than the prediction value. For example, the experimental ε is 0.357 (Δε10) less than the prediction value at α ¼ 135� . Based on the above analysis, the Gaussian distribution curves do not accurately fit the ε - α experimental curves for β ¼ 30� , 60� , 90� , 120� , and 150� . Therefore, modified Gaussian distribution function fitting experimental results is provided: � �2 8 α μ 1 > > m1 2 > < He þ 1 ðα < μÞ ε¼ ; (11) � �2 > > > α μ 1 : m2 2 þ 1 ðα � μÞ He where H, μ, m1, and m2 are coefficients, which are provided in Table 1, for β ¼ 30� , 60� , 90� , 120� , and 150� . H determines the peak value of the function; μ, the symmetry axis, determines the horizontal shift of the function; and m1 and m2 determine the peak width of the function. The ε - β experimental curves for α ¼ 15� , 30� , 45� , 60� , 75� , 90� , 105� , 120� , 135� , 150� , and 165� are shown in Fig. 13(a). To make the growth trend of ε obvious, the ordinate range of the ε - β experimental curves is smaller than that of the ε - α experimental curves. The figure shows that ε continues increasing as β decreases whenever α is fixed. ε continues increasing sharply for β � 30� because as β decreases, the angle patterns gradually transform into model II from model III and then into model I, which results in an increase of ε. Furthermore, when it approaches the α - axis, the most significant hindrance from the bulges at the fracture intersection results in the largest deviation of ε. 3.4. Effects of the angle patterns on Δε The mechanism analysis based on the three models shows that the 8

L.F. Fan et al.

International Journal of Rock Mechanics and Mining Sciences 124 (2019) 104136

Fig. 10. Streamlines to sketch the experimental phenomena for a α ¼ 90� & β ¼ 30� ; b α ¼ 112.5� & β ¼ 75� ; c α ¼ 157.5� & β ¼ 105� .

9

L.F. Fan et al.

International Journal of Rock Mechanics and Mining Sciences 124 (2019) 104136

Fig. 11. Relationship between ε and α with a J gradient for β ¼ 30� , 60� , 90� , 120� , and 150� . a β ¼ 30� ; b β ¼ 60� ; c β ¼ 90� ; d β ¼ 120� ; e β ¼ 150� .

effects of the intersection shapes on the outlet flow rate distribution are different. The effects of the angle patterns on ε is shown in Fig. 14, where the blank circles represent angle patterns that make ε smaller than that in natural cases, the black circles represent angle patterns that make ε larger than that in natural cases, and the blank and black circles represent angle patterns that have no effect on ε. The distribution areas of the different circles correspond to the distribution areas of the three models. The sizes of the circles, which represent the magnitude of the effects of the angle patterns on ε, are based on the values of Δε. The magnitude of the effects of the model I angle patterns on ε is greater than that of the angle patterns of model II and model III (i.e., Δε2 ¼ 21.541, Δε3 ¼ 34.898, Δε1 ¼ 7.318, Δε4 ¼ 12.625, Δε7 ¼ 5.418, Δε8 ¼ 1.820, Δε10 ¼ 0.357). The results suggest that cases in which the effects of the angle patterns on Δε are larger are mainly distributed in β � 60� and 45� � α � 105� among the α - β plane. It can be inferred that as α approaches 90� and β approaches 0� , the effects of the angle patterns on Δε become increasingly significant because the nonlinear flow can completely flow into the advantageous fracture instead of the disadvantageous fracture when α ¼ 90� and β approaches 0� due to the hindrance from the bulges at the fracture intersection under model I angle patterns. In contrast, under model II and model III angle patterns, the nonlinear flow cannot flow into either of the effluent fractures. In fact, cases where β � 60� and 45� � α � 105� belong to model I, and the bulges have a significant effect on Δε. As shown in Fig. 12(b), the experimental ε of only 0.701 (Δε5) is less than the prediction value for α ¼ 45� & β ¼ 90� (see Fig. 7(f)), which belongs to model II. The reason is that under model II angle patterns, the overlap between the advantageous fracture and the extended orienta­ tion enhances the effect of the inertia force and weakens the hindrance from the bulges, which leads to small values of Δε. This is also why Δε reaches the minimum of 0.146 (Δε9) for α ¼ 75� & β ¼ 150� (see Fig. 7 (d)) in Fig. 12(c). As shown in Fig. 11, the fact that the experimental ε is close to the prediction value for α ¼ 15� & β ¼ 30� (see Fig. 7(h)), α ¼ 30� & β ¼ 60� (see Fig. 7(g)) and α ¼ 60� & β ¼ 120� (see Fig. 7(e)) confirms this explanation.

The angle patterns have no effect on ε for α ¼ 0� , 180� and 0� < β < 180� . 3.5. Effects of the angle patterns on the peak values of ε As shown in Fig. 13(b), the peak values of ε decrease (129.123, 10.784, 9.841, 5.128, and 4.546) for β ¼ 30� , 60� , 90� , 120� , 150� . As a result, the Gaussian distribution fitting curve changes from “concentrate type” to “discrete type”. For β � 60� , the peak values of ε increase slowly with decreasing β, whereas for β < 60� , the peak values of ε increase sharply with decreasing β, which leads to the “steep peak” for β ¼ 30� (see Fig. 13(b)). On the other hand, for different values of β, ε reaches its peak value at different values of α, i.e., for β � 90� , α ¼ 90� ; for β ¼ 120� , α ¼ 60� ; and for β ¼ 150� , α ¼ 75� . Therefore, the distribution of the peak values of ε is determined by the angle patterns, especially α. As shown in Fig. 14, the distribution fitting curve of the peak values of ε in the α - β plane is: � � � � β 120 β 120 þ1 (12) α ¼ 90 30 � exp exp 17:5 17:5 Based on the distribution fitting curve of the peak values of ε, for β � 90� , ε reaches its peak value at α ¼ 90� ; for 90� < β < 120� , as β in­ creases, the value of α where ε reaches its peak value decreases; for β ¼ 120� , the value of α for which ε reaches its peak value is at its minimum of 60� ; and for 120� < β < 180� , as β increases, the value of α for which ε reaches its peak value increases back towards 90� . The reason for this change is that for β � 90� , the shapes of the fracture in­ tersections have a substantial effect on ε, especially for α ¼ 120� , whereas for 90� < β < 180� , the effects of the fracture intersection shapes on ε decrease and the effects of the inertial force on ε increase. When the advantageous fracture and the extended orientation overlap, the inertial force has the strongest effect on ε. For example, for α ¼ 60� & β ¼ 120� (see Fig. 7(e)) and α ¼ 75� & β ¼ 150� (see Fig. 7(d)), at high Re, nonlinear flow can flow straight into the advantageous fracture with minimal effects of the fracture intersection shapes. Moreover, α goes 10

L.F. Fan et al.

International Journal of Rock Mechanics and Mining Sciences 124 (2019) 104136

Table 1 The parameters of modified Gaussian distribution function (Eq. (11)). β 30 60� 90� 120� 150� �

H

μ

m1

m2

128.123 9.784 8.841 4.128 3.546

90 90 90 60 75

24.354 32.095 22.167 20.491 22.400

6.207 13.496 8.220 44.666 44.944

Fig. 13. Experimental curves for a ε – β with J ¼ 17; b ε – α with J ¼ 17.

Fig. 12. Comparisons of the experimental and real results of the relationship between ε and α for a J ¼ 17, β ¼ 30� ; b J ¼ 17, β ¼ 90� ; c J ¼ 17, β ¼ 150� .

back towards 90� because the advantageous fracture and the extended orientation are nearly overlapped for the cases where β approaches 180� and α approaches 90� .

11

L.F. Fan et al.

International Journal of Rock Mechanics and Mining Sciences 124 (2019) 104136

shows that ε is larger for cases where 45� � α � 105� and β � 60� . This area corresponds to the larger upheaval of the unclosed curved surface of ε for 45� � α � 105� . Furthermore, it can be inferred that ε went to infinity for the cases where α ¼ 90� and β approaches 0� . 4. Conclusions In this paper, a series of experiments including fifty-one angle pat­ terns with varying hydraulic gradients (J ¼ 2, 7, 12, 17, 22, and 30) are conducted on one-inlet-two-outlet specimens at high Reynolds numbers. Based on the experimental results, for high hydraulic gradients and high Reynolds numbers, fluid flow among fractures with intersections shows strong nonlinearity. In addition, the ideal intersection shapes, with regular bulges and corners, promote the occurrences of eddies by hin­ dering fluid flow, which can narrow the effective flow path and conse­ quently enhance the nonlinearity of the fluid flow. Furthermore, the angle patterns due to the shapes of the fracture intersections affect ε. The effects of the angle patterns on the outlet flow rate distribution can be divided into three models based on the different mechanisms of the ef­ fects of the intersection shape on the fluid flow: model I, model II and model III. Furthermore, the modified Gaussian distribution function is formulated to fit the ε – α experimental curves well for β ¼ 30� , 60� , 90� , 120� , and 150� . The distribution fitting curve of the ε peak values is provided in the α - β plane. The main conclusions are as follows. For J (100 - 101) and Re (104 - 105), J is quadratically correlated with Q (or q). The intersection shapes of model I angle patterns make ε larger than that in natural cases. However, the intersection shapes of model II and

Fig. 14. The influence of α and β on ε (the size of circle represents the influ­ ence degree).

3.6. The relationship between ε, α and β The 3D relationship between ε, α and β is shown in Fig. 15 for e ¼ 2.108–2.126 mm, l ¼ 40 mm and J ¼ 17. The unclosed curved sur­ face of ε is an irregular symmetric distribution about α ¼ 90� with a larger upheaval for 45� � α � 105� , especially for α ¼ 90� .The figure

Fig. 15. 3D relationship between ε, α and β for J ¼ 17. 12

L.F. Fan et al.

International Journal of Rock Mechanics and Mining Sciences 124 (2019) 104136

model III angle patterns make ε smaller than that in natural cases. The effect of the angle patterns of model I on Δε is larger than that of the angle patterns of model II and model III. For β � 90� , the value of α where ε reaches its peak value is 90� ; for 90� < β < 120� , as β increases, the value of α for which ε reaches its peak value decreases; for β ¼ 120� , the value of α for which ε reaches its peak value is at its minimum of 60� ; and for 120� < β < 180� , as β increases, the value of α for which ε reaches its peak value increases back towards 90� . The relationship between ε and α is an irregular symmetric distri­ bution for β ¼ 30� , 60� , 90� , 120� , and 150� , which can be fitted well by the modified Gaussian distribution function. For 45� � α � 105� and β � 60� , ε and Δε are both larger than in other cases among the α - β plane. Moreover, for cases where α approaches 90� and β approaches 0� , the angle patterns have an increasingly significant effect on ε and Δε. However, for cases where α ¼ 0� , 180� and 0� < β < 180� , the angle patterns have no effect on ε. The specific geometries caused by the real intersection shapes with rounded bulges and corners can strongly influence the nonlinear flow behaviour and the outlet flow rate distribution. Further experimental studies are required to fully elucidate the nonlinear flow behaviour and the outlet flow rate distribution for a broad spectrum of intersection geometries. Moreover, the effects of fracture length, intersection angle shape and the roughness of fracture surfaces will also be considered in the future study.

14 Chen Y, Ma GW, Wang HD, Li T. Evaluation of geothermal development in fractured hot dry rock based on three dimensional unified pipe-network method. Appl Therm Eng. 2018;136:219–228. 15 Ranjith PG, Darlington W. Nonlinear single-phase flow in real rock joints. Water Resour Res. 2007;43(9):146–156. 16 Chen YF, Liu MM, Hu SH, Zhou CB. Non-Darcy’s law-based analytical models for data interpretation of high-pressure packer tests in fractured rocks. Eng Geol. 2015; 199:91–106. 17 Chen YF, Zhou JQ, Hu SH, Hu R, Zhou CB. Evaluation of Forchheimer equation coefficients for non-Darcy flow in deformable rough-walled fractures. J Hydrol. 2015; 529:993–1006. 18 Kosakowski G, Berkowitz B. Flow pattern variability in natural fracture intersections. Geophys Res Lett. 1999;26(12):1765–1768. 19 Li B, Liu RC, Jiang YJ. Influences of hydraulic gradient, surface roughness, intersecting angle, and scale effect on nonlinear flow behavior at single fracture intersections. J Hydrol. 2016;538:440–453. 20 Liu RC, Zhu TT, Jiang YJ, et al. A predictive model correlating permeability to twodimensional fracture network parameters. Bull Eng Geol Environ. 2018. https://doi. org/10.1007/s10064-018-1231-8. 21 Hull LC, Koslow KN. Streamline routing through fracture junctions. Water Resour Res. 1986;22(12):1731–1734. 22 Lomize GM. Flow in Fractured Rocks (In Russian). Moscow: Gosenergoizdat; 1951. 23 Louis C. A study of groundwater flow in jointed rock and its influence on the stability of rock masses. Rock Mech Res Rep. 1969;10:1–90. 24 Iwai K. Fundamental Studies of Fluid Flow through a Single Fracture. Berkeley: Dissertation, University of California; 1976. 25 Tsang YW, Witherspoon PA. The dependence of fracture mechanical and fluid flow properties on fracture roughness and sample size. J Geophys Res. 1983;88(B3): 2359–2366. 26 Tsang YW. The effect of tortuosity on fluid flow through a single fracture. Water Resour Res. 1984;20(9):1209–1215. 27 Tsang YW, Tsang CF. Channel model of flow through fracture media. Water Resour Res. 1987;23(3):467–479. 28 Barton N, Bandis S, Bakhtaret K. Strength, deformation and conductivity coupling of rock joints. Int J Rock Mech Min Sci Geomech Abstr. 1985;22(3):121–140. 29 Barton N, Quadros EFD. Joint aperture and roughness in the prediction of flow and groutability of rock masses. Int J Rock Mech Min Sci. 1997;34(3-4), 252.e1-252.e14. 30 Brown SR. Fluid flow through rock joints: the effect of surface roughness. J Geophys Res. 1987;92(B2):1337–1347. 31 Lee YH, Carr JR, Barr DJ, Haas CJ. The fractal dimension as a measure of the roughness of rock discontinuity profiles. Int J Rock Mech Min Sci. 1990;27(6): 395–404. 32 Zimmerman RW, Al-Yaarubi AA, Pain CC, Grattoni CA. Nonlinear regimes of fluid flow in rock fractures. Int J Rock Mech Min Sci. 2004;41(3):163–169. 33 Klimczak C, Schultz RA, Parashar R, Reeves DM. Cubic law with aperture-length correlation: implications for network scale fluid flow. Hydrogeol J. 2010;18(4): 851–862. 34 Zhang Z, Nemcik J, Ma S. Micro-and macro-behaviour of fluid flow through rock fractures: an experimental study. Hydrogeol J. 2013;21(8):1717–1729. 35 Scesi L, Gattinoni P. Roughness control on hydraulic conductivity in fractured rocks. Hydrogeol J. 2007;15(2):201–211. 36 Ju Y, Zhang QG, Yang YM, Xie HP, Gao F, Wang HJ. An experimental investigation on the mechanism of fluid flow through single rough fracture of rock. Sci China Technol Sci. 2013;56(8):2070–2080. 37 Olsson WA, Brown SR. Hydromechanical response of a fracture undergoing compression and shear. Int J Rock Mech Min Sci. 1993;30(7):845–851. 38 Olsson R, Barton N. An improved model for hydromechanical coupling during shearing of rock joints. Int J Rock Mech Min Sci. 1993;38(3):317–329. 39 Hans J, Boulon M. A new device for investigating the hydro-mechanical properties of rock joints. Int J Numer Anal Methods Geomech. 2003;27(6):513–548. 40 Auradou H, Drazer G, Hulin JP, Koplik J. Permeability anisotropy induced by the shear displacement of rough fracture walls. Water Resour Res. 2005;41(9):477–487. 41 Auradou H, Drazer G, Boschan A, Hulin JP, Koplik J. Flow channeling in a single fracture induced by shear displacement. Geothermics. 2006;35(5-6):576–588. 42 Koyama T, Fardin N, Jing LR, Stephansson O. Numerical simulation of shear-induced flow anisotropy and scale-dependent aperture and transmissivity evolution of rock fracture replicas. Int J Rock Mech Min Sci. 2006;43(1):89–106. 43 Koyama T, Li B, Jiang Y, Jing LR. Numerical modelling of fluid flow tests in a rock fracture with a special algorithm for contact areas. Comput Geotech. 2009;36(1-2): 291–303. 44 Li B, Jiang Y, Koyama T, Jing LR, Tanabashi Y. Experimental study of the hydromechanical behavior of rock joints using a parallel-plate model containing contact areas and artificial fractures. Int J Rock Mech Min Sci. 2008;45(3):362–375. 45 Giger S, Clennell MB, Harbers C, et al. Design, operation and validation of a new fluid-sealed direct shear apparatus capable of monitoring fault-related fluid flow to large displacements. Int J Rock Mech Min Sci. 2011;48(7):1160–1172. 46 Chen YF, Zhou JQ, Hu SH, Hu R, Zhou CB. Evaluation of Forchheimer equation coefficients for non-Darcy flow in deformable rough-walled fractures. J. Hydrol. 2015;529:993–1006. 47 Zhou JQ, Hu SH, Fang S, Chen YF, Zhou CB. 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–218. 48 Wilson CR, Witherspoon PA. Flow interference effects at fracture intersections. Water Resour Res. 1976;12(1):102–104. 49 Tian KM. The hydraulic properties of crossing-flow in an intersected fracture. Geol Rev. 1983;60(2):202–214.

Declaration of competing interest No Conflict of Interest. Acknowledgements The research work is supported by the National Natural Science Foundation of China (Grant No. 41772309), National Key Research and Development Program of China (2017YFC1501302). References 1 Tao QF. Numerical Modeling of Fracture Permeability Change in Naturally Fractured Reservoirs Using a Fully Coupled Displacement Discontinuity Method. Dissertation. Texas: A&M University; 2010. 2 Chen Y, Ma GW, Wang HD. Heat extraction mechanism in a geothermal reservoir with rough-walled fracture networks. Int J Heat Mass Transf. 2018;126:1083–1093. 3 Chen Y, Ma GW, Wang HD. The simulation of thermo-hydro-chemical coupled heat extraction process in fractured geothermal reservoir. Appl Therm Eng. 2018;143: 859–870. 4 Ma GW, Chen Y, Jin Y, Wang HD. Modelling temperature-influenced acidizing process in fractured carbonate rocks. Int J Rock Mech Min. 2018;105:73–84. 5 Fan LF, Wu ZJ, Wan Z, Gao JW. Experimental investigation of thermal effects on dynamic behavior of granite. Appl Therm Eng. 2017;125:94–103. 6 Jiao YY, Zhang XL, Zhang HQ, Li HB, Yang SQ, Li JC. A coupled thermo-mechanical discontinuum model for simulating rock cracking induced by temperature stresses. Comput Geotech. 2015;67:142–149. 7 Jiang QH, Ye ZY, Zhou CB. A numerical procedure for transient free surface seepage through fracture networks. J Hydrol. 2014;519:881–891. 8 Ma GW, Wang HD, Fan LF, Wang B. Simulation of two-phase flow in horizontal fracture networks with numerical manifold method. Adv Water Resour. 2017;108: 293–309. 9 Ma GW, Wang HD, Fan LF, Chen Y. Segmented two-phase flow analysis in fractured geological medium based on the numerical manifold method. Adv Water Resour. 2018;121:112–129. 10 Dai F, Wei MD, Xu NW, Ma Y, Yang DS. Numerical assessment of the progressive rock fracture mechanism of cracked chevron notched Brazilian discspecimens. Rock Mech Rock Eng. 2015;48(2):463–479. 11 Ma GW, Wang HD, Fan LF, Chen Y. A unified pipe-network-based numerical manifold method for simulating immiscible two-phase flow in geological media. J. Hydrol. 2019;568:119–134. 12 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–1804. 13 Singh KK, Singh DN, Ranjith PG. Laboratory simulation of flow through single fractured granite. Rock Mech Rock Eng. 2015;48(3):987–1000.

13

L.F. Fan et al.

International Journal of Rock Mechanics and Mining Sciences 124 (2019) 104136 61 Witherspoon PA, Wang JSY, Iwai K, Gale JE. Validity of cubic law for fluid flow in a deformable rock fracture. Water Resour Res. 1980;16(6):1016–1024. 62 Bear J. Dynamics of Fluids in Porous Media. New York: Am. Elsevier; 1972. 63 Jung R. Hydraulic in situ investigations of an artificial fracture in the Falkenberg granite. Int J Rock Mech Min Sci. 1989;26(3-4):301–308. 64 Wen Z, Huang GH, Zhan HB. Non-Darcian flow in a single confined vertical fracture toward a well. J Hydrol. 2006;330(3-4):698–708. 65 Cherubini C, Giasi CI, Pastore N. Bench scale laboratory tests to analyze non-linear flow in fractured media. Hydrol Earth Syst Sci Discuss. 2012;9(4):5575–5609. 66 Javadi M, Sharifzadeh M, Shahriar K. A new geometrical model for non-linear fluid flow through rough fractures. J Hydrol. 2010;389(1-2):1830. 67 Whitaker S. The Forchheimer equation: a theoretical development. Transp Porous Media. 1996;25(1):27–61. 68 Holditch SA, Morse RA. The effects of non-Darcy flow on the behavior of hydraulically fractured gas wells. J Pet Technol. 1976;28(10):1169–1179. 69 Kohl T, Evans KF, Hopkirk RJ, Jung R, Rybach L. Observation and simulation of nonDarcian flow transients in fractured rock. Water Resour Res. 1997;33(3):407–418. 70 Moutsopoulos KN. Exact and approximate analytical solutions for unsteady fully developed turbulent flow in porous media and fractures for time dependent boundary conditions. J Hydrol. 2009;369(1-2):78–89. 71 Qian JZ, Zhan HB, Chen Z, Ye H. Experimental study of solute transport under nonDarcian flow in a single fracture. J Hydrol. 2011;399(3):246–254. 72 Liu RC, Li B, Jiang YJ, Huang N. Review: mathematical expressions for estimating equivalent permeability of rock fracture networks. Hydrogeol J. 2016;24:1623–1649. 73 Wu ZJ, Ma HY, Zhou MD. Vorticity and Vortex Dynamics. Berlin: Springer; 2006.

50 Philip JR. The fluid mechanics of fracture and other junctions. Water Resour Res. 1988;24(2):239–246. 51 Johnson J, Brown S. Experimental mixing variability in intersecting natural fractures. Geophys Res Lett. 2001;28(22):4303–4306. 52 Johnson J, Brown S, Stockman H. Fluid flow and mixing in rough-walled fracture intersections. J Geophys Res. 2006;111(B12):97–112. 53 Hu YJ, Mao GH, Cheng WP, Zhang JJ. Theoretical and experimental study on flow distribution at fracture intersections. J Hydraul Res. 2005;43(3):321–327. 54 Ji S, Nicholl MJ, Glass RJ, Lee K. Influence of simple fracture intersections with differing aperture on density-driven immiscible flow: wetting versus nonwetting flows. Water Resour Res. 2006;42(10):730–732. 55 Zafarani A, Detwiler RL. An efficient time-domain approach for simulating Pedependent transport through fracture intersections. Adv Water Resour. 2013;53(1): 198–207. 56 Liu RC, Li B, Jiang YJ. Critical hydraulic gradient for nonlinear flow through rock fracture networks: the roles of aperture, surface roughness, and number of intersections. Adv Water Resour. 2016;88:53–65. 57 Neuville A, Flekkøy EG, Toussaint R. Influence of asperities on fluid and thermal flow in a fracture: a coupled lattice Boltzmann study. J Geophys Res. 2013;118(7): 3394–3407. 58 Xie LZ, Gao C, RenL, Li CB. Numerical investigation of geometrical and hydraulic properties in a single rock fracture during shear displacement with the Navier–Stokes equations. Environ Earth Sci. 2015;73(11):7061–7074. 59 Zou LC, Jing LR, Cvetkovic V. Roughness decomposition and nonlinear fluid flow in a single rock fracture. Int J Rock Mech Min Sci. 2015;75:102–118. 60 Xiong X, Li B, Jiang YJ, Koyama T, Zhang CH. Experimental and numerical study of the geometrical and hydraulic characteristics of a single rock fracture during shear. Int J Rock Mech Min Sci. 2011;48(8):1292–1302.

14