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

#### 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 , fracture networks essentially consist of very small fractures and the percolation theory applies. On the other hand, when , flow is mostly channelled into longer fractures and fracture networks can be considered as a superposition of long fractures. Between these ranges (), 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 () versus flow rate () increases with 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 ( and 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 , intersecting angle (), aperture (), 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 , , , and* Rr* on the nonlinearity of fluid flow in both effluent fractures at large Re and large . 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 aperture [24, 25]. At sufficiently low , this calculation yields the well-known “cubic law” [26, 27] and “Darcy’s law” [28, 29].where is the total volumetric flow rate, which represents the flow rate of the influent fracture in this study, is the pressure gradient, is the width of the fracture, is the uniform aperture of the idealized smooth fracture, is the dynamic viscosity,* k* is the intrinsic permeability defined as [5, 27], and is the cross-sectional area equal to . Equation (1) predicts that is proportional to the cube of and that is linearly correlated with when Re is sufficiently small. The cubic law has been widely applied to rough-walled fractures in rock; however, in these applications, is replaced by the so-called hydraulic aperture [24, 30, 31]. In this paper, is equal to 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 when is low. Nonlinear flow occurs when, with an incremental increase in Re, 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:where and are coefficients, and , is the fluid density, and is called either the non-Darcy coefficient or the inertial resistance [5, 39, 40].

The hydraulic gradient is proportional to :where and are the hydraulic heads at each end of a fracture, and are the pressures at each end of a fracture, and is the effective path under the hydraulic gradient. Equation (2) can be written in a Forchheimer’s law form as follows:The term represents the linear gradient, and the term represents the nonlinear gradient.

Coefficients and in (4) are commonly written as follows:Equation (5a) states that is linearly proportional to the reciprocal of the cube of . If the values of and are measurable and constant, the value of can be calculated directly. In addition, in some studies, the value of is considered the intrinsic permeability of the fractured media [29, 41, 42]. The coefficient 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]: where 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 velocity can be reduced to a one-dimensional flow velocity along a fracture intersection:

##### 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°, = 2, 3, 5, 10, and 20 mm, and* Rr* = 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 ,* Rr*, and were precisely measured by a digital Vernier caliper. Note that the 5.7 mm width was controlled by the length of the drill employed to cut the PMMA. The error of was ±0.1 mm (i.e., = 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.

**(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 Figure 2(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} m^{3}/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/m^{3} and a dynamic viscosity of 0.001 Pa·s at a room temperature of 20°C.

**(a)**

**(b)**

#### 3. Results and Discussions

As shown in Figure 3, 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.

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

The tests in Figure 4 were conducted on specimens with = 2 mm,* Rr* = 40 mm, = 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 is linearly correlated with , which can be expressed as . The relationship between and agrees with Forchheimer’s law. In the tests, as varies from 1.184 to 31.270, varies from to m^{3}/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 is 10^{−7}–10^{−6} and the magnitude of Re is 10^{1}–10^{2}. In addition, in this paper, it is demonstrated that Forchheimer’s law also applies for larger magnitudes of (–) and Re (10^{4}–10^{5}).

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

The variations of , , and as well as the cubic law-based with increasing for the cases of = 30°, 60°, 90°, 120°, and 150° are shown in Figures 4(b)–4(f). is quadratically correlated with , , and with a correlation coefficient > 0.955. When* J* = 10^{0}–10^{1}, as shown in Figure 4 the -, -, and - curves deviate from straight lines (cubic law) and the deviations increase with increasing . 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 -, -, and - curves from cubic law line are all in a decreasing order. In this group of tests, and are considered to be constant in (5a) and (5b); therefore, the coefficient of the linear term in (4) is calculated as the constant value 2.1465 × 10^{4}. Table 1 provides the coefficient of the nonlinear term in (4) for the cases of = 30°, 60°, 90°, 120°, and 150°. The coefficient can quantify the deviations of the -, -, and - curves. The coefficient of fracture_1 decreases with increasing *θ*, especially when ≥ 90° (i.e., as increases from 30° to 90°, decreases from 2.9477 × 10^{9} to 2.2524 × 10^{9}). In contrast, the coefficient of fracture_2 increases with increasing , especially when ≥ 90° (i.e., as increases from 30° to 90°, increases from 2.1500 × 10^{10} to 5.0517 × 10^{10}). However, the coefficient of the influent fracture only slightly changes with increasing (i.e., as increases from 30° to 150°, varies from 1.4966 × 10^{9} to 1.4981 × 10^{9}). This emphasizes that the magnitude of the coefficient of fracture_1 (or fracture_2) remains on the order of 10^{9} (or 10^{10}) as *θ* changes*. *, , and are measurable and assumed to be constant in (5a), so that the coefficient 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 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 with different specimens. Hence, for test cases with different , the - curves overlap each other, as demonstrated in Figure 4(a). The fixed 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 . of fracture_2 increases with increasing* J*, whereas *β* of fracture_1 decreases. As varies, the change in is clearly similar to that of .

##### 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 of 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 is based on (4) and declines with increasing . Furthermore, the change in represents the degree of nonlinearity. Therefore, as increases, 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 is between 10^{−3} and 10^{−2} for the case of* e* = 2 mm. In this paper, the fluid flow within the fractures is certainly in the nonlinear regime, since* J* = 10^{1}–10^{2}. The relative change in , , is defined as follows:where is the value of calculated by solving the cubic law, is the value of calculated by analysing data from the experiments, and expresses how the true value of () deviates from the prediction made by the cubic law (). Similarly, and , the relative discrepancies of fracture_1 and fracture_2, respectively, are calculated as follows:where and are the volumetric flow rates within fracture_1 and fracture_2, respectively. In this fluid flow test system,By substituting (10) into (8), (9a), and (9b), the relationship among , , and can be expressed as follows:Equation (11) suggests that if *δ* remains a constant value, and 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 Figure 5 were conducted on the specimens with = 2 mm,* Rr* = 40 mm, = 5.7 mm, and = 30°, 60°, 90°, 120°, and 150°. There was a full flow in the influent fracture and both effluent fractures. As increases, the variations in and for the cases of = 30°, 60°, 90°, 120°, and 150° are indicated in Figure 5. As provided in Table 2, and vary within only a small range as varies over a large range (0° < *θ* < 180°) for the cases of* J* = 7, 12, 17, and 22. For* J* = 7, with *θ* increasing from 30° to 150°, 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 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, decreases and increases with increasing *θ*, especially when < 90°; as increases, and have opposite changes in magnitude compared with those of and . and are negatively and symmetrically correlated, similar to and . Thus, increasing *θ* enhances the nonlinearity of the fluid flow in fracture_2 and decreases that of fracture_1. Note that > > . 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, for* J* = 7, decreases from 85.300% to 83.235% and increases from 94.518% to 96.336%; for* J* = 22, decreases from 91.993% to 90.913% and increases from 96.979% to 97.967%. In fact, , , , and only slightly vary with varying *θ*. The fluid flow nonlinearities in fracture_2 and fracture_1 are considerably greater because and are both approaching 100% for the cases of* J* = 10^{1}–10^{2}. Moreover, the changes in and diminish when *θ* > 90° (i.e., when* J* = 7, for *θ* = 90°, = 83.335% and = 96.368%; for *θ* = 120°, = 83.638% and = 96.368%; and for *θ* = 150°, = 83.235% and = 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.

**(a)**

**(b)**

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

The tests in Figure 6 were conducted on the specimens with = 2 mm,* Rr* = 40 mm, = 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 , based on (1), is a constant value of 4.637 × 10^{−5}. With the changes of , 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 -*J* and *δ*-*J* for the cases of* J* = 10^{−5}–10^{0} from Li et al. [23]. For the cases of* J* = 2–32, the changes in the -*J* and *δ*-*J* curves are exactly equal to the trendlines of Li et al., and the trendlines of -*J* and *δ*-*J* are extended to the cases of* J* = 10^{0}–10^{1}. As shown in Figures 6(a) and 6(b), when 2 ≤ ≤ 35 and increases, , , and decrease whereas , , and *δ* increase. Therefore, increasing strengthens the nonlinearity of fluid flow in all fractures. An increase in 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 10^{0} <* J* < 10^{1}. However, the magnitude of the change decreases when* J* > 10^{1}. For example, for the case of *θ* = 120°, when 10^{0} <* J* < 10^{1}, *δ* increases from 65.959 to 81.557% as increases from 2.736 to 8.167; when* J* > 10^{1}, *δ* increases from 85.872 to 91.061% as increases from 13.512 to 34.201. As indicated in the enlarged views of Figures 6(a) and 6(b), as increases, increases (or decreases) and decreases (or increases). Note that > > *δ*. 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 -*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 -*J* or *δ*-*J* influent fracture curves with different *θ*.

**(a)**

**(b)**

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

The tests in Figure 7 were conducted on the specimens with* e* = 2, 3, 5, 10, and 20 mm,* Rr* = 100 mm, = 5.7 mm, and = 60°, 90°, and 120°. There was a full flow within all fractures for the cases with* e* = 2, 3, and 5 mm. However, for the cases of* e* = 10 and 20 mm, the fluid flow in fracture_2 did not completely fill the open fracture space. For the case of* e* = 10 mm, in fracture_2, there were many bubbles that had diameters greater than 1 mm. For the case of* J* ≈ 10, as *θ* increases from 60° to 120°, decreased from 2.831 × 10^{−5} to 1.930 × 10^{−5} m^{3}/s; 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 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 = 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 = 20 mm, 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 of influent fracture was 0.723 m/s (for *θ* = 60°), 0.660 m/s (for *θ* = 90°), and 0.729 m/s (for *θ* = 120°). Hence, of fracture_2 sharply decreases with increasing *θ*, especially at = 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 of* e* = 20 mm and *θ* = 120°. In fact, for* e* = 20 mm and *θ* = 120°, a small volume of fluid turns back into fracture_2 due to eddy_3 (see Figure 8(c)).

**(a)**

**(b)**

**(a)**

**(b)**

**(c)**

In this paper, the relationship between and agrees with the Forchheimer’s law, though fluid flow in fracture_2 is in a non-full-flow state. With changing* e*, and *δ* for the cases of *θ* = 60°, 90°, and 120° are shown in Figure 7. As shown in Figure 7(b), , , and *δ* increase with increasing* e*, whereas , , and vary in different ways (i.e., as increases, remains constant, increases, and decreases). In fact, and have a negative and symmetric correlation. From the enlarged view given in Figure 7(b), is still greater than , which is still greater than . Based on (5a), the coefficient is linearly proportional to the reciprocal of the cube of . Thus, the coefficient rapidly decreases with widening , causing the broad changes in the cubic law-based . Considering that the basic value of the cubic law-based is not constant, in this group of tests, the changes in , , and in Figure 7(a) cannot describe the degree of nonlinearity of the fluid flow. , , and approach 100% for > 5 mm. For example, when *θ* = 60°, as widens from 5 to 20 mm, increases from 99.668 to 99.996%, increases from 99.008 to 99.982%, and increases from 98.627 to 99.977%. Note that varying has a negligible effect on because the syringe pump in the tests can provide a stable flow rate for each hydraulic gradient. Unfortunately, the differences in the -*e*, -*e*, -*e*, and -*e* curves caused by the varying *θ* have not yet been presented and require more attention in future experimental studies.

##### 3.5. Effects of* Rr* on the Nonlinearity of Fluid Flow

The tests in Figure 9 were conducted on the specimens with = 2 mm,* Rr* = 40, 100, 200, 300, and 400 mm, = 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 based on (1) is still the constant value 4.637 × 10^{−5}. As* Rr* increases, the changes in 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), when* J* = 10, increases and declines with increasing* Rr*, especially for* Rr* < 100 mm. However, the magnitudes of and have different responses as* Rr* increases; as one parameter increases, the other decreases. In fact, and have a negative and symmetric correlation, similar to and . Thus, increasing* Rr* 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 and diminishes when* Rr* > 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. As* Rr* 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, increases (or decreases) and decreases (or increases). In addition, varying* Rr* and *θ* have a negligible effect on the nonlinearity of fluid flow in the influent fracture, resulting in a constant value of *δ*.

**(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°, = 2, 3, 5, 10, and 20 mm,* Rr* = 40, 100, 200, 300, and 400 mm, and = 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 was stable. As a result, the effects of , , , and* Rr* 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 gradient , intersecting angle *θ*, aperture , and fracture length* Rr* significantly affect the nonlinearity of fluid flow through the fracture specimens with a single intersection.

(2) Forchheimer’s law also applies for large magnitudes of (–) and large Re (10^{4}–10^{5}).

(3) The inertial resistance of fracture_2 increases with increasing* J*, 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) As increases, the nonlinearities of fluid flow in fracture_1, fracture_2, and the influent fracture all increase, especially for the cases of 10^{0} < < 10^{1}. However, when* J* > 10^{1}, the increase in nonlinearity stops.

(6) The nonlinearities of fluid flow in fracture_1, fracture_2, and the influent fracture sharply increase as widens, especially for < 5 mm.

(7) As* Rr* increases, the nonlinearity of fluid flow in fracture_1 increases, whereas that of fracture_2 decreases.

(8) With an increase in *θ*,* J*, , or* Rr*, 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 of , , , and* Rr* on the nonlinearity of fluid flow at large* Re* and large . 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.

#### Conflicts of Interest

The authors declare that they have no conflicts of interest.