Nonlinear Flow Characteristics of a System of Two Intersecting Fractures with Different Apertures

The nonlinear flow regimes of a crossed fracture model consisting of two fractures have been investigated, in which the influences of hydraulic gradient, surface roughness, intersecting angle, and scale effect have been taken into account. However, in these attempts, the aperture of the two crossed fractures is the same and effects of aperture ratio have not been considered. This study aims to extend their works, characterizing nonlinear flow through a system of two intersecting fractures with different apertures. First, three experiment models with two fractures having different apertures were established and flow tests were carried out. Then, numerical simulations by solving the Navier-Stokes equations were performed and the results compared with the experiment results. Finally, the effects of fracture aperture on the critical pressure difference and the ratio of hydraulic aperture to mechanical aperture were systematically analyzed. The results show that the numerical simulation results agree well with those of the fluid flow tests, which indicates that the visualization techniques and the numerical simulation code are reliable. With the increment of flow rate, the pressure difference increases first linearly and then nonlinearly, which can be best fitted using Forchheimer’s law. The two coefficients in Forchheimer’s law decrease with the increasing number of outlets. When increasing fracture aperture from 3 mm to 5 mm, the critical pressure difference increases significantly. However, when continuously increasing fracture aperture from 5 mm to 7 mm, the critical pressure difference changes are negligibly small. The ratio of hydraulic aperture to mechanical aperture decreases more significantly for a fracture that has a larger aperture. Increasing fracture aperture from 5 mm to 7 mm, that has a negligibly small effect on the critical pressure difference will however significantly influence the ratio of hydraulic aperture to mechanical aperture.


Introduction
The assessment of hydraulic properties of deep fractured rock masses has attracted attention in many fields, such as CO 2 sequestration [1,2], enhanced oil recovery [3,4], and groundwater use [5][6][7][8][9]. In tight rocks (i.e., granite and basalt), the permeability of fractures is much larger than that of the rock matrix, and as a result, the rock matrix is assumed to be impermeable [10][11][12][13]. This assumption is commonly adopted when calculating the permeability of highly fractured rock masses using modelling techniques such as discrete fracture network (DFN) models [14][15][16]. Another assumption is that the fluid flow through fractures is in the linear regime and follows the cubic law in which the flow rate is linearly proportional to the pressure difference [10,[17][18][19][20]. However, in the high-pressure packer tests and/or karst systems, the fluid flow is in the nonlinear regime where the cubic law is not suitable [21,22]. In such a case, the permeability/conductivity of a fractured rock mass will be underestimated when still using the cubic law [21]. Therefore, the nonlinear flow characteristics of fractured rock masses should be investigated.
The nonlinear flow regimes through single fractures that are the simplest element of the complex fracture networks, have been extensively studied in the previous works. It has been found that rougher fractures can more significantly increase the length of flow paths and induce energy losses, which subsequently decreases the permeability [23,24]. Some studies have shown that the contact including contact ratio and contact shape can robustly change the tortuosity and flow paths of fluid, and then influence the nonlinear flow properties [25,26]. Since the rock masses are commonly moved along a fault when an earthquake occurs, the effect of shear displacement on the onset of nonlinear flow has been estimated both experimentally and numerically [27][28][29][30]. Besides the nonlinear flow through single fractures, the fracture interaction is another main cause for the nonlinear flow because the fluid redistributes at the interaction. Kosakowski and Berkowitz [21] depicted that the fluid flow transits from the linear regime to the nonlinear regime when Re > 10, by numerically solving the complex Navier-Stokes equations. The influences of hydraulic gradient, surface roughness, intersecting angle, and scale effect on the nonlinear flow characteristics at the intersections have been both experimentally and numerically investigated [31][32][33]. They found that the nonlinearity can be negligible when the ratio of fracture aperture to fracture length is smaller than 0.001, and the hydraulic gradient that is less than 0.001 plays a sufficiently small role on the nonlinear flow. However, in these attempts, the aperture of the two crossed fractures is the same and effects of aperture ratio have not been considered, if any. According to the cubic law, the flow rate through the fracture is linearly proportional to the cube of aperture [34]. A small variation in fracture aperture can induce a significant change in the flow rate. Many previous studies have contributed to the investigations of effect of fracture aperture distribution on the macroscopic flow rate through complex fracture networks [10,35,36], while other studies focused on the fracture aperture distribution [37][38][39] and linked it with fracture length [40][41][42][43][44][45][46]. Therefore, fracture aperture is a very important parameter that dominates the fluid flow. For a system of two intersecting fractures that is a simple fracture network, the aperture ratio of the two fractures would also be of special importance for investigating the hydraulic properties such as nonlinear flow regimes.
The present study extends the works of Li et al. [32] and Wu et al. [33] to study the nonlinear flow properties of a crossed fracture model consisting of two fractures with different apertures. First, three experiment models, with two fractures having different apertures, were established and flow tests were carried out. Then, numerical simulations by solving the Navier-Stokes equations were performed and the results compared with the experiments, which verified the validity of the simulation. Finally, the effects of fracture aperture ratio on the critical pressure difference and the ratio of hydraulic aperture to mechanical aperture were systematically analyzed.

Theoretical Background
The flow of incompressible steady Newtonian fluid is governed by the Navier-Stokes equations, written as [47]: where u is the flow velocity tensor, ρ is the fluid density, P is the hydraulic pressure, T is the shear stress tensor, t is the time, and f is the body force tensor. For a laminar flow in single fractures, the convective acceleration terms, (u · ∇)u, are the only source of nonlinearity, which are strongly affected by the geometries of void spaces resulting from complex surface roughness of fractures [48]. When the flow rate is sufficiently small, the Navier-Stokes equations can reduce to the cubic law, in which the flow rate is linearly proportional to the pressure gradient and is commonly adopted when simulating fluid flow in complex fracture networks [14,49,50], as follows: where Q is the flow rate, w is the fracture width, e is the hydraulic aperture, µ is the dynamic viscosity, and ∇P is the pressure gradient. Note that Equation (2) is applicable when the flow rate is sufficiently small. When the flow rate is large (i.e., Re > 10), the Q is no longer linearly proportional to ∇P [51]. For describing the nonlinear relationship between Q and ∇P, the famous Forchheimer's law is proposed as [34,52,53]: where A and B are two coefficients that represent the linear and nonlinear terms, respectively. To estimate the nonlinearity of fluid flow, a nonlinear factor E is defined as [54]: The E contributes to the ratio of pressure drop caused by the nonlinear term to the total pressure drop. Many critical values have been assigned to E to quantify the onset of nonlinear flow, such as E = 0.01 [55], 0.05 [56], and 0.1 [28,51], among which E = 0.1 is the most adopted and is used in the present study.

Specimen Preparation
One glass plate that has a size of 50 cm × 50 cm × 0.5 cm was cut into two parts in which there are two intersected fractures as shown in Figure 1. The intersecting angle of the two fractures is 60 degrees. Since the glass is a brittle material, the fractures deviated from the pre-designed sketch, resulting in rough fracture surfaces. During the cutting, some parts of the glass were damaged and removed from the model, making the two walls of a fracture ill-mated. This is also commonly observed in nature, for example, in shallow rock masses where mineral dissolution occurs [57]. The mean mechanical aperture (E 0 ) of each fracture was designed by shifting a separated glass plate. However, the shifting was handled by hand. Therefore, the accurate E 0 should be calculated again after the model is manufactured. The dye solute was injected into the model and then all the inlet and outlets were closed. A high-resolution charged coupled device (CCD) camera (Optem Inc., New York, NY, USA) was then utilized to capture each part of the model. Finally, the geometries of the fractures were re-plotted in AutoCAD. As a result, three experimental models (Cases 1-3) were made; their geometries included tortuous segment length (L t ), aperture, and surface roughness and are shown in Table 1, in which the mechanical aperture is represented by E 0 and fracture surface roughness is represented by root-mean-square of first derivative of asperity height, Z 2 , and is calculated using the following equation [58,59]: where x i and z i are the coordinates of the fracture surface profile, and M is the number of sampling points along the length (i.e., x coordinate) of a fracture.     Figure 2 shows the experimental system of the fluid flow test. The fluid was first stored in a syringe pump with an accuracy of ±0.001 mL/min. Then the valve in the left side of the syringe pump was closed and that in the right side of the syringe pump was opened. The fluid was then injected into the model with a fixed flow rate. The pressure difference between the inlet and outlet was measured using a differential pressure gauge with an accuracy of ±10 Pa. The fluid out of the model was collected using a collector; its weight was measured using an electronic balance with an accuracy of ±0.01 g and was recorded in a computer that was linked with the electronic balance using a R232 data line.

Test System
It is assumed that the fluid is pure water and the density is 998.2 kg/m 3 at a room temperature 20 • C. The dynamic viscosity is 0.001 Pa·s. The glass is impermeable and fluid only flows along the connected fractures.  Figure 2 shows the experimental system of the fluid flow test. The fluid was first stored in a syringe pump with an accuracy of ±0.001 mL/min. Then the valve in the left side of the syringe pump was closed and that in the right side of the syringe pump was opened. The fluid was then injected into the model with a fixed flow rate. The pressure difference between the inlet and outlet was measured using a differential pressure gauge with an accuracy of ±10 Pa. The fluid out of the model was collected using a collector; its weight was measured using an electronic balance with an accuracy of ±0.01 g and was recorded in a computer that was linked with the electronic balance using a R232 data line.

Test System
It is assumed that the fluid is pure water and the density is 998.2 kg/m 3 at a room temperature 20 °C. The dynamic viscosity is 0.001 Pa•s. The glass is impermeable and fluid only flows along the connected fractures.

Test Procedure
Before the test, all the inflow and outflow tanks, tubes, and fractures were filled with water. The air bubbles within them were removed using a vacuum pump. Although there were four tanks, only one tank was used as the inflow tank and the other three tanks were used as the outflow tanks. For each model, six cases were tested in which the combinations of outlets were different, including outlet 2, 3, 4, 2 and 3, 2 and 4, 2 and 3 and 4. The flow rate of fluid from the syringe pump was set varying from 0 to 100 mL/min. When the fluid flow was in the steady state, the pressure difference between the inlet and outlet was recorded. Finally, the flow rate-pressure difference (Q~ΔP) curves were plotted.

Numerical Models and Meshing
To verify the validity of the visualization technique used for obtaining fracture geometries, three numerical models were established by replotting the geometries of experimental models corresponding to Cases 1-3 in Table 1. The three numerical models have the same model lengths and

Test Procedure
Before the test, all the inflow and outflow tanks, tubes, and fractures were filled with water. The air bubbles within them were removed using a vacuum pump. Although there were four tanks, only one tank was used as the inflow tank and the other three tanks were used as the outflow tanks. For each model, six cases were tested in which the combinations of outlets were different, including outlet 2, 3, 4, 2 and 3, 2 and 4, 2 and 3 and 4. The flow rate of fluid from the syringe pump was set varying from 0 to 100 mL/min. When the fluid flow was in the steady state, the pressure difference between the inlet and outlet was recorded. Finally, the flow rate-pressure difference (Q~∆P) curves were plotted.

Numerical Models and Meshing
To verify the validity of the visualization technique used for obtaining fracture geometries, three numerical models were established by replotting the geometries of experimental models corresponding to Cases 1-3 in Table 1. The three numerical models have the same model lengths and fracture geometries such as aperture and fracture surface roughness with those of experimental models. However, as shown in Table 1, the apertures of the segments of fractures for Cases 1-3 are different and do not show some rules, which therefore cannot be utilized for parametric study. To solve this issue, we subsequently established another three numerical models and their enlarged views of intersection are shown in Figure 3. The fracture surface roughness and segment length are the same as those in the Cases 1-3, but the apertures of the fractures are different as tabulated in Table 1. The apertures of segments 1 and 3 are fixed as 3 mm, yet the apertures of segments 2 and 4 are changed from 3 mm (see Figure 3a), to 5 mm (see Figure 3b), and to 7 mm (see Figure 3c), respectively. fracture geometries such as aperture and fracture surface roughness with those of experimental models. However, as shown in Table 1, the apertures of the segments of fractures for Cases 1-3 are different and do not show some rules, which therefore cannot be utilized for parametric study. To solve this issue, we subsequently established another three numerical models and their enlarged views of intersection are shown in Figure 3. The fracture surface roughness and segment length are the same as those in the Cases 1-3, but the apertures of the fractures are different as tabulated in Table  1. The apertures of segments 1 and 3 are fixed as 3 mm, yet the apertures of segments 2 and 4 are changed from 3 mm (see Figure 3a), to 5mm (see Figure 3b), and to 7 mm (see Figure 3c), respectively. The numerical models were plotted using software AutoCAD (Autodesk Inc., California, United States, AutoCAD2010), and were then exported in the sat file format. Software Gambit (Fluent Inc., New York, NY, USA, Gambit2.4.6) was used for meshing using hexahedral cells. An interval of approximately 0.5 mm was set. As a result, the total number of hexahedral cells for the six cases ranged from 212,320 to 323,060 due to the difference in fracture aperture. Finally, the stored file from Gambit was imported into finite volume method (FVM) software FLUENT for flow calculation, in which the Navier-Stokes equations were solved. For the fracture aperture, we chose values of 3-7 mm, which is much larger than that in practical engineering. This is because the meshing number (or number of hexahedral cells) should be one hundred times larger than the current models (with number of hexahedral cells varying from 212,320 to 323,060 as shown in Table 1), if the fracture aperture range from 0.3-0.7 mm and the same number (i.e., about 6 in the present study) of layer is kept along the aperture direction. Using a common computer, this meshing cannot be generated and the fluid flow cannot be efficiently modeled by solving the Naver-Stokes equations. Considering this, the magnitude of fracture aperture that is several millimeters should be accepted.

Boundary Conditions and Flow Calculation
Assumptions were made that the rock matrix would be impermeable with fluid flowing along connected fractures. Therefore, a constant hydraulic pressure was applied on the inlet and the outlet connected to the air, and the flow rate from the outlet was recorded. For each case in Table 1, the fluid flowed through the model from inlet (inlet 1) to outlet (i.e., outlet 2, 3, 4, 2 and 3, 2 and 4, and 2, 3 and 4). The calculated Q~ΔP curves for Cases 1-3 were compared with those of experimental results, and the calculated Q~ΔP curves for Cases 4-6 were used for deep analysis.

Comparison between Flow Tests and Numerical Simulations
Since the geometries of the models used in flow tests and numerical simulations for Cases 1-3 in Table 1 were the same, their results, such as the Q~ ΔP curves with different combinations of outlets (i.e., outlet 2, 3, 4, 2 and 3, 2 and 4, and 2 and 3 and 4) were compared as shown in Figures 4-6. Both The numerical models were plotted using software AutoCAD (Autodesk Inc., California, United States, AutoCAD2010), and were then exported in the sat file format. Software Gambit (Fluent Inc., New York, NY, USA, Gambit2.4.6) was used for meshing using hexahedral cells. An interval of approximately 0.5 mm was set. As a result, the total number of hexahedral cells for the six cases ranged from 212,320 to 323,060 due to the difference in fracture aperture. Finally, the stored file from Gambit was imported into finite volume method (FVM) software FLUENT for flow calculation, in which the Navier-Stokes equations were solved. For the fracture aperture, we chose values of 3-7 mm, which is much larger than that in practical engineering. This is because the meshing number (or number of hexahedral cells) should be one hundred times larger than the current models (with number of hexahedral cells varying from 212,320 to 323,060 as shown in Table 1), if the fracture aperture range from 0.3-0.7 mm and the same number (i.e., about 6 in the present study) of layer is kept along the aperture direction. Using a common computer, this meshing cannot be generated and the fluid flow cannot be efficiently modeled by solving the Naver-Stokes equations. Considering this, the magnitude of fracture aperture that is several millimeters should be accepted.

Boundary Conditions and Flow Calculation
Assumptions were made that the rock matrix would be impermeable with fluid flowing along connected fractures. Therefore, a constant hydraulic pressure was applied on the inlet and the outlet connected to the air, and the flow rate from the outlet was recorded. For each case in Table 1, the fluid flowed through the model from inlet (inlet 1) to outlet (i.e., outlet 2, 3, 4, 2 and 3, 2 and 4, and 2, 3 and 4). The calculated Q~∆P curves for Cases 1-3 were compared with those of experimental results, and the calculated Q~∆P curves for Cases 4-6 were used for deep analysis.

Comparison between Flow Tests and Numerical Simulations
Since the geometries of the models used in flow tests and numerical simulations for Cases 1-3 in Table 1 were the same, their results, such as the Q~∆P curves with different combinations of outlets (i.e., outlet 2, 3, 4, 2 and 3, 2 and 4, and 2 and 3 and 4) were compared as shown in Figures 4-6. Both results show that with the increment of Q, ∆P increases first linearly and then nonlinearly. This is because when the Q is small, the inertial effects can be negligible and fluid flow obeys the cubic law (see Equation (2)), in which Q is linearly proportional to ∆P. When the Q is large, the inertial force is significantly large with respect to the viscous force and it cannot be negligible. As a result, with continuously increasing Q, ∆P increases nonlinearly versus Q (see Equation (3)) and the increasing rate of ∆P increases due to energy losses induced by eddies. With the increase in the number of outlets, for examples from 1 outlet (outlet 2, 3, or 4) to 3 outlets (outlets 2 and 3 and 4), for the same Q, ∆P decreases. A possible reason is that the more outlets correspond to the larger flow volume, and as a result, the ∆P decreases according to Darcy's law. Even though the simulated results agree well with the tested results in a whole such as Figure 4c,e,f, Figure 5b,f and Figure 6b,c, their results deviate for some cases such as Figure 4a,b and Figure 6a,c, in which the calculated ∆P is either underestimated (Figure 4a) or overestimated (Figure 4b).
The deviations arise from the model manufacture in which the glass was not cut strictly along the z-direction as shown in Figure 1 and there existed an angle between the destroyed direction and the z-direction. Thus, in AutoCAD, only the surface of the fractures can be replotted and the void space of fractures along the z-direction cannot be accurately characterized. As a result, the geometries of the fractures in the numerical models are slightly different from those in the experimental models.
Processes 2018, 6, x FOR PEER REVIEW 7 of 24 results show that with the increment of Q, ΔP increases first linearly and then nonlinearly. This is because when the Q is small, the inertial effects can be negligible and fluid flow obeys the cubic law (see Equation (2)), in which Q is linearly proportional to ΔP. When the Q is large, the inertial force is significantly large with respect to the viscous force and it cannot be negligible. As a result, with continuously increasing Q, ΔP increases nonlinearly versus Q (see Equation (3)) and the increasing rate of ΔP increases due to energy losses induced by eddies. With the increase in the number of outlets, for examples from 1 outlet (outlet 2, 3, or 4) to 3 outlets (outlets 2 and 3 and 4), for the same Q, ΔP decreases. A possible reason is that the more outlets correspond to the larger flow volume, and as a result, the ΔP decreases according to Darcy's law. Even though the simulated results agree well with the tested results in a whole such as Figures 4c,e,f, 5b,f and 6b,c, their results deviate for some cases such as Figures 4a,b and 6a,c, in which the calculated ΔP is either underestimated (Figure 4a) or overestimated (Figure 4b). The deviations arise from the model manufacture in which the glass was not cut strictly along the z-direction as shown in Figure 1 and there existed an angle between the destroyed direction and the z-direction. Thus, in AutoCAD, only the surface of the fractures can be replotted and the void space of fractures along the z-direction cannot be accurately characterized. As a result, the geometries of the fractures in the numerical models are slightly different from those in the experimental models.

Parametric Studies
In this section, parametric studies are performed using the numerical models as shown for Cases 4~6 in Table 1, in which the aperture of one fracture is fixed as 3 mm and the aperture of another fracture varies from 3 mm, to 5 mm, and to 7 mm, respectively, but their intersecting angle is still fixed as 60 degrees.  (3) decrease. For the cases with the same number of outlets (i.e., outlet 2, 3, and 4), the Q~ΔP curves and coefficients A and B are close to each other, indicating that when the flow volume is in the same order of magnitude, the hydraulic properties are close. When increasing apertures of segments 2 and 4 from 3 mm to 5 mm (from Figure 7 to Figure 8), the Q~ΔP curves move rightwards significantly. However, when increasing apertures of segments 2 and 4 from 5 mm to 7 mm (from Figure 8 to Figure 9), the Q~ΔP curve variation is negligibly small. This indicates that when the Q~ΔP curves are more sensitive to fractures that have smaller apertures.

Effect of fracture Aperture on E~ΔP Curves
The nonlinear factor E defined in Equation (3) is adopted to quantify the nonlinearity of fluid flow in rock fractures [51,60]. Figures 10-12 Figure 11a) to 2 outlets (i.e., outlet 2 and 4 in Figure 11e), or to 3 outlets (i.e., outlet 2 and 3 and 4 in Figure 11f), for the same ΔP, E decreases. This is because a greater number of outlets indicate the larger flow volume, and the larger fracture aperture (or flow volume) corresponds to a weaker nonlinearity of fluid flow that is represented by a smaller E [55]. When increasing apertures of segments 2 and 4 from 3 mm to 5 mm (from Figure 10 to Figure 11), the E~ΔP curves move downwards significantly. However, when increasing apertures of segments 2 and 4 from 5 mm to 7 mm (from Figure 11 to Figure 12), the Q~ΔP curve variation is negligibly small. This indicates that

Parametric Studies
In this section, parametric studies are performed using the numerical models as shown for Cases 4~6 in Table 1, in which the aperture of one fracture is fixed as 3 mm and the aperture of another fracture varies from 3 mm, to 5 mm, and to 7 mm, respectively, but their intersecting angle is still fixed as 60 degrees. For the cases with the same number of outlets (i.e., outlet 2, 3, and 4), the Q~∆P curves and coefficients A and B are close to each other, indicating that when the flow volume is in the same order of magnitude, the hydraulic properties are close. When increasing apertures of segments 2 and 4 from 3 mm to 5 mm (from Figure 7 to Figure 8), the Q~∆P curves move rightwards significantly. However, when increasing apertures of segments 2 and 4 from 5 mm to 7 mm (from Figure 8 to Figure 9), the Q~∆P curve variation is negligibly small. This indicates that when the Q~∆P curves are more sensitive to fractures that have smaller apertures.

Effect of fracture Aperture on E~∆P Curves
The nonlinear factor E defined in Equation (3) is adopted to quantify the nonlinearity of fluid flow in rock fractures [51,60]. Figures 10-12 exhibit the variations of E versus ∆P hat ranges from 0 to 1000 Pa. The results show that with the increment ∆P, E increases first fast and then gently, and finally approaches constant values. With increasing the number of outlets, for examples from 1 outlet (i.e., outlet 2 in Figure 11a) to 2 outlets (i.e., outlet 2 and 4 in Figure 11e), or to 3 outlets (i.e., outlet 2 and 3 and 4 in Figure 11f), for the same ∆P, E decreases. This is because a greater number of outlets indicate the larger flow volume, and the larger fracture aperture (or flow volume) corresponds to a weaker nonlinearity of fluid flow that is represented by a smaller E [55]. When increasing apertures of segments 2 and 4 from 3 mm to 5 mm (from Figure 10 to Figure 11), the E~∆P curves move downwards significantly. However, when increasing apertures of segments 2 and 4 from 5 mm to 7 mm (from Figure 11 to Figure 12), the Q~∆P curve variation is negligibly small. This indicates that the E~∆P curves are more sensitive to fractures that have smaller apertures, which have the same variation trend with those of Q~∆P curves as illustrated in Section 5.2.1.
Processes 2018, 6, x FOR PEER REVIEW 11 of 24 the E~ΔP curves are more sensitive to fractures that have smaller apertures, which have the same variation trend with those of Q~ΔP curves as illustrated in Section 5.2.1.           Figure 13 depicts the relationships between the critical hydraulic difference (ΔPc) that corresponds to E = 0.1 and aperture of segments 2 and 4. For the six cases considered, with the increment of aperture of segments 2 and 4 from 3 mm to 5 mm, the average ΔPc increases significantly from 38.84 Pa to 489.38 Pa. However, with continuously increasing aperture of segments 2 and 4 from 5 mm to 7 mm, the average ΔPc changes slightly, or even decreases to some extent. This indicates that  Figure 13 depicts the relationships between the critical hydraulic difference (∆P c ) that corresponds to E = 0.1 and aperture of segments 2 and 4. For the six cases considered, with the increment of aperture of segments 2 and 4 from 3 mm to 5 mm, the average ∆P c increases significantly from 38.84 Pa to 489.38 Pa. However, with continuously increasing aperture of segments 2 and 4 from 5 mm to 7 mm, the average ∆P c changes slightly, or even decreases to some extent. This indicates that ∆P c is more sensitive to a smaller aperture of segments 2 and 4. Before selecting the governing equation for modeling fluid flow in rock fractures, the applied pressure difference should be compared with ∆P c . When the applied pressure difference is less than ∆P c , the cubic law that can be simply solved should be selected. When the applied pressure difference is larger than ∆P c , the Navier-Stokes equations or the Forchheimer's law should be adopted to solve the nonlinear flow through rock fractures/fracture networks.

Processes 2018, 6, x FOR PEER REVIEW 17 of 24
ΔPc is more sensitive to a smaller aperture of segments 2 and 4. Before selecting the governing equation for modeling fluid flow in rock fractures, the applied pressure difference should be compared with ΔPc. When the applied pressure difference is less than ΔPc, the cubic law that can be simply solved should be selected. When the applied pressure difference is larger than ΔPc, the Navier-Stokes equations or the Forchheimer's law should be adopted to solve the nonlinear flow through rock fractures/fracture networks.

Effect of Fracture Aperture on e/E0~ΔP Curves
In Equation (2), e is the hydraulic aperture, and in Table 1, E0 is the mechanical aperture. Generally, to quantify the effect of parameters such as fracture surface roughness, Reynolds number and contact on the hydraulic properties of rock fractures, the ratio of hydraulic aperture to mechanical aperture is utilized [61][62][63]. Therefore, the influence of aperture of segments 2 and 4 on the variations of e/E0 was analyzed in this section and the results are presented in Figures 14-16. For each case, with increasing ΔP when ΔP is large, the fluid flow enters the nonlinear flow regime and energy losses due to the formation of eddies. As a result, the value of ΔQ/ΔP decreases and then the e that is backcalculated according to Equation (2) decreases. Thus, e/E0 decreases, and the decreasing rate decreases as ΔP increases. In Figure 14, since the apertures of the two crossed fractures are both 3 mm, the variations of e/E0 versus ΔP have the same trends and orders of magnitude. When increasing the aperture of segments 2 and 4 to 5 mm, the results show that the e/E0 of segment 2, which is represented by outlet 2 in Figure 15a, decreases more significantly than that of segment 3, and the decreasing rate of e/E0 of segment 2 is smaller than that of segment 3. Since the apertures of segments 2 and 4 are the same, that is 5 mm, the variations of e/E0 of segments 2 and 4 are close to each other as shown in Figure 15b. Figure 15c shows the same results as that in Figure 15, in which the variations of e/E0 of segment that has a larger aperture (i.e., segments 2 and 4 that have an aperture of 5 mm) are more robust than one with a smaller aperture (i.e., segment 3 that has an aperture of 3 mm). With continuously increasing the aperture of segments 2 and 4 to 7 mm, the variations of e/E0 as shown in Figure 16 are close to those shown in Figure 15. The difference is that the decreasing degree of e/E0 of a segment that has a larger aperture is more significant. Therefore, even though increasing aperture of segments 2 and 4 from 5 mm to 7 had a negligibly small effect on the ΔPc, its influence on e/E0 cannot be negligible and should be considered. In Equation (2), e is the hydraulic aperture, and in Table 1, E 0 is the mechanical aperture. Generally, to quantify the effect of parameters such as fracture surface roughness, Reynolds number and contact on the hydraulic properties of rock fractures, the ratio of hydraulic aperture to mechanical aperture is utilized [61][62][63]. Therefore, the influence of aperture of segments 2 and 4 on the variations of e/E 0 was analyzed in this section and the results are presented in Figures 14-16. For each case, with increasing ∆P when ∆P is large, the fluid flow enters the nonlinear flow regime and energy losses due to the formation of eddies. As a result, the value of ∆Q/∆P decreases and then the e that is back-calculated according to Equation (2) decreases. Thus, e/E 0 decreases, and the decreasing rate decreases as ∆P increases. In Figure 14, since the apertures of the two crossed fractures are both 3 mm, the variations of e/E 0 versus ∆P have the same trends and orders of magnitude. When increasing the aperture of segments 2 and 4 to 5 mm, the results show that the e/E 0 of segment 2, which is represented by outlet 2 in Figure 15a, decreases more significantly than that of segment 3, and the decreasing rate of e/E 0 of segment 2 is smaller than that of segment 3. Since the apertures of segments 2 and 4 are the same, that is 5 mm, the variations of e/E 0 of segments 2 and 4 are close to each other as shown in Figure 15b. Figure 15c shows the same results as that in Figure 15, in which the variations of e/E 0 of segment that has a larger aperture (i.e., segments 2 and 4 that have an aperture of 5 mm) are more robust than one with a smaller aperture (i.e., segment 3 that has an aperture of 3 mm). With continuously increasing the aperture of segments 2 and 4 to 7 mm, the variations of e/E 0 as shown in Figure 16 are close to those shown in Figure 15. The difference is that the decreasing degree of e/E 0 of a segment that has a larger aperture is more significant. Therefore, even though increasing aperture of segments 2 and 4 from 5 mm to 7 had a negligibly small effect on the ∆P c , its influence on e/E 0 cannot be negligible and should be considered.

Conclusions
In the present study, fluid flow tests on three crossed fracture models were carried out, and the geometries of the fractures were accurately obtained using visualization techniques. Three numerical

Conclusions
In the present study, fluid flow tests on three crossed fracture models were carried out, and the geometries of the fractures were accurately obtained using visualization techniques. Three numerical models that corresponded to the three experimental models were established, and fluid flow was modelled by solving the Navier-Stokes equations. Finally, another three numerical models were established to perform parametric studies to quantify the onset of nonlinear flow and to systematically investigate the variations of the ratio of hydraulic aperture to mechanical aperture.
The results show that the numerical simulation results agree well with those of the fluid flow tests, indicating that the visualization techniques and the numerical simulation code are reliable. With the increment of flow rate, the pressure difference first increases linearly and then nonlinearly, which can be best fitted using Forchheimer's law. The two coefficients in Forchheimer's law decrease with the increasing number of outlets. With the increment of pressure difference, the nonlinear factor first increases significantly and then gently, and finally approaches constant values. The nonlinear factor-pressure difference curves are more sensitive to fractures that have smaller apertures. With the increment of pressure difference, the ratio of hydraulic aperture to mechanical aperture decrease and the decreasing rate decreases. The ratio of hydraulic aperture to mechanical aperture decreases more significantly for a fracture that has a larger aperture.
Note that the present study assumes that the rock matrix is impermeable. However, in practical engineering, the rock matrix is a porous media and fluid flows through it. Therefore, the interaction between a fracture and a rock matrix may play a significant role in the exchange of fluid. The present study adopts an aperture of 3-7 mm for rock fractures, which is much larger than that of fractures deep underground. Further, the influence of fracture length ratio and combinations of different inlets and outlets (i.e., two inlets and one outlet, two inlets and two outlets) on the nonlinear flow regimes are not taken into account. All of these should be the focus of our future works.

Data Availability
All the data used in the present study are available by connecting the corresponding author and/or the first author. Their emails are jiang@nagasaki-u.ac.jp (Yujing Jiang) and liuricheng@cumt.edu.cn (Richeng Liu), respectively.
Author Contributions: R.L. and Y.J. conceived and designed the experiments; R.L. performed the experiments; R.L. and H.J. analyzed the data; Y.J. and L.Y. contributed reagents/materials/analysis tools; R.L. wrote the paper.
Funding: This research received no external funding.