# Effects of Hydraulic Gradient, Intersecting Angle, Aperture, and Fracture Length on the Nonlinearity of Fluid Flow in Smooth Intersecting Fractures: An Experimental Investigation

Authors: Zhijun Wu; Lifeng Fan; Shihe Zhao
Journal: Geofluids
 (2018)

---

## Abstract

This study experimentally investigated the nonlinearity of fluid flow in smooth intersecting fractures with a high Reynolds number and high hydraulic gradient. A series of fluid flow tests were conducted on one-inlet-two-outlet fracture patterns with a single intersection. During the experimental tests, the syringe pressure gradient was controlled and varied within the range of 0.20–1.80 MPa/m. Since the syringe pump used in the tests provided a stable flow rate for each hydraulic gradient, the effects of hydraulic gradient, intersecting angle, aperture, and fracture length on the nonlinearities of fluid flow have been analysed for both effluent fractures. The results showed that as the hydraulic gradient or aperture increases, the nonlinearities of fluid flow in both the effluent fractures and the influent fracture increase. However, the nonlinearity of fluid flow in one effluent fracture decreased with increasing intersecting angle or increasing fracture length, as the nonlinearity of fluid flow in the other effluent fracture simultaneously increased. In addition, the nonlinearities of fluid flow in each of the effluent fractures exceed that of the influent fracture.


---

## Body

## 1. Introduction

The nonlinear behaviour of fluid flow in natural fractured rock masses is a critical hotspot in numerous branches of geoscience and rock engineering fields, such as geological disposal of radioactive waste [1, 2], reservoir storage [3, 4], underground mining [5, 6], and geothermal extraction [7, 8]. Considering that an intact rock matrix in a deep formation has a very low permeability (which is usually assumed to be impermeable) [9–11], the nonlinear behaviour of fluid flow in fractured rock masses is heavily controlled by the flow behaviour in a single fracture or fracture network. In recent decades, the effects of fracture and hydraulic characteristics, such as fracture length, fracture density, aperture (or hydraulic aperture), fracture orientation, fracture connectivity ratio, hydraulic gradient, intersecting angle, surface roughness, scale effect, and number of intersections and dead-ends, on the nonlinear behaviour of fluid flow in rock fractures have been extensively and systematically studied, forming the basis for understanding nonlinear fluid flow behaviour in natural fractured rock masses. However, due to the complexity of fracture distribution and fracture characteristics in naturally fractured rock masses, clearly determining the effects of all the fracture and hydraulic parameters on the nonlinear behaviour of fluid flow in naturally fractured rock masses is still challenging.Many researchers have focused on the effects of the fracture and hydraulic characteristics on the fluid flow behaviour in a single fracture, intersecting fractures, and fracture networks. Long and Billaux [10] presented a fracture network model that accounted for the observed spatial variability by generating a network in subregions, where the properties of each subregion were predicted through geostatistics. They found that approximately 0.1% of the fractures primarily controlled the permeability of the system. De Dreuzy et al. [11, 12] conducted a numerical and theoretical study on the permeability variation by assuming that the fracture length obeys a power law distribution and concluded that the hydraulic properties of fracture networks with a power law length distribution can be classified into three simplified types. When the power law exponent a>3, fracture networks essentially consist of very small fractures and the percolation theory applies. On the other hand, when a<2, flow is mostly channelled into longer fractures and fracture networks can be considered as a superposition of long fractures. Between these ranges (2<a<3), fracture length distributions cannot be restricted to a unique length, even though the smaller fractures have a small contribution to flow. Olson [13] focused on a nonlinear aperture-to-length relationship and suggested that fracture apertures scale with their lengths to the 1/2 power; this result was obtained by using linear elastic fracture mechanics in a homogeneous body, with subcritical and critical (equilibrium law) fracture propagation criteria. In addition, Olson [13] determined that a fracture aspect ratio (aperture/length) decreases with increasing fracture length to the negative 1/2 power. Min and Jing [14] conducted a series of numerical simulations of the mechanical deformation of fractured rock masses at different scales with many realizations of Discrete Fracture Networks (DFNs) and concluded that a representative elementary volume can be defined and the elastic properties of the fractured rock mass can be approximately represented by the elastic compliance tensor. Min et al. [15] investigated the stress-dependent permeability issue in fractured rock masses, considering the effects of nonlinear normal deformation and shear dilation of fractures. They found that the permeability of fractured rocks decreases with increasing stress magnitudes when the stress ratio is not sufficiently large to cause shear dilation of the fractures, whereas the permeability increases with increasing stress when the stress ratio is sufficiently large. In addition, permeability changes at low stress levels are more substantial than those at high stress levels due to the nonlinear relation between fracture normal stress and displacement. Based on a newly developed correlation equation, Baghbanan and Jing [16] investigated the permeability of fractured rocks by considering the correlation between the distributions of aperture and trace length of the fractures. Their results showed that there is a significant difference between correlated and noncorrelated apertures and fracture length distributions, which demonstrated that the hydromechanical behaviour of fractured rocks is considerably scale- and stress-dependent when the aperture and fracture length distributions are correlated. By solving the Navier-Stokes equations without linearization, using a self-developed 2D finite volume method, Zou et al. [17] investigated the effects of wall surface roughness on fluid flow through rock fractures. Their results indicated that even with the same total flow rate, the flow patterns and velocity fields have significant differences, which are caused by the secondary roughness and can even create time-dependent dynamic flow, with moving and changing eddies, when the Reynolds number (Re) is high. Furthermore, Zou et al. [18, 19] studied the fluid flow and solute transport in a 3D rock fracture-matrix system, which had two rough-walled fractures with an orthogonal intersection. Zhou et al. [20] experimentally investigated the nonlinear flow characteristics of fluid flow through the rough-walled fractures subjected to a wide range of confining pressures (1.0–30.0 MPa/m) at low Re. They found that the obtained critical Reynolds number versus confining pressure curves generally displays a nonlinear weakening stage (I) in the early stage of confining pressure loading, which is followed by a nonlinear enhancement stage (II) as the confining pressure further increases. Using numerical simulations based on DFNs, Liu et al. [21, 22] estimated the effects of fracture intersections and dead-ends on nonlinear flow and particle transport in 2D DFNs and demonstrated that wider fracture apertures, rougher fracture surfaces, and greater numbers of fracture intersections in a DFN may result in the onset of nonlinear flow at a lower critical hydraulic gradient. In addition, they found that the effects of fracture dead-ends on fluid flow are negligible (<1.5%); however, fracture dead-ends have a strong impact on the breakthrough curves of particles in DFNs with a relative time deviation rate in the range of 5–35%. Based on fluid flow tests and numerical simulations, Li et al. [23] found that the nonlinear degree of hydraulic gradient (J) versus flow rate (Q) increases with J and the joint roughness coefficient (JRC), whereas the nonlinear degree of hydraulic gradient versus flow rate increases as the radius of the truncating circle decreases (the radius is equivalent to fracture length). In addition, they found that the intersecting angle affects the effluent fracture flow rate (q1 and q2 stand for two effluent fractures) for three different intersecting angle patterns.Previous studies show that a numerical model based on discrete fracture networks (DFNs), which usually consists of thousands of fracture elements and nodes [14, 21, 22], was mainly used to investigate the nonlinear fluid flow in fractured rock masses. However, a reasonable description of the fluid flow in the basic elements and nodes is the key prerequisite to ensure the accuracy of the numerical results. Reasonable modelling of nonlinear fluid flow behaviour through an entire fractured rock mass depends on realistically describing the nonlinear fluid flow in both a single fracture and an intersecting fracture network. However, the fluid flow through one-inlet-two-outlet fracture patterns with a single intersection, which is a fundamental element of DFNs, has not yet been extensively studied. Though a few studies have been carried out by numerical modelling [23] to investigate the effects of parameters such as J, intersecting angle (θ), aperture (e), and fracture length (Rr) on the nonlinear behaviour of fluid flow, the reliability of the numerical results is low without experimental validation. In this study, a fluid flow test system was built to conduct a series of laboratory tests on one-inlet-two-outlet specimens with particular attention paid to the effects of J, θ, e, andRr on the nonlinearity of fluid flow in both effluent fractures at large Re and large J. During the tests, the syringe pressure gradient was controlled and varied in the range of 0.20–1.80 MPa/m, resulting in a high flow velocity at the inlet. The outlet of each effluent was open, where the discharged fluid was collected in a bucket and later measured by an electronic balance. Based on the collected fluid weight, the flux of each effluent was calculated at the end of the experiment by dividing the effluent quantity by the collecting time.


## 2. Methodology

### 2.1. Theoretical Background

#### 2.1.1. Linear Darcy Flow Zone

The simple parallel plate model is the only fracture model available to calculate the hydraulic conductivity of a fracture. This model assumes a steady incompressible Newtonian fluid flow in a single fracture under a one-dimensional pressure gradient between two smooth parallel plates separated by an aperturee [24, 25]. At sufficiently low Q, this calculation yields the well-known “cubic law” [26, 27] and “Darcy’s law” [28, 29].(1)Q=-e3w12μ∇P=-kAμ∇P,where Q is the total volumetric flow rate, which represents the flow rate of the influent fracture in this study, ∇P is the pressure gradient, w is the width of the fracture, e is the uniform aperture of the idealized smooth fracture, μ is the dynamic viscosity,k is the intrinsic permeability defined as e2/12[5, 27], and A is the cross-sectional area equal to ew. Equation (1) predicts that Q is proportional to the cube of e and that Q is linearly correlated with ∇P when Re is sufficiently small. The cubic law has been widely applied to rough-walled fractures in rock; however, in these applications, e is replaced by the so-called hydraulic aperture (eH) [24, 30, 31]. In this paper, eH is equal to e since the walls of the fractures in the test specimen are smooth.


#### 2.1.2. Non-Darcy Flow Zone

Darcy’s law and the cubic law are only valid when the inertial force is negligible, compared with the magnitude of the viscous force, and this condition is anticipated only whenQ is low. Nonlinear flow occurs when, with an incremental increase in Re, Q increases more than a proportional incremental amount [26, 32–35]. Forchheimer’s law [8, 36–38] has been most widely used to describe the nonlinear flow in fractures and porous media, especially in strong inertial regimes:(2)-∇P=a′Q+b′Q2,where a′ and b′ are coefficients, a′=12μ/we3 and b′=βρ/w2e2, ρ is the fluid density, and β is called either the non-Darcy coefficient or the inertial resistance [5, 39, 40].The hydraulic gradientJ is proportional to ∇P:(3)J=h1-h2ΔL=1ρgP1-P2ΔL=1ρg∇P∇L=1ρg∇P,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. Equation (2) can be written in a Forchheimer’s law form as follows:(4)J=aQ+bQ2.The term aQ represents the linear gradient, and the term bQ2 represents the nonlinear gradient.Coefficientsa and b in (4) are commonly written as follows:(5a)a=μρgkA=12μρgwe3,(5b)b=βgA2=βgw2e2.Equation (5a) states that a is linearly proportional to the reciprocal of the cube of e. If the values of w and e are measurable and constant, the value of a can be calculated directly. In addition, in some studies, the value of a is considered the intrinsic permeability of the fractured media [29, 41, 42]. The coefficient b in (5b) is determined from the geometrical properties of the fractures and media in the experiments [34, 40, 43].For fluid flow within fractures, the ratio of the inertial force to the viscous force is defined as the Reynolds number (Re), which is expressed as follows [44, 45]: (6)Re=ρveμ=ρQμw,where v is the bulk flow velocity in the fractures. Critical values of Re can distinguish between the linear and nonlinear flow conditions [20, 22, 24, 31, 46, 47].The vector flow velocityv can be reduced to a one-dimensional flow velocity along a fracture intersection:(7)v=QA.


### 2.2. Experimental Setup

#### 2.2.1. Specimen Preparation

First, specimens containing one inlet and two outlets were designed by AutoCAD, which can accurately create a regular single intersection of three fractures. Based on the control variate method, the specimens were designed withθ = 30°, 60°, 90°, 120°, and 150°, e = 2, 3, 5, 10, and 20 mm, andRr = 40, 100, 200, 300, and 400 mm. The specimens were made with PMMA (polymethyl methacrylate) material, which is transparent and easy to process. Then, the specimen design diagrams, in EPS format, were transferred to the computer of the dedicated PMMA machine. Next, the specimens were created at a series of scales with a small error: ±0.02 mm. The sizes of e,Rr, and w were precisely measured by a digital Vernier caliper. Note that the 5.7 mm width w was controlled by the length of the drill employed to cut the PMMA. The error of w was ±0.1 mm (i.e., w = 5.648–5.745 mm) because the surface of the initial PMMA slab was not flat. Next, a high-pressure flame gun was used to polish the fracture walls, so that the initial condition of smooth fracture walls was met in this study. Figure 1 shows the pictures of the specimens. Finally, the specimens were glued to plates with the same side area as the specimens using the PMMA adhesive. This procedure was repeated to ensure absolute sealing of the fractures.Figure 1
The pictures of specimens for the cases ofθ = 30°, 60°, 90°, 120°, and 150°.
(a)
θ = 30° (b)
θ = 60° (c)
θ = 90° (d)
θ = 120° (e)
θ = 150°


#### 2.2.2. Test System and Experimental Procedure

A schematic view of the fluid flow test system is shown in Figure2(a). Filtered tap water was introduced into the inlet of the specimen by a high-pressure syringe pump; the syringe pressure had an accuracy of ±0.01 MPa/m. The syringe pressure gradient range of 0.20–1.80 MPa/m resulted in the range of 2.00–9.00 L/min (3.33 × 10−5–1.50 × 10−4 m3/s) for the syringe flow rates. A digital pressure gauge with an accuracy of ±0.1 kPa/m was used to measure the pressure between the inlet and the outlets. The specimens were placed on a horizontal platform fixed by transparent tape so that the gravitational pressure difference between the inlet and the outlets was negligible. Two buckets, each of which had a large hole in the side, were employed to collect the discharged water. Electronic scales with an accuracy of ±0.02 g were used to measure the weight of the effluent. Fluid flow through the transparent PMMA specimen was recorded by a SLR (single lens reflex) camera mounted above the experimental setup. To connect the tube of the syringe pump to the inflow tank of the model, at the end of each fracture, a square slot with a size of 5 × 8 × 8 mm was left to fit to be connected with the pressure-proof PVC tubule of a 6 mm inner diameter and 8 mm outer diameter. The connection between the PVC tubule and the inflow tank of the model was then sealed using a hot glue gun as shown in Figure 2(b) in the revised manuscript. In addition, the other end of the PVC inlet tubule of the inflow tank was quickly connected with the tube of syringe pump via an inflow butt joint. The entire flow path, including the connecting tubes, inflow butt joint, and all three fractures, was sealed and watertight. Filtered tap water was used as the fluid for the experiments, with a density of 998.2 kg/m3 and a dynamic viscosity μ of 0.001 Pa·s at a room temperature of 20°C.Figure 2
(a) The fluid flow test system; (b) the connection between the PVC tubule and the inflow tank of the model.
(a) (b)


## 2.1. Theoretical Background

### 2.1.1. Linear Darcy Flow Zone

The simple parallel plate model is the only fracture model available to calculate the hydraulic conductivity of a fracture. This model assumes a steady incompressible Newtonian fluid flow in a single fracture under a one-dimensional pressure gradient between two smooth parallel plates separated by an aperturee [24, 25]. At sufficiently low Q, this calculation yields the well-known “cubic law” [26, 27] and “Darcy’s law” [28, 29].(1)Q=-e3w12μ∇P=-kAμ∇P,where Q is the total volumetric flow rate, which represents the flow rate of the influent fracture in this study, ∇P is the pressure gradient, w is the width of the fracture, e is the uniform aperture of the idealized smooth fracture, μ is the dynamic viscosity,k is the intrinsic permeability defined as e2/12[5, 27], and A is the cross-sectional area equal to ew. Equation (1) predicts that Q is proportional to the cube of e and that Q is linearly correlated with ∇P when Re is sufficiently small. The cubic law has been widely applied to rough-walled fractures in rock; however, in these applications, e is replaced by the so-called hydraulic aperture (eH) [24, 30, 31]. In this paper, eH is equal to e since the walls of the fractures in the test specimen are smooth.


### 2.1.2. Non-Darcy Flow Zone

Darcy’s law and the cubic law are only valid when the inertial force is negligible, compared with the magnitude of the viscous force, and this condition is anticipated only whenQ is low. Nonlinear flow occurs when, with an incremental increase in Re, Q increases more than a proportional incremental amount [26, 32–35]. Forchheimer’s law [8, 36–38] has been most widely used to describe the nonlinear flow in fractures and porous media, especially in strong inertial regimes:(2)-∇P=a′Q+b′Q2,where a′ and b′ are coefficients, a′=12μ/we3 and b′=βρ/w2e2, ρ is the fluid density, and β is called either the non-Darcy coefficient or the inertial resistance [5, 39, 40].The hydraulic gradientJ is proportional to ∇P:(3)J=h1-h2ΔL=1ρgP1-P2ΔL=1ρg∇P∇L=1ρg∇P,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. Equation (2) can be written in a Forchheimer’s law form as follows:(4)J=aQ+bQ2.The term aQ represents the linear gradient, and the term bQ2 represents the nonlinear gradient.Coefficientsa and b in (4) are commonly written as follows:(5a)a=μρgkA=12μρgwe3,(5b)b=βgA2=βgw2e2.Equation (5a) states that a is linearly proportional to the reciprocal of the cube of e. If the values of w and e are measurable and constant, the value of a can be calculated directly. In addition, in some studies, the value of a is considered the intrinsic permeability of the fractured media [29, 41, 42]. The coefficient b in (5b) is determined from the geometrical properties of the fractures and media in the experiments [34, 40, 43].For fluid flow within fractures, the ratio of the inertial force to the viscous force is defined as the Reynolds number (Re), which is expressed as follows [44, 45]: (6)Re=ρveμ=ρQμw,where v is the bulk flow velocity in the fractures. Critical values of Re can distinguish between the linear and nonlinear flow conditions [20, 22, 24, 31, 46, 47].The vector flow velocityv can be reduced to a one-dimensional flow velocity along a fracture intersection:(7)v=QA.


## 2.1.1. Linear Darcy Flow Zone

The simple parallel plate model is the only fracture model available to calculate the hydraulic conductivity of a fracture. This model assumes a steady incompressible Newtonian fluid flow in a single fracture under a one-dimensional pressure gradient between two smooth parallel plates separated by an aperturee [24, 25]. At sufficiently low Q, this calculation yields the well-known “cubic law” [26, 27] and “Darcy’s law” [28, 29].(1)Q=-e3w12μ∇P=-kAμ∇P,where Q is the total volumetric flow rate, which represents the flow rate of the influent fracture in this study, ∇P is the pressure gradient, w is the width of the fracture, e is the uniform aperture of the idealized smooth fracture, μ is the dynamic viscosity,k is the intrinsic permeability defined as e2/12[5, 27], and A is the cross-sectional area equal to ew. Equation (1) predicts that Q is proportional to the cube of e and that Q is linearly correlated with ∇P when Re is sufficiently small. The cubic law has been widely applied to rough-walled fractures in rock; however, in these applications, e is replaced by the so-called hydraulic aperture (eH) [24, 30, 31]. In this paper, eH is equal to e since the walls of the fractures in the test specimen are smooth.


## 2.1.2. Non-Darcy Flow Zone

Darcy’s law and the cubic law are only valid when the inertial force is negligible, compared with the magnitude of the viscous force, and this condition is anticipated only whenQ is low. Nonlinear flow occurs when, with an incremental increase in Re, Q increases more than a proportional incremental amount [26, 32–35]. Forchheimer’s law [8, 36–38] has been most widely used to describe the nonlinear flow in fractures and porous media, especially in strong inertial regimes:(2)-∇P=a′Q+b′Q2,where a′ and b′ are coefficients, a′=12μ/we3 and b′=βρ/w2e2, ρ is the fluid density, and β is called either the non-Darcy coefficient or the inertial resistance [5, 39, 40].The hydraulic gradientJ is proportional to ∇P:(3)J=h1-h2ΔL=1ρgP1-P2ΔL=1ρg∇P∇L=1ρg∇P,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. Equation (2) can be written in a Forchheimer’s law form as follows:(4)J=aQ+bQ2.The term aQ represents the linear gradient, and the term bQ2 represents the nonlinear gradient.Coefficientsa and b in (4) are commonly written as follows:(5a)a=μρgkA=12μρgwe3,(5b)b=βgA2=βgw2e2.Equation (5a) states that a is linearly proportional to the reciprocal of the cube of e. If the values of w and e are measurable and constant, the value of a can be calculated directly. In addition, in some studies, the value of a is considered the intrinsic permeability of the fractured media [29, 41, 42]. The coefficient b in (5b) is determined from the geometrical properties of the fractures and media in the experiments [34, 40, 43].For fluid flow within fractures, the ratio of the inertial force to the viscous force is defined as the Reynolds number (Re), which is expressed as follows [44, 45]: (6)Re=ρveμ=ρQμw,where v is the bulk flow velocity in the fractures. Critical values of Re can distinguish between the linear and nonlinear flow conditions [20, 22, 24, 31, 46, 47].The vector flow velocityv can be reduced to a one-dimensional flow velocity along a fracture intersection:(7)v=QA.


## 2.2. Experimental Setup

### 2.2.1. Specimen Preparation

First, specimens containing one inlet and two outlets were designed by AutoCAD, which can accurately create a regular single intersection of three fractures. Based on the control variate method, the specimens were designed withθ = 30°, 60°, 90°, 120°, and 150°, e = 2, 3, 5, 10, and 20 mm, andRr = 40, 100, 200, 300, and 400 mm. The specimens were made with PMMA (polymethyl methacrylate) material, which is transparent and easy to process. Then, the specimen design diagrams, in EPS format, were transferred to the computer of the dedicated PMMA machine. Next, the specimens were created at a series of scales with a small error: ±0.02 mm. The sizes of e,Rr, and w were precisely measured by a digital Vernier caliper. Note that the 5.7 mm width w was controlled by the length of the drill employed to cut the PMMA. The error of w was ±0.1 mm (i.e., w = 5.648–5.745 mm) because the surface of the initial PMMA slab was not flat. Next, a high-pressure flame gun was used to polish the fracture walls, so that the initial condition of smooth fracture walls was met in this study. Figure 1 shows the pictures of the specimens. Finally, the specimens were glued to plates with the same side area as the specimens using the PMMA adhesive. This procedure was repeated to ensure absolute sealing of the fractures.Figure 1
The pictures of specimens for the cases ofθ = 30°, 60°, 90°, 120°, and 150°.
(a)
θ = 30° (b)
θ = 60° (c)
θ = 90° (d)
θ = 120° (e)
θ = 150°


### 2.2.2. Test System and Experimental Procedure

A schematic view of the fluid flow test system is shown in Figure2(a). Filtered tap water was introduced into the inlet of the specimen by a high-pressure syringe pump; the syringe pressure had an accuracy of ±0.01 MPa/m. The syringe pressure gradient range of 0.20–1.80 MPa/m resulted in the range of 2.00–9.00 L/min (3.33 × 10−5–1.50 × 10−4 m3/s) for the syringe flow rates. A digital pressure gauge with an accuracy of ±0.1 kPa/m was used to measure the pressure between the inlet and the outlets. The specimens were placed on a horizontal platform fixed by transparent tape so that the gravitational pressure difference between the inlet and the outlets was negligible. Two buckets, each of which had a large hole in the side, were employed to collect the discharged water. Electronic scales with an accuracy of ±0.02 g were used to measure the weight of the effluent. Fluid flow through the transparent PMMA specimen was recorded by a SLR (single lens reflex) camera mounted above the experimental setup. To connect the tube of the syringe pump to the inflow tank of the model, at the end of each fracture, a square slot with a size of 5 × 8 × 8 mm was left to fit to be connected with the pressure-proof PVC tubule of a 6 mm inner diameter and 8 mm outer diameter. The connection between the PVC tubule and the inflow tank of the model was then sealed using a hot glue gun as shown in Figure 2(b) in the revised manuscript. In addition, the other end of the PVC inlet tubule of the inflow tank was quickly connected with the tube of syringe pump via an inflow butt joint. The entire flow path, including the connecting tubes, inflow butt joint, and all three fractures, was sealed and watertight. Filtered tap water was used as the fluid for the experiments, with a density of 998.2 kg/m3 and a dynamic viscosity μ of 0.001 Pa·s at a room temperature of 20°C.Figure 2
(a) The fluid flow test system; (b) the connection between the PVC tubule and the inflow tank of the model.
(a) (b)


## 2.2.1. Specimen Preparation

First, specimens containing one inlet and two outlets were designed by AutoCAD, which can accurately create a regular single intersection of three fractures. Based on the control variate method, the specimens were designed withθ = 30°, 60°, 90°, 120°, and 150°, e = 2, 3, 5, 10, and 20 mm, andRr = 40, 100, 200, 300, and 400 mm. The specimens were made with PMMA (polymethyl methacrylate) material, which is transparent and easy to process. Then, the specimen design diagrams, in EPS format, were transferred to the computer of the dedicated PMMA machine. Next, the specimens were created at a series of scales with a small error: ±0.02 mm. The sizes of e,Rr, and w were precisely measured by a digital Vernier caliper. Note that the 5.7 mm width w was controlled by the length of the drill employed to cut the PMMA. The error of w was ±0.1 mm (i.e., w = 5.648–5.745 mm) because the surface of the initial PMMA slab was not flat. Next, a high-pressure flame gun was used to polish the fracture walls, so that the initial condition of smooth fracture walls was met in this study. Figure 1 shows the pictures of the specimens. Finally, the specimens were glued to plates with the same side area as the specimens using the PMMA adhesive. This procedure was repeated to ensure absolute sealing of the fractures.Figure 1
The pictures of specimens for the cases ofθ = 30°, 60°, 90°, 120°, and 150°.
(a)
θ = 30° (b)
θ = 60° (c)
θ = 90° (d)
θ = 120° (e)
θ = 150°


## 2.2.2. Test System and Experimental Procedure

A schematic view of the fluid flow test system is shown in Figure2(a). Filtered tap water was introduced into the inlet of the specimen by a high-pressure syringe pump; the syringe pressure had an accuracy of ±0.01 MPa/m. The syringe pressure gradient range of 0.20–1.80 MPa/m resulted in the range of 2.00–9.00 L/min (3.33 × 10−5–1.50 × 10−4 m3/s) for the syringe flow rates. A digital pressure gauge with an accuracy of ±0.1 kPa/m was used to measure the pressure between the inlet and the outlets. The specimens were placed on a horizontal platform fixed by transparent tape so that the gravitational pressure difference between the inlet and the outlets was negligible. Two buckets, each of which had a large hole in the side, were employed to collect the discharged water. Electronic scales with an accuracy of ±0.02 g were used to measure the weight of the effluent. Fluid flow through the transparent PMMA specimen was recorded by a SLR (single lens reflex) camera mounted above the experimental setup. To connect the tube of the syringe pump to the inflow tank of the model, at the end of each fracture, a square slot with a size of 5 × 8 × 8 mm was left to fit to be connected with the pressure-proof PVC tubule of a 6 mm inner diameter and 8 mm outer diameter. The connection between the PVC tubule and the inflow tank of the model was then sealed using a hot glue gun as shown in Figure 2(b) in the revised manuscript. In addition, the other end of the PVC inlet tubule of the inflow tank was quickly connected with the tube of syringe pump via an inflow butt joint. The entire flow path, including the connecting tubes, inflow butt joint, and all three fractures, was sealed and watertight. Filtered tap water was used as the fluid for the experiments, with a density of 998.2 kg/m3 and a dynamic viscosity μ of 0.001 Pa·s at a room temperature of 20°C.Figure 2
(a) The fluid flow test system; (b) the connection between the PVC tubule and the inflow tank of the model.
(a) (b)


## 3. Results and Discussions

As shown in Figure3, θ is defined as the intersection angle between the effluent fractures in a 2D plane and can vary from 0° to 180°. Pictures of the specimens for the cases of θ = 30°, 60°, 90°, 120°, and 150° are presented in Figure 1. The effluent fracture located in the orientation of influent fracture is called fracture_1, while the effluent fracture located in the deflected orientation is called fracture_2. Moreover, in this paper, in the cases of θ = 30°, 60°, 90°, 120°, and 150°, fracture_1 is called fracture_11, fracture_12, fracture_13, fracture_14, and fracture_15, respectively. Similarly, in the cases of θ = 30°, 60°, 90°, 120°, and 150°, fracture_2 is called fracture_21, fracture_22, fracture_23, fracture_24, and fracture_25.Figure 3
Definition ofθ and fractures.


### 3.1. The Nonlinear Behaviour of Fluid Flow in Intersecting Fractures

The tests in Figure4 were conducted on specimens with e = 2 mm,Rr = 40 mm, w = 5.7 mm, and θ = 30°, 60°, 90°, 120°, and 150°. The specimens were full of fluid and flow was continuous in the influent fracture as well as both effluent fractures. The fitting result of the test data in Figure 4(a) shows that J/Q is linearly correlated with Q, which can be expressed as J/Q=a+bQ. The relationship between J and Q agrees with Forchheimer’s law. In the tests, as J varies from 1.184 to 31.270, Q varies from 3.296×10-5 to 1.354×10-4 m3/s and Re of the influent fracture varies from 5782.456 to 23759.649. Researchers [5, 20, 22, 23, 48] have concluded that Forchheimer’s law applies for the cases in which the magnitude of Q is 10−7–10−6 and the magnitude of Re is 101–102. In addition, in this paper, it is demonstrated that Forchheimer’s law also applies for larger magnitudes of Q (10-5–10-4) and Re (104–105).Figure 4
Relationship betweenJ/Q and Q as well as relationship between J and Q; (a) relationship between J/Q and Q; (b) relationship between J and Q for θ = 30°; (c) relationship between J and Q for θ = 60°; (d) relationship between J and Q for θ = 90°; (e) relationship between J and Q for θ = 120°; (f) relationship between J and Q for θ = 150°.
(a) (b) (c) (d) (e) (f)The variations ofq1, q2, and Q as well as the cubic law-based Q with increasing J for the cases of θ = 30°, 60°, 90°, 120°, and 150° are shown in Figures 4(b)–4(f). J is quadratically correlated with q1, q2, and Q with a correlation coefficient R2 > 0.955. WhenJ = 100–101, as shown in Figure 4 the J-q1, J-q2, and J-Q curves deviate from straight lines (cubic law) and the deviations increase with increasing J. These results agree with the experimental and numerical results of Koyama et al. [49] and Li et al. [23, 50]. As shown in Figures 4(b)–4(f), the deviations of the J-q2, J-q1, and J-Q curves from cubic law line are all in a decreasing order. In this group of tests, e and w are considered to be constant in (5a) and (5b); therefore, the coefficient a of the linear term in (4) is calculated as the constant value 2.1465 × 104. Table 1 provides the coefficient b of the nonlinear term in (4) for the cases of θ = 30°, 60°, 90°, 120°, and 150°. The coefficient b can quantify the deviations of the J-q1, J-q2, and J-Q curves. The coefficient b of fracture_1 decreases with increasing θ, especially when θ ≥ 90° (i.e., as θ increases from 30° to 90°, b decreases from 2.9477 × 109 to 2.2524 × 109). In contrast, the coefficient b of fracture_2 increases with increasing θ, especially when θ ≥ 90° (i.e., as θ increases from 30° to 90°, b increases from 2.1500 × 1010 to 5.0517 × 1010). However, the coefficient b of the influent fracture only slightly changes with increasing θ (i.e., as θ increases from 30° to 150°, b varies from 1.4966 × 109 to 1.4981 × 109). This emphasizes that the magnitude of the coefficient b of fracture_1 (or fracture_2) remains on the order of 109 (or 1010) as θ changes.g, w, and eH are measurable and assumed to be constant in (5a), so that the coefficient b is determined by the inertial resistance β. Ruth and Ma [40] and Cherubini et al. [34] suggested that if the structure of a porous medium is such that microscopic inertial effects are rare (i.e., essentially straight fractures), then β will be small unless v and Re are high. Furthermore, for high values of Re, Re and turbulence are closely connected. The direction of the fluid flow within fracture_2 sharply changes because orientation of fracture_2 diverges from that of the influent fracture. Thus, the inertial resistance β of fracture_2 has a significant influence on the inertial force of the fluid flow. Additionally, the inertial resistance β of fracture_1 also affects the inertial force of the fluid flow due to the expansion of the fracture cross-sectional area at the intersection. The inertial resistance β of fracture_2 is always greater than that of fracture_1. In addition, the inertial resistance β of both fracture_2 and fracture_1 exceeds that of the influent fracture. As θ increases from 0° to 180°, the inertial resistance β of fracture_2 intensifies because fracture_2 turns away from direction of the influent fracture and thus the direction of fluid flow. Due to the limit of the syringe pump used in the tests, the influent flow rate is fixed for an equivalent J with different specimens. Hence, for test cases with different θ, the J/Q-Q curves overlap each other, as demonstrated in Figure 4(a). The fixed Q provides a constant inertial force of the fluid flow at the inlet. Therefore, the inertial resistance β of fracture_2 and fracture_1 has different effects on the inertial force with increasing J. β of fracture_2 increases with increasingJ, whereas β of fracture_1 decreases. As J varies, the change in β is clearly similar to that of b.Table 1
The coefficientB of J-q1, q2, and Q for the cases of θ=30°, 60°, 90°, 120°, and 150°.
Coefficientb θ = 30 ° θ = 60 ° θ = 90 ° θ = 120 ° θ = 150 ° J-q1 (×109) 2.9477 2.5986 2.2524 2.4162 2.2645 J-q2 (×1010) 2.1500 2.9410 5.0517 5.2248 5.2181 J-Q (×109) 1.4966 1.4785 1.4816 1.5864 1.4981


### 3.2. Effects ofθon the Nonlinearity of Fluid Flow

When fluid flow is in a linear regime, in this group of tests, the predicted value ofQ/J based on (1) is a constant value of 4.637 × 10−5. On the other hand, in a nonlinear flow regime, the predicted value of Q/J is based on (4) and declines with increasing J. Furthermore, the change in Q/J represents the degree of nonlinearity. Therefore, as J increases, Q/J changes from a constant value to decreasing values, which indicates the transition of fluid flow from the linear regime to the nonlinear regime. Liu et al. [22], using numerical simulation, presented that the critical value of J is between 10−3 and 10−2 for the case ofe = 2 mm. In this paper, the fluid flow within the fractures is certainly in the nonlinear regime, sinceJ = 101–102. The relative change in Q/J, δ, is defined as follows:(8)δ=Q/Jcubic-Q/JexperimentQ/Jcubic×100%,where Q/Jcubic is the value of Q/J calculated by solving the cubic law, Q/Jexperiment is the value of Q/J calculated by analysing data from the experiments, and δ expresses how the true value of Q/J (Q/Jexperiment) deviates from the prediction made by the cubic law (Q/Jcubic). Similarly, δ1 and δ2, the relative discrepancies of fracture_1 and fracture_2, respectively, are calculated as follows:(9a)δ1=q1/Jcubic-q1/Jexperimentq1/Jcubic×100%,(9b)δ2=q2/Jcubic-q2/Jexperimentq2/Jcubic×100%,where q1 and q2 are the volumetric flow rates within fracture_1 and fracture_2, respectively. In this fluid flow test system,(10)Q=q1+q2.By substituting (10) into (8), (9a), and (9b), the relationship among δ1, δ2, and δ can be expressed as follows:(11)δ+1=δ1+δ2.Equation (11) suggests that if δ remains a constant value, δ1 and δ2 are negatively correlated. In other words, if the fluid flow in the influent fracture is nonlinear, the decrease in the nonlinearity of fluid flow in one of the two effluent fractures will be offset by an increase in the other.The tests in Figure5 were conducted on the specimens with e = 2 mm,Rr = 40 mm, w = 5.7 mm, and θ = 30°, 60°, 90°, 120°, and 150°. There was a full flow in the influent fracture and both effluent fractures. As J increases, the variations in Q/J and δ for the cases of θ = 30°, 60°, 90°, 120°, and 150° are indicated in Figure 5. As provided in Table 2, Q/J and δ vary within only a small range as θ varies over a large range (0° < θ < 180°) for the cases ofJ = 7, 12, 17, and 22. ForJ = 7, with θ increasing from 30° to 150°, Q/J of the influent fracture varies from 9.449 × 10−6 to 9.267 × 10−6 and δ varies from 79.719% to 80.108%. The nearly stable values of Q/J and δ indicate the nearly invariable nonlinearity of the fluid flow within the influent fracture for the cases of θ = 30°, 60°, 90°, 120°, and 150°. As shown in Figures 5(a) and 5(b), for the cases of 7 ≤J ≤ 22, q2/J decreases and δ2 increases with increasing θ, especially when θ < 90°; as θ increases, q1/J and δ1 have opposite changes in magnitude compared with those of q2/J and δ2. q1/J and q2/J are negatively and symmetrically correlated, similar to δ1 and δ2. Thus, increasing θ enhances the nonlinearity of the fluid flow in fracture_2 and decreases that of fracture_1. Note that δ2 > δ1 > δ. Clearly, the nonlinearity of fluid flow in fracture_2 is greatest while that of the influent fracture is the least of the three fractures. As θ increases, forJ = 7, δ1 decreases from 85.300% to 83.235% and δ2 increases from 94.518% to 96.336%; forJ = 22, δ1 decreases from 91.993% to 90.913% and δ2 increases from 96.979% to 97.967%. In fact, q1/J, q2/J, δ1, and δ2 only slightly vary with varying θ. The fluid flow nonlinearities in fracture_2 and fracture_1 are considerably greater because δ2 and δ1 are both approaching 100% for the cases ofJ = 101–102. Moreover, the changes in δ2 and δ1 diminish when θ > 90° (i.e., whenJ = 7, for θ = 90°, δ1 = 83.335% and δ2 = 96.368%; for θ = 120°, δ1 = 83.638% and δ2 = 96.368%; and for θ = 150°, δ1 = 83.235% and δ2 = 96.336%). The inertial resistance caused by the fracture intersection changes considerably before the direction of fluid flow in fracture_2 is at a right angle to that in the influent fracture; when the fluid flow direction changes more (i.e., fluid flow turns back), the inertial resistance is constant.Table 2
Experimental results ofQ/J and δ in influent fracture for the cases of θ = 30°, 60°, 90°, and 120° with J = 7, 12, 17, and 22.
θ Q / J (×10−6 m3/s) δ (%) J = 7 J = 12 J = 17 J = 22 J = 7 J = 12 J = 17 J = 22 30° 9.449 7.065 5.999 5.177 79.719 84.837 87.125 88.889 60° 9.424 7.117 5.973 5.221 79.771 84.723 87.180 88.794 90° 9.454 7.054 5.978 5.203 79.708 84.858 87.169 88.833 120° 9.342 7.181 5.879 5.193 79.947 84.587 87.381 88.853 150° 9.267 6.979 5.943 5.187 80.108 85.019 87.243 88.866Figure 5
(a) Relationship betweenQ/J and θ for J = 7, 12, 17, and 22; (b) relationship between δ and θ for J = 7, 12, 17, and 22.
(a) (b)


### 3.3. Effects ofJ on the Nonlinearity of Fluid Flow

The tests in Figure6 were conducted on the specimens with e = 2 mm,Rr = 40 mm, w = 5.7 mm, and θ = 30°, 60°, 90°, 120°, and 150°. The fluid flow among the influent fractures and effluent fractures filled the fractures. In this group of tests, the predicted value of Q/J, based on (1), is a constant value of 4.637 × 10−5. With the changes of J, Q/J and δ for the cases of θ = 30°, 60°, 90°, 120°, and 150° are shown in Figure 6. Figures 6(a) and 6(b) show the trendlines of Q/J-J and δ-J for the cases ofJ = 10−5–100 from Li et al. [23]. For the cases ofJ = 2–32, the changes in the Q/J-J and δ-J curves are exactly equal to the trendlines of Li et al., and the trendlines of Q/J-J and δ-J are extended to the cases ofJ = 100–101. As shown in Figures 6(a) and 6(b), when 2 ≤ J ≤ 35 and J increases, Q/J, q1/J, and q2/J decrease whereas δ2, δ1, and δ increase. Therefore, increasing J strengthens the nonlinearity of fluid flow in all fractures. An increase in J would increase the number and volume of eddies in the fractures, narrowing the effective flow paths, and consequently enhance the nonlinearity of fluid flow [22]. A large magnitude of the change in the nonlinearity arises for the especial cases of 100 <J < 101. However, the magnitude of the change decreases whenJ > 101. For example, for the case of θ = 120°, when 100 <J < 101, δ increases from 65.959 to 81.557% as J increases from 2.736 to 8.167; whenJ > 101, δ increases from 85.872 to 91.061% as J increases from 13.512 to 34.201. As indicated in the enlarged views of Figures 6(a) and 6(b), as J increases, q1/J increases (or q2/J decreases) and δ1 decreases (or δ2 increases). Note that δ2 > δ1 > δ. According to 3.2, when θ ≥ 90°, varying θ has a negligible effect on the nonlinearity of fluid flow in fracture_1 and fracture_2. Thus, the Q/J-J or δ-J curves of fracture_1 and fracture_2 overlap for θ = 90°, 120°, and 150°. Furthermore, varying θ has a negligible effect on the nonlinearity of fluid flow in the influent fracture, resulting in the overlapping Q/J-J or δ-J influent fracture curves with different θ.Figure 6
(a) Relationship betweenQ/J and J for θ = 30°, 60°, 90°, 120°, and 150°; (b) relationship between δ and J for θ = 30°, 60°, 90°, 120°, and 150°.
(a) (b)


### 3.4. Effects ofe on the Nonlinearity of Fluid Flow

The tests in Figure7 were conducted on the specimens withe = 2, 3, 5, 10, and 20 mm,Rr = 100 mm, w = 5.7 mm, and θ = 60°, 90°, and 120°. There was a full flow within all fractures for the cases withe = 2, 3, and 5 mm. However, for the cases ofe = 10 and 20 mm, the fluid flow in fracture_2 did not completely fill the open fracture space. For the case ofe = 10 mm, in fracture_2, there were many bubbles that had diameters greater than 1 mm. For the case ofJ ≈ 10, as θ increases from 60° to 120°, q2 decreased from 2.831 × 10−5 to 1.930 × 10−5 m3/s; v of fracture_2 decreased to 0.412 m/s (for θ = 60°), 0.301 m/s (for θ = 90°), and 0.189 m/s (for θ = 120°), although v of the influent fracture was 1.442 m/s (for θ = 60°), 1.409 m/s (for θ = 90°), and 1.424 m/s (for θ = 120°). For the case of e = 20 mm, when θ = 60° and 90°, there were more bubbles with diameters greater than 15 mm in fracture_2 (see Figures 8(a) and 8(b)); when θ = 120°, there was only a small volume of fluid flow through fracture_2 (see Figure 8(c)). For the case of e = 20 mm, v of fracture_2 decreased to 0.141 m/s (for θ = 60°), 0.054 m/s (for θ = 90°), and 0.001 m/s (for θ = 120°), although v of influent fracture was 0.723 m/s (for θ = 60°), 0.660 m/s (for θ = 90°), and 0.729 m/s (for θ = 120°). Hence, v of fracture_2 sharply decreases with increasing θ, especially at e = 20 mm. Clearly, each bubble is caused by an eddy that occurs because of the sudden decrease in velocity of the fluid flow and the retroflex direction of the flow. As shown in the experiment diagrams (see Figure 8), there was a curved surface boundary encompassing one or more eddies (i.e., eddy_1, eddy_2, and eddy_3) at the fracture intersection, which reduced to a curvilinear boundary in the 2D diagram. The curvilinear boundary is the separatrix of the fluid flow, which means that only a little fluid flow can cross it. In addition, as θ increases from 60° to 120°, the curvilinear boundary moves from fracture_2 to fracture_1. Therefore, with a change in θ, the flow decreases, crosses the curvilinear boundary, and moves into fracture_2, intensifying the non-full-flow state within fracture_2. Therefore, eddy_2 is larger than eddy_1 and the eddy diminishes in fracture_2 with a small fluid flow for the cases ofe = 20 mm and θ = 120°. In fact, fore = 20 mm and θ = 120°, a small volume of fluid turns back into fracture_2 due to eddy_3 (see Figure 8(c)).Figure 7
(a) Relationship betweenQ/J and e for θ = 60°, 90°, and 120°; (b) relationship between δ and e for θ = 60°, 90°, and 120°.
(a) (b)Figure 8
The experimental phenomena for the cases of (a)e = 20 mm,J = 10.701, and θ = 60°; (b) e = 20 mm,J = 10.529, and θ = 90°; and (c) e = 20 mm,J = 10.046, and θ = 120°.
(a) (b) (c)In this paper, the relationship betweenJ and q2 agrees with the Forchheimer’s law, though fluid flow in fracture_2 is in a non-full-flow state. With changinge, Q/J and δ for the cases of θ = 60°, 90°, and 120° are shown in Figure 7. As shown in Figure 7(b), δ2, δ1, and δ increase with increasinge, whereas Q/J, q1/J, and q2/J vary in different ways (i.e., as e increases, Q/J remains constant, q1/J increases, and q2/J decreases). In fact, q1/J and q2/J have a negative and symmetric correlation. From the enlarged view given in Figure 7(b), δ2 is still greater than δ1, which is still greater than δ. Based on (5a), the coefficient a is linearly proportional to the reciprocal of the cube of e. Thus, the coefficient a rapidly decreases with widening e, causing the broad changes in the cubic law-based Q/J. Considering that the basic value of the cubic law-based Q/J is not constant, in this group of tests, the changes in Q/J, q1/J, and q2/J in Figure 7(a) cannot describe the degree of nonlinearity of the fluid flow. δ2, δ1, and δ approach 100% for e > 5 mm. For example, when θ = 60°, as e widens from 5 to 20 mm, δ2 increases from 99.668 to 99.996%, δ1 increases from 99.008 to 99.982%, and δ increases from 98.627 to 99.977%. Note that varying e has a negligible effect on Q/J because the syringe pump in the tests can provide a stable flow rate for each hydraulic gradient. Unfortunately, the differences in the q1/J-e, q2/J-e, δ1-e, and δ2-e curves caused by the varying θ have not yet been presented and require more attention in future experimental studies.


### 3.5. Effects ofRr on the Nonlinearity of Fluid Flow

The tests in Figure9 were conducted on the specimens with e = 2 mm,Rr = 40, 100, 200, 300, and 400 mm, w = 5.7 mm, and θ = 60°, 90°, and 120°. There was a full flow in the influent fracture as well as both effluent fractures. In this group of tests, the predicted value of Q/J based on (1) is still the constant value 4.637 × 10−5. AsRr increases, the changes in Q/J and δ for all three fractures for the cases of θ = 60°, 90°, and 120° are shown in Figure 9. As shown in Figures 9(a) and 9(b), whenJ = 10, q2/J increases and δ2 declines with increasingRr, especially forRr < 100 mm. However, the magnitudes of q1/J and δ1 have different responses asRr increases; as one parameter increases, the other decreases. In fact, q1/J and q2/J have a negative and symmetric correlation, similar to δ1 and δ2. Thus, increasingRr enhances the nonlinearity of fluid flow in fracture_1 but decreases that in fracture_2. The consequence is that the magnitude of these changes in q2/J and δ2 diminishes whenRr > 100 mm, similar to the result obtained by Li et al. [23]. Therefore, if the fracture length is long, the disturbance on the fluid flow resulting from the intersection would become negligibly small. Furthermore, the on-way resistance resulting from the long effluent fractures decreases the dominant effect of the intersection. AsRr increases, the inertial resistances of fracture_1 and fracture_2 caused by the fracture length and the intersection decrease. In addition, as the on the way resistance becomes more dominant, the nonlinearity of fluid flow approaches a constant value. Note that the nonlinearity of fluid flow in fracture_2 is still greater than fracture_1, which is still greater than that in the influent fracture. As indicated in the enlarged views of Figures 9(a) and 9(b), as θ increases, q1/J increases (or q2/J decreases) and δ1 decreases (or δ2 increases). In addition, varyingRr and θ have a negligible effect on the nonlinearity of fluid flow in the influent fracture, resulting in a constant value of δ.Figure 9
(a) Relationship betweenQ/J andRr for θ = 60°, 90°, and 120°; (b) relationship between δ andRr for θ = 60°, 90°, and 120°.
(a) (b)


## 3.1. The Nonlinear Behaviour of Fluid Flow in Intersecting Fractures

The tests in Figure4 were conducted on specimens with e = 2 mm,Rr = 40 mm, w = 5.7 mm, and θ = 30°, 60°, 90°, 120°, and 150°. The specimens were full of fluid and flow was continuous in the influent fracture as well as both effluent fractures. The fitting result of the test data in Figure 4(a) shows that J/Q is linearly correlated with Q, which can be expressed as J/Q=a+bQ. The relationship between J and Q agrees with Forchheimer’s law. In the tests, as J varies from 1.184 to 31.270, Q varies from 3.296×10-5 to 1.354×10-4 m3/s and Re of the influent fracture varies from 5782.456 to 23759.649. Researchers [5, 20, 22, 23, 48] have concluded that Forchheimer’s law applies for the cases in which the magnitude of Q is 10−7–10−6 and the magnitude of Re is 101–102. In addition, in this paper, it is demonstrated that Forchheimer’s law also applies for larger magnitudes of Q (10-5–10-4) and Re (104–105).Figure 4
Relationship betweenJ/Q and Q as well as relationship between J and Q; (a) relationship between J/Q and Q; (b) relationship between J and Q for θ = 30°; (c) relationship between J and Q for θ = 60°; (d) relationship between J and Q for θ = 90°; (e) relationship between J and Q for θ = 120°; (f) relationship between J and Q for θ = 150°.
(a) (b) (c) (d) (e) (f)The variations ofq1, q2, and Q as well as the cubic law-based Q with increasing J for the cases of θ = 30°, 60°, 90°, 120°, and 150° are shown in Figures 4(b)–4(f). J is quadratically correlated with q1, q2, and Q with a correlation coefficient R2 > 0.955. WhenJ = 100–101, as shown in Figure 4 the J-q1, J-q2, and J-Q curves deviate from straight lines (cubic law) and the deviations increase with increasing J. These results agree with the experimental and numerical results of Koyama et al. [49] and Li et al. [23, 50]. As shown in Figures 4(b)–4(f), the deviations of the J-q2, J-q1, and J-Q curves from cubic law line are all in a decreasing order. In this group of tests, e and w are considered to be constant in (5a) and (5b); therefore, the coefficient a of the linear term in (4) is calculated as the constant value 2.1465 × 104. Table 1 provides the coefficient b of the nonlinear term in (4) for the cases of θ = 30°, 60°, 90°, 120°, and 150°. The coefficient b can quantify the deviations of the J-q1, J-q2, and J-Q curves. The coefficient b of fracture_1 decreases with increasing θ, especially when θ ≥ 90° (i.e., as θ increases from 30° to 90°, b decreases from 2.9477 × 109 to 2.2524 × 109). In contrast, the coefficient b of fracture_2 increases with increasing θ, especially when θ ≥ 90° (i.e., as θ increases from 30° to 90°, b increases from 2.1500 × 1010 to 5.0517 × 1010). However, the coefficient b of the influent fracture only slightly changes with increasing θ (i.e., as θ increases from 30° to 150°, b varies from 1.4966 × 109 to 1.4981 × 109). This emphasizes that the magnitude of the coefficient b of fracture_1 (or fracture_2) remains on the order of 109 (or 1010) as θ changes.g, w, and eH are measurable and assumed to be constant in (5a), so that the coefficient b is determined by the inertial resistance β. Ruth and Ma [40] and Cherubini et al. [34] suggested that if the structure of a porous medium is such that microscopic inertial effects are rare (i.e., essentially straight fractures), then β will be small unless v and Re are high. Furthermore, for high values of Re, Re and turbulence are closely connected. The direction of the fluid flow within fracture_2 sharply changes because orientation of fracture_2 diverges from that of the influent fracture. Thus, the inertial resistance β of fracture_2 has a significant influence on the inertial force of the fluid flow. Additionally, the inertial resistance β of fracture_1 also affects the inertial force of the fluid flow due to the expansion of the fracture cross-sectional area at the intersection. The inertial resistance β of fracture_2 is always greater than that of fracture_1. In addition, the inertial resistance β of both fracture_2 and fracture_1 exceeds that of the influent fracture. As θ increases from 0° to 180°, the inertial resistance β of fracture_2 intensifies because fracture_2 turns away from direction of the influent fracture and thus the direction of fluid flow. Due to the limit of the syringe pump used in the tests, the influent flow rate is fixed for an equivalent J with different specimens. Hence, for test cases with different θ, the J/Q-Q curves overlap each other, as demonstrated in Figure 4(a). The fixed Q provides a constant inertial force of the fluid flow at the inlet. Therefore, the inertial resistance β of fracture_2 and fracture_1 has different effects on the inertial force with increasing J. β of fracture_2 increases with increasingJ, whereas β of fracture_1 decreases. As J varies, the change in β is clearly similar to that of b.Table 1
The coefficientB of J-q1, q2, and Q for the cases of θ=30°, 60°, 90°, 120°, and 150°.
Coefficientb θ = 30 ° θ = 60 ° θ = 90 ° θ = 120 ° θ = 150 ° J-q1 (×109) 2.9477 2.5986 2.2524 2.4162 2.2645 J-q2 (×1010) 2.1500 2.9410 5.0517 5.2248 5.2181 J-Q (×109) 1.4966 1.4785 1.4816 1.5864 1.4981


## 3.2. Effects ofθon the Nonlinearity of Fluid Flow

When fluid flow is in a linear regime, in this group of tests, the predicted value ofQ/J based on (1) is a constant value of 4.637 × 10−5. On the other hand, in a nonlinear flow regime, the predicted value of Q/J is based on (4) and declines with increasing J. Furthermore, the change in Q/J represents the degree of nonlinearity. Therefore, as J increases, Q/J changes from a constant value to decreasing values, which indicates the transition of fluid flow from the linear regime to the nonlinear regime. Liu et al. [22], using numerical simulation, presented that the critical value of J is between 10−3 and 10−2 for the case ofe = 2 mm. In this paper, the fluid flow within the fractures is certainly in the nonlinear regime, sinceJ = 101–102. The relative change in Q/J, δ, is defined as follows:(8)δ=Q/Jcubic-Q/JexperimentQ/Jcubic×100%,where Q/Jcubic is the value of Q/J calculated by solving the cubic law, Q/Jexperiment is the value of Q/J calculated by analysing data from the experiments, and δ expresses how the true value of Q/J (Q/Jexperiment) deviates from the prediction made by the cubic law (Q/Jcubic). Similarly, δ1 and δ2, the relative discrepancies of fracture_1 and fracture_2, respectively, are calculated as follows:(9a)δ1=q1/Jcubic-q1/Jexperimentq1/Jcubic×100%,(9b)δ2=q2/Jcubic-q2/Jexperimentq2/Jcubic×100%,where q1 and q2 are the volumetric flow rates within fracture_1 and fracture_2, respectively. In this fluid flow test system,(10)Q=q1+q2.By substituting (10) into (8), (9a), and (9b), the relationship among δ1, δ2, and δ can be expressed as follows:(11)δ+1=δ1+δ2.Equation (11) suggests that if δ remains a constant value, δ1 and δ2 are negatively correlated. In other words, if the fluid flow in the influent fracture is nonlinear, the decrease in the nonlinearity of fluid flow in one of the two effluent fractures will be offset by an increase in the other.The tests in Figure5 were conducted on the specimens with e = 2 mm,Rr = 40 mm, w = 5.7 mm, and θ = 30°, 60°, 90°, 120°, and 150°. There was a full flow in the influent fracture and both effluent fractures. As J increases, the variations in Q/J and δ for the cases of θ = 30°, 60°, 90°, 120°, and 150° are indicated in Figure 5. As provided in Table 2, Q/J and δ vary within only a small range as θ varies over a large range (0° < θ < 180°) for the cases ofJ = 7, 12, 17, and 22. ForJ = 7, with θ increasing from 30° to 150°, Q/J of the influent fracture varies from 9.449 × 10−6 to 9.267 × 10−6 and δ varies from 79.719% to 80.108%. The nearly stable values of Q/J and δ indicate the nearly invariable nonlinearity of the fluid flow within the influent fracture for the cases of θ = 30°, 60°, 90°, 120°, and 150°. As shown in Figures 5(a) and 5(b), for the cases of 7 ≤J ≤ 22, q2/J decreases and δ2 increases with increasing θ, especially when θ < 90°; as θ increases, q1/J and δ1 have opposite changes in magnitude compared with those of q2/J and δ2. q1/J and q2/J are negatively and symmetrically correlated, similar to δ1 and δ2. Thus, increasing θ enhances the nonlinearity of the fluid flow in fracture_2 and decreases that of fracture_1. Note that δ2 > δ1 > δ. Clearly, the nonlinearity of fluid flow in fracture_2 is greatest while that of the influent fracture is the least of the three fractures. As θ increases, forJ = 7, δ1 decreases from 85.300% to 83.235% and δ2 increases from 94.518% to 96.336%; forJ = 22, δ1 decreases from 91.993% to 90.913% and δ2 increases from 96.979% to 97.967%. In fact, q1/J, q2/J, δ1, and δ2 only slightly vary with varying θ. The fluid flow nonlinearities in fracture_2 and fracture_1 are considerably greater because δ2 and δ1 are both approaching 100% for the cases ofJ = 101–102. Moreover, the changes in δ2 and δ1 diminish when θ > 90° (i.e., whenJ = 7, for θ = 90°, δ1 = 83.335% and δ2 = 96.368%; for θ = 120°, δ1 = 83.638% and δ2 = 96.368%; and for θ = 150°, δ1 = 83.235% and δ2 = 96.336%). The inertial resistance caused by the fracture intersection changes considerably before the direction of fluid flow in fracture_2 is at a right angle to that in the influent fracture; when the fluid flow direction changes more (i.e., fluid flow turns back), the inertial resistance is constant.Table 2
Experimental results ofQ/J and δ in influent fracture for the cases of θ = 30°, 60°, 90°, and 120° with J = 7, 12, 17, and 22.
θ Q / J (×10−6 m3/s) δ (%) J = 7 J = 12 J = 17 J = 22 J = 7 J = 12 J = 17 J = 22 30° 9.449 7.065 5.999 5.177 79.719 84.837 87.125 88.889 60° 9.424 7.117 5.973 5.221 79.771 84.723 87.180 88.794 90° 9.454 7.054 5.978 5.203 79.708 84.858 87.169 88.833 120° 9.342 7.181 5.879 5.193 79.947 84.587 87.381 88.853 150° 9.267 6.979 5.943 5.187 80.108 85.019 87.243 88.866Figure 5
(a) Relationship betweenQ/J and θ for J = 7, 12, 17, and 22; (b) relationship between δ and θ for J = 7, 12, 17, and 22.
(a) (b)


## 3.3. Effects ofJ on the Nonlinearity of Fluid Flow

The tests in Figure6 were conducted on the specimens with e = 2 mm,Rr = 40 mm, w = 5.7 mm, and θ = 30°, 60°, 90°, 120°, and 150°. The fluid flow among the influent fractures and effluent fractures filled the fractures. In this group of tests, the predicted value of Q/J, based on (1), is a constant value of 4.637 × 10−5. With the changes of J, Q/J and δ for the cases of θ = 30°, 60°, 90°, 120°, and 150° are shown in Figure 6. Figures 6(a) and 6(b) show the trendlines of Q/J-J and δ-J for the cases ofJ = 10−5–100 from Li et al. [23]. For the cases ofJ = 2–32, the changes in the Q/J-J and δ-J curves are exactly equal to the trendlines of Li et al., and the trendlines of Q/J-J and δ-J are extended to the cases ofJ = 100–101. As shown in Figures 6(a) and 6(b), when 2 ≤ J ≤ 35 and J increases, Q/J, q1/J, and q2/J decrease whereas δ2, δ1, and δ increase. Therefore, increasing J strengthens the nonlinearity of fluid flow in all fractures. An increase in J would increase the number and volume of eddies in the fractures, narrowing the effective flow paths, and consequently enhance the nonlinearity of fluid flow [22]. A large magnitude of the change in the nonlinearity arises for the especial cases of 100 <J < 101. However, the magnitude of the change decreases whenJ > 101. For example, for the case of θ = 120°, when 100 <J < 101, δ increases from 65.959 to 81.557% as J increases from 2.736 to 8.167; whenJ > 101, δ increases from 85.872 to 91.061% as J increases from 13.512 to 34.201. As indicated in the enlarged views of Figures 6(a) and 6(b), as J increases, q1/J increases (or q2/J decreases) and δ1 decreases (or δ2 increases). Note that δ2 > δ1 > δ. According to 3.2, when θ ≥ 90°, varying θ has a negligible effect on the nonlinearity of fluid flow in fracture_1 and fracture_2. Thus, the Q/J-J or δ-J curves of fracture_1 and fracture_2 overlap for θ = 90°, 120°, and 150°. Furthermore, varying θ has a negligible effect on the nonlinearity of fluid flow in the influent fracture, resulting in the overlapping Q/J-J or δ-J influent fracture curves with different θ.Figure 6
(a) Relationship betweenQ/J and J for θ = 30°, 60°, 90°, 120°, and 150°; (b) relationship between δ and J for θ = 30°, 60°, 90°, 120°, and 150°.
(a) (b)


## 3.4. Effects ofe on the Nonlinearity of Fluid Flow

The tests in Figure7 were conducted on the specimens withe = 2, 3, 5, 10, and 20 mm,Rr = 100 mm, w = 5.7 mm, and θ = 60°, 90°, and 120°. There was a full flow within all fractures for the cases withe = 2, 3, and 5 mm. However, for the cases ofe = 10 and 20 mm, the fluid flow in fracture_2 did not completely fill the open fracture space. For the case ofe = 10 mm, in fracture_2, there were many bubbles that had diameters greater than 1 mm. For the case ofJ ≈ 10, as θ increases from 60° to 120°, q2 decreased from 2.831 × 10−5 to 1.930 × 10−5 m3/s; v of fracture_2 decreased to 0.412 m/s (for θ = 60°), 0.301 m/s (for θ = 90°), and 0.189 m/s (for θ = 120°), although v of the influent fracture was 1.442 m/s (for θ = 60°), 1.409 m/s (for θ = 90°), and 1.424 m/s (for θ = 120°). For the case of e = 20 mm, when θ = 60° and 90°, there were more bubbles with diameters greater than 15 mm in fracture_2 (see Figures 8(a) and 8(b)); when θ = 120°, there was only a small volume of fluid flow through fracture_2 (see Figure 8(c)). For the case of e = 20 mm, v of fracture_2 decreased to 0.141 m/s (for θ = 60°), 0.054 m/s (for θ = 90°), and 0.001 m/s (for θ = 120°), although v of influent fracture was 0.723 m/s (for θ = 60°), 0.660 m/s (for θ = 90°), and 0.729 m/s (for θ = 120°). Hence, v of fracture_2 sharply decreases with increasing θ, especially at e = 20 mm. Clearly, each bubble is caused by an eddy that occurs because of the sudden decrease in velocity of the fluid flow and the retroflex direction of the flow. As shown in the experiment diagrams (see Figure 8), there was a curved surface boundary encompassing one or more eddies (i.e., eddy_1, eddy_2, and eddy_3) at the fracture intersection, which reduced to a curvilinear boundary in the 2D diagram. The curvilinear boundary is the separatrix of the fluid flow, which means that only a little fluid flow can cross it. In addition, as θ increases from 60° to 120°, the curvilinear boundary moves from fracture_2 to fracture_1. Therefore, with a change in θ, the flow decreases, crosses the curvilinear boundary, and moves into fracture_2, intensifying the non-full-flow state within fracture_2. Therefore, eddy_2 is larger than eddy_1 and the eddy diminishes in fracture_2 with a small fluid flow for the cases ofe = 20 mm and θ = 120°. In fact, fore = 20 mm and θ = 120°, a small volume of fluid turns back into fracture_2 due to eddy_3 (see Figure 8(c)).Figure 7
(a) Relationship betweenQ/J and e for θ = 60°, 90°, and 120°; (b) relationship between δ and e for θ = 60°, 90°, and 120°.
(a) (b)Figure 8
The experimental phenomena for the cases of (a)e = 20 mm,J = 10.701, and θ = 60°; (b) e = 20 mm,J = 10.529, and θ = 90°; and (c) e = 20 mm,J = 10.046, and θ = 120°.
(a) (b) (c)In this paper, the relationship betweenJ and q2 agrees with the Forchheimer’s law, though fluid flow in fracture_2 is in a non-full-flow state. With changinge, Q/J and δ for the cases of θ = 60°, 90°, and 120° are shown in Figure 7. As shown in Figure 7(b), δ2, δ1, and δ increase with increasinge, whereas Q/J, q1/J, and q2/J vary in different ways (i.e., as e increases, Q/J remains constant, q1/J increases, and q2/J decreases). In fact, q1/J and q2/J have a negative and symmetric correlation. From the enlarged view given in Figure 7(b), δ2 is still greater than δ1, which is still greater than δ. Based on (5a), the coefficient a is linearly proportional to the reciprocal of the cube of e. Thus, the coefficient a rapidly decreases with widening e, causing the broad changes in the cubic law-based Q/J. Considering that the basic value of the cubic law-based Q/J is not constant, in this group of tests, the changes in Q/J, q1/J, and q2/J in Figure 7(a) cannot describe the degree of nonlinearity of the fluid flow. δ2, δ1, and δ approach 100% for e > 5 mm. For example, when θ = 60°, as e widens from 5 to 20 mm, δ2 increases from 99.668 to 99.996%, δ1 increases from 99.008 to 99.982%, and δ increases from 98.627 to 99.977%. Note that varying e has a negligible effect on Q/J because the syringe pump in the tests can provide a stable flow rate for each hydraulic gradient. Unfortunately, the differences in the q1/J-e, q2/J-e, δ1-e, and δ2-e curves caused by the varying θ have not yet been presented and require more attention in future experimental studies.


## 3.5. Effects ofRr on the Nonlinearity of Fluid Flow

The tests in Figure9 were conducted on the specimens with e = 2 mm,Rr = 40, 100, 200, 300, and 400 mm, w = 5.7 mm, and θ = 60°, 90°, and 120°. There was a full flow in the influent fracture as well as both effluent fractures. In this group of tests, the predicted value of Q/J based on (1) is still the constant value 4.637 × 10−5. AsRr increases, the changes in Q/J and δ for all three fractures for the cases of θ = 60°, 90°, and 120° are shown in Figure 9. As shown in Figures 9(a) and 9(b), whenJ = 10, q2/J increases and δ2 declines with increasingRr, especially forRr < 100 mm. However, the magnitudes of q1/J and δ1 have different responses asRr increases; as one parameter increases, the other decreases. In fact, q1/J and q2/J have a negative and symmetric correlation, similar to δ1 and δ2. Thus, increasingRr enhances the nonlinearity of fluid flow in fracture_1 but decreases that in fracture_2. The consequence is that the magnitude of these changes in q2/J and δ2 diminishes whenRr > 100 mm, similar to the result obtained by Li et al. [23]. Therefore, if the fracture length is long, the disturbance on the fluid flow resulting from the intersection would become negligibly small. Furthermore, the on-way resistance resulting from the long effluent fractures decreases the dominant effect of the intersection. AsRr increases, the inertial resistances of fracture_1 and fracture_2 caused by the fracture length and the intersection decrease. In addition, as the on the way resistance becomes more dominant, the nonlinearity of fluid flow approaches a constant value. Note that the nonlinearity of fluid flow in fracture_2 is still greater than fracture_1, which is still greater than that in the influent fracture. As indicated in the enlarged views of Figures 9(a) and 9(b), as θ increases, q1/J increases (or q2/J decreases) and δ1 decreases (or δ2 increases). In addition, varyingRr and θ have a negligible effect on the nonlinearity of fluid flow in the influent fracture, resulting in a constant value of δ.Figure 9
(a) Relationship betweenQ/J andRr for θ = 60°, 90°, and 120°; (b) relationship between δ andRr for θ = 60°, 90°, and 120°.
(a) (b)


## 4. Conclusions

In this study, a series of tests including on 182 one-inlet-two-outlet specimens that including different fracture patterns withθ = 30°, 60°, 90°, 120°, and 150°, e = 2, 3, 5, 10, and 20 mm,Rr = 40, 100, 200, 300, and 400 mm, and w = 5.7 mm were carried out with a fluid flow testing system. Since the syringe pump in the tests provided a stable flow rate for each hydraulic gradient, the inertial force of fluid flow at the inlet remained constant and the nonlinearity of the fluid flow in the influent fracture was fixed when e was stable. As a result, the effects of J, θ, e, andRr on the nonlinearities of fluid flow in effluent fracture_1 and fracture_2 were independently investigated with one-inlet-two-outlet specimens containing a single fracture intersection. The main conclusions are as follows.(1) The hydraulic gradientJ, intersecting angle θ, aperture e, and fracture lengthRr significantly affect the nonlinearity of fluid flow through the fracture specimens with a single intersection.(2) Forchheimer’s law also applies for large magnitudes ofQ (10-5–10-4) and large Re (104–105).(3) The inertial resistanceβ of fracture_2 increases with increasingJ, whereas β of fracture_1 decreases.(4) The nonlinearity of fluid flow in fracture_1 decreases with the increase inθ, whereas that of fracture_2 increases. The nonlinearities of fluid flow in fracture_1 and fracture_2 have a negative and symmetric correlation.(5) AsJ increases, the nonlinearities of fluid flow in fracture_1, fracture_2, and the influent fracture all increase, especially for the cases of 100 < J < 101. However, whenJ > 101, the increase in nonlinearity stops.(6) The nonlinearities of fluid flow in fracture_1, fracture_2, and the influent fracture sharply increase ase widens, especially for e < 5 mm.(7) AsRr increases, the nonlinearity of fluid flow in fracture_1 increases, whereas that of fracture_2 decreases.(8) With an increase inθ,J, e, orRr, the nonlinearity of fluid flow in fracture_2 is always greater than fracture_1, which is always greater than that of the influent fracture.However, in this study, only one-inlet-two-outlet specimens were tested with particular attention paid to the effects ofJ, θ, e, andRr on the nonlinearity of fluid flow at largeRe and large J. To clearly figure out the mechanism of the nonlinear fluid flow behaviour through an entire fractured rock mass, further experimental tests on other fundamental elements of DFNs, especially the two inlets and two outlets type and two inlets and one outlet type, are needed to be carried out. In addition, since the specific geometries caused by the real intersection shapes with round bulges and corners can also strongly affect the nonlinear flow behaviour, further experimental studies on a broad spectrum of intersection geometries are required to fully elucidate the nonlinear flow behaviour through a fractured intersection. Based on the test results obtained from fundamental elements of DFNs, numerical model which can more realistically reflect the nonlinear fluid flow behaviour through an entire fractured rock mass by reasonably realizing the fracture system information with DFNs can then be developed in the future.



---

*Source: 9352608-2018-03-07.xml*