Experimental Investigation on the Hydraulic Characteristics of the Two-Phase Flow in the Intersecting Rock Fractures

,


Introduction
In the subsurface space, multiple phases exist in the rock fracture network, such as geothermal energy reservoirs, natural gas-oil reservoirs, coal seams, and shale gas reservoirs [1][2][3][4][5].The fluid-conducting capacity of the fracture network is the main concern in such engineering applications.Since the hydraulic properties of one phase of fluid in fractures will be totally different in the presence of another phase [6], a comprehensive understanding of the two-phase flow in fractures is necessary.Different from the single-phase flow, the two-phase flow in the fracture is influenced not only by the intrinsic properties of fractures, such as the fracture roughness and aperture, but also by the interactions between two phases [7][8][9][10].The two-phase flow in a single fracture is the first step to understanding the flow characteristics in the fracture network [11][12][13][14].The conventional approach to predict the two-phase flow in a single fracture is to extend Darcy's law for the single-phase flow to a generalized Darcy's equation for the two-phase flow [15][16][17], in which relative permeability is an important parameter.Relative permeability is defined as the ratio of the effective permeability of a certain phase in the two-phase flow to the permeability in the single-phase flow.It is the parameter to indicate the flow interference between two phases.With an exact model of relative permeability and permeability, the relationship between the pressure field and the velocity field in the two-phase flow can be established.The main concern of the researchers is how the relative permeability evolves, and different models are proposed to describe its evolution [7,[18][19][20][21][22].
In the subsurface, single fractures intersect with each other and constitute the fracture network.Besides the hydraulic characteristics of fluid in single fractures, insight into that of fracture intersections is also essential.Studies on the single-phase flow indicate that the fluid flow in the intersecting fractures shows quite different characteristics from that in a single fracture due to the energy loss and interference induced by the fracture intersection [23,24].Wilson and Withspoon [25] studied the nonlinearity of fluid transport at the intersection of circular tubes with experiment, and their results indicate that the head loss induced by the intersection is about five times of the tube diameter when the Reynolds number is 100.Kosakowski and Berkowitz [26] numerically investigated the hydraulic properties of the intersecting fractures with Navier-Stokes equation, and their results indicate that when the Reynolds number is within 1 to 100, the inertial effect occurs and leads to a nonlinear flow.Su and Zhan [27] conducted water flow tests in crossed fractures and established the theoretical model for local loss of head at the intersection, and they proposed the correction factor for the theoretical model.All the abovementioned studies indicate that the hydraulic properties in the intersecting fractures have a nonlinearity and additional pressure drop and Darcy's law is generally not effective in such situations.However, the two-phase flow characteristics in the intersecting fractures are still not well understood.The generalized Darcy's equation is usually used to describe the two-phase flow, in which the relative permeability is the parameter to account for the flow interference between two phases, but the relative permeability does not account for the nonlinearity of the flow induced by the inertial effect or turbulence [10,28], which cannot be neglected in the intersecting fractures.
With a developed system, an experiment is conducted to investigate the two-phase flow characteristics in a 3D fracture model with smooth intersecting fractures.The transparency of the 3D fracture model makes it possible to capture the flow structures.The pressure drop characteristics are analyzed with the Lockhart-Martinelli model and Delhaye's model, which is initially proposed for the two-phase flow in pipes.This study presents a basis for further studies to understand the two-phase flow properties in the rock fracture network and the mechanism of groundwater inflow in engineering applications.

Experiment in the Intersecting Fractures
2.1.Experiment System.The experimental system is shown in Figure 1, which is composed of four subsystems, as introduced in the following: (1) The water supply subsystem.In this subsystem, the peristaltic pump is used for injecting water at a specified flow rate within 0-2000 mL/min.Since the pressure and flow rate of the injected water is always fluctuating, a pulse damper is connected to the pump to decrease the pulse in order to uniformly inject water into the 3D fracture model (2) The gas supply subsystem.In this subsystem, nitrogen is supplied from a gas cylinder.The gas pressure in the cylinder is as large as 10 MPa, so a pressure regulator is connected to the cylinder to decrease the gas pressure to be under 0.3 MPa in order to protect the mass flow controller.The mass flow controller can control the gas flow rate to be a constant value within 0-5000 mL/min (3) The 3D fracture model.This model is made of acrylic material and transparent, as shown in Figure 2. The schematic of the model from a vertical view (indicated by the blue arrow in Figure 2) is shown in Figure 3.The model has five smooth fractures: the apertures of fractures 1-5 are 0.6 mm, 0.3 mm, 0.3 mm, 0.6 mm, and 1 mm, respectively.The length of fracture 5 is 10 cm and the length of other four fractures is 5.78 cm.The model has two inlets connected to fracture 1 and fracture 2 and two outlets connected to fracture 3 and fracture 4. By opening or closing different inlets and outlets, different combinations of fractures can be established.In this experiment, two cases are adopted for testing, as shown in Figure 3, in which the red lines indicate the activated fractures, meaning that the flow of fluids only happens in these fractures, while the black lines indicate the fractures where inlets or outlets are closed, meaning that no flow occurs.The red arrow in Figure 2 shows the direction from which the photos were taken.
The presented configuration of fractures in this 3D fracture model is a simplified configuration of the real ones.In real cases, the fractures are with roughness and may be continuous after the intersection.These multiple factors influence the twophase flow in complicated ways.In order to avoid the effect of multiple factors and quantitatively study the effect of intersecting fractures, this simplified model is used for experiment (4) The measurement subsystem.This subsystem includes a camera, three pressure sensors, a separation bottle, an electronic balance, and a flowmeter.The flow structures are captured by the camera.One pressure sensor is connected to the gas inlet, one is connected to the water inlet, and one is connected to the outlet.Pressure data are recorded by a digital recorder (not shown in Figure 1).The separation bottle, which is put on the electronic balance, is used to separate water and gas.Gas flows out of the separation bottle and water remains in the bottle.The water mass is recorded by the electronic balance in real time, and the data are transmitted to the computer simultaneously to calculate the flow rate of water It should be noted that the water flow rate can be indicated in the peristaltic water pump, but without enough accuracy.That is why we set up an electronic balance to measure the water flow rate.The gas flow rate can be indicated by both the mass flow controller and the flowmeter.It shows that the gas flow rate indicated by the mass flow controller is identical to that of the flowmeter.Since the gas flow rate displayed by the mass flow controller is accurate enough, the flowmeter can be removed.
2 Geofluids 2.2.Experiment Procedures.As mentioned above, this experiment was conducted with two cases.In each case, the tests were conducted in the following procedures: (1) Connect all the units in the experiment system and adjust the levelness of the 3D fracture model with a level gauge (2) Conduct the preliminary test to check the tightness of the system and the accuracy of the measurement subsystem.The flowmeter connected to the separation bottle is used to check the accuracy of the gas flow rate displayed by the mass flow controller, and the electronic balance is used to check the accuracy of the water flow rate displayed by the peristaltic pump.The test results show that the accuracy of the peristaltic pump is not good, but the gas flow rate indicated by the flowmeter is identical to that displayed by the mass flow controller.However, we found that the flowmeter brought a large resistance to the flow of fluids.Since this would increase burden to the pump, the flowmeter was removed in the formal experiment in steps (3) and (4) (3) Conduct a single-phase flow test of water to obtain the single-phase hydraulic properties of the 3D model.Water was injected at different flow rates within 0-600 mL/min (4) Conduct the two-phase flow test.In each case, five groups of two-phase flow tests were conducted.In each group, the water flow rate was kept constant while the gas flow rate was increased step by step.
In each test, when the injection of water and gas started, wait for five minutes for the flow to reach a stable state.Then, the pressure and the water mass in the separation bottle were recorded with the pressure sensors and electronic balance.At the same time, the flow structures were captured by the camera.Since the pressure and the flow rates were always fluctuating, in each test, the data recording lasted one minute in order to obtain the average values; ten photos were taken to acquire the flow structures and average saturation in each test round

Model for Describing the Two-Phase Flow in Intersecting Fractures
As is known, in the single-phase flow, Darcy's law is not effective if the flow nonlinearity exists.This nonlinearity is caused by the inertial effect or turbulence, which can be induced by the fracture roughness or fracture intersection [23].Consequently, revised laws such as Forchheimer's law are proposed to account for the nonlinearity [29][30][31].Forchheimer's law can be described as equation ( 1): in which a is the parameter related to the aperture of the fracture and b is a parameter that indicates the nonlinearity of the flow.The equation is found to be applicable to the flow in porous media [32,33] and in rough fractures [34,35].It has been demonstrated by literature [23] that 3 Geofluids the single-phase flow in this 3D model used in our experiment is nonlinear and can be described by Forchheimer's law.
As for the two-phase flow, the Lockhart-Martinelli model considers the flow nonlinearity.This model is originally established to describe the two-phase flow in pipes [36].It introduced two important parameters: the multiplier of water phase Φ w and the multiplier of gas phase Φ g , in order to assess the flow resistivity induced by two-phase interactions.They are defined as equations ( 2) and ( 3), respectively: in which dP/dx refers to the two-phase pressure gradient of both fluids, ðdP/dxÞ w refers to the pressure gradient in the virtual single-phase flow of water with a same flow rate as the water flow rate in the two-phase flow, and ðdP/dxÞ g refers to that of gas.The Martinelli parameter χ is used to analyze the evolution of the phase multipliers.It is defined as equation ( 4), which represents the fractions of two fluids, namely, the relative importance of the flow of water to gas.If χ is 0, it indicates the single-phase flow of gas; when χ increases, it means the increase of the fraction of water and the decrease of the fraction of gas; when χ approaches infinity, it means that the fraction of gas approaches 0 and the two-phase flow approaches the single-phase flow of water.
It should be noted that Φ w and Φ g are originally defined as the square root of the ratio between the two-phase pressure gradient and the single-phase one.For the convenience of comparison with the definition of relative permeability, Fourar and Bories [37] redefined them as equations ( 2) and (3).In this paper, we followed this revised definition by Fourar and Bories.
It should be noted that the stratified flow has formed due to the difference of density between gas and water, meaning that gas and water flowed in separate layers.In this situation, the capillary pressure does not contribute to the total pressure drop, so the gas pressure drop and the water pressure drop in the two-phase flow are assumed to be equal to each other.
As shown in Figure 1, two pressure sensors are connected to the water injection tube and gas injection tube.The water injection tube and gas injection tube are both connected to the inlet (the black parts of the 3D model).Theoretically, they should be identical but there are small differences between the gas pressure and the water pressure at the inlet.This may be induced by the compressibility of gas, which led to severe fluctuation of the pressure.The adopted inlet pressure for calculation is the average of the inlet pressure of water and that of gas.Both water and gas flows out from the same outlet, and one pressure sensor is connected to the outlet and its average value is adopted as the pressure of water and gas at the outlet.The two-phase pressure drop dP is calculated as the difference between the pressure at the inlet and the pressure at the outlet.
With the Lockhart-Martinelli model, the two-phase pressure field can be calculated by the following steps in engineering applications (assuming that in a two-phase flow case, the flow rate of gas Q g and the flow rate of water Q w are given): (1) Based on Forchheimer's law, the corresponding pressure gradient in the virtual single-phase flow of gas (which has the same flow rate as the gas flow rate in the two-phase flow) can be calculated as  4 Geofluids (2) With the same method, the corresponding pressure gradient in the virtual single-phase flow of water (which has the same flow rate as the water flow rate in the two-phase flow) can be calculated: The Lockhart-Martinelli model evaluates the two-phase flow by tracing the evolution of Φ w and Φ g with respect to χ.Since at the same flow rate, the pressure drop of each phase in the two-phase flow is always larger than that in the single-phase flow; these two parameters represent the extra resistivity induced by the two-phase interference.Actually, the Martinelli-Lockhart model is similar to the generalized Darcy's law, because 1/Φ w is similar to k rw (the relative permeability of water in generalized Darcy's law): 1/Φ w can be expressed as ðdP/dxÞ w /ðdP/dxÞ = ½ðdP/dxÞ w / Q w /½ðdP/dxÞ/Q w = ðthe two − phase conductivity of waterÞ/ ðthe single − phase conductivity of waterÞ, while k rw is defined as the ratio of effect permeability of water in the two-phase flow and the permeability in the single-phase flow.The difference lies in that the relative permeability does not consider the nonlinearity of flow but the Martinelli-Lockhart model does (the pressure gradients of the single-phase flow ðdP/dxÞ w and ðdP/dxÞ g are calculated by Forchheimer's law).Consequently, the Lockhart-Martinelli model can be seen as a generalization of the two-phase Darcy's law in order to make it fit the non-Darcy two-phase flow [37].
The Lockhart-Martinelli model and Delhaye's model are initially proposed for describing the two-phase flow in pipes.Since these two models consider the flow nonlinearity, they have the potential to be effective for describing the twophase flow in intersecting fractures.In this study, it is investigated whether the experiment data of the two-phase flow in the intersecting fractures can fit these models by the following steps: If the calculated Φ w , Φ g , and χ prove to be in accordance with Delhaye's model, the pressure field of the two-phase flow in intersecting fractures can be calculated with the five steps mentioned above in engineering applications.

Hydraulic Characteristics of the Two-Phase
Flow in the Intersecting Fractures 6 Geofluids and pressure drop ΔP is shown in Figure 4.Here the pressure drop ΔP is the pressure difference between the inlet and outlet.The results show that the pressure drop increases nonlinearly with respect to the flow rate, which is in accordance with the Forchheimer's law, which means that the flow nonlinearity cannot be neglected.Based on the results of single-phase flow test, the two-phase flow hydraulic characteristics can be calculated, as introduced in the next section.5 and 6 show the evolution of the two-phase pressure drop with respect to the gas flow rate of case 1 and case 2, respectively.Each curve represents the tests with a constant water flow rate.The pressure drop refers to the pressure difference between the inlet and outlet.Case 1 is composed of fractures 2, 5, and 3, with an aperture of 0.3 mm-1 mm-0.3 mm; case 2 is composed of fractures 1, 5, and 3, with an aperture of 0.6 mm-1 mm-0.3 mm.Case 1 has a smaller hydraulic aperture, and consequently, the pressure drop of case 1 is larger than that of case 2.

Hydraulic
In the two-phase flow in pipes, there are usually sharp variations in the P − Q (pressure drop-flow rate) curves if the flow structure changes.In our test, no sharp variations are observed in Figures 5 and 6.The reason can be sought from the flow structures.In our tests, the fractures in the 3D model were vertical to the horizontal plane, so gas and water flowed in separate layers and gas flowed above the water due to gravity, consequently forming a stratified flow.In the stratified flow, there is not an obvious transfer of flow 7 Geofluids structures from the bubble flow to the continuous flow, so there are no sharp variations in the P − Q curves.
As is indicated earlier, the single-phase flow in this 3D model used in our experiment is nonlinear and can be described by Forchheimer's law.As for the two-phase flow, the generalized Darcy's law does not account for the nonlinearity, either.Due to the existence of the fracture intersections in this 3D fracture model, which will contribute to flow nonlinearity, the generalized Darcy's law is not expected to fit the pressure drop data very well.However, we still analyzed the flow structures and measured the saturation as a reference.
Figure 7 shows the flow structures of case 1 at W r = 600 mL/min (water flow rate) and G r = 1500 mL/min (gas flow rate).The photos were taken from the lateral direction, as indicated by the red arrow in Figure 2. In the photos, the dark color represents water and the light color represents gas.It shows that the saturation of gas was fluctuating with respect to time, which is similar to the stratified wavy flow in pipes.Since fractures 1-4 were not vertical to the camera, to avoid refraction in taking photos, only the flow structures in fracture 5 are analyzed and the saturations are calculated.
The photos of fracture 5 are extracted from the original photos, as shown in Figures 8 and 9. Figure 8 shows the flow structures at W r = 500 mL/min and G r = 500 mL/min in case 1, while Figure 9 shows the flow structures at W r = 500 mL /min and G r = 1000 mL/min in case 1.The red pictures are the original photos taken by the camera, and the blackwhite pictures were converted from the original photos by binarization in order to calculate the water saturation S w .It shows that the saturation was fluctuating with respect to time, though the flow rates of water and gas were kept constant.When the gas flow rate increased from 500 mL/ min to 1000 mL/min, there was not an obvious variation in the flow structures.
Figure 10 shows the evolution of water saturation in case 2. Each saturation value is the average of ten pictures.It shows that the water saturation tends to decrease with respect to the gas flow rate but some exceptions exist, for example, in the curve of W r = 400 mL/min, when G r increases from 250 mL/min to 500 mL/min, the saturation increases.In addition, the variation of water saturation is not obvious.When W r is 700 mL/min, G r increases from 250 mL/min to 1500 mL/min, the water saturation decreased only by 10%.Besides, the curves at different water flow rates may intersect with each other.This also indicates that the generalized Darcy's law may not work well in describing the two-phase flow pressure drop in the intersecting fractures since the saturation cannot indicate the fractions of the two fluids very well.

Interpretation Model.
Since the nonlinearity of the flow is confirmed in Figure 4, the Lockhart-Martinelli model is adopted to analyze the hydraulic characteristics of the twophase flow in this experiment.The evolution of 1/Φ w and 1/Φ g in the two cases is shown in Figures 11 and 12. χð1 + χÞ is selected as the x-axis instead of χ for convenience, since in this way, the x-axis varies from 0 to unity.As shown in Figures 11 and 12, the sum of 1/Φ w and 1/Φ g are less than 1, indicating a quite significant interference between two phases.When χ/ð1 + χÞ approaches 1, it means that χ approaches infinity and the fraction of gas approaches 0, meaning that the two-phase flow approaches the singlephase flow of water.Consequently, in this condition, 1/Φ w approaches 1 and 1/Φ g approaches 0. All the experiment data approximately fall on the same curve, indicating that 1/Φ w and 1/Φ g are independent of the flow rates to some extent, but dependent on χ.
It should be noted that in Figures 11 and 12, 1/Φ g seems to be decreasing linearly and 1/Φ w is increasing nonlinearly.This is induced by the physical property difference between two phases.1/Φ w and 1/Φ g are defined as ðdP/dxÞ w /ðdP/d xÞ and ðdP/dxÞ g /ðdP/dxÞ, respectively, in which the twophase pressure gradients (dP/dx) are the same for either gas or water, but the single-phase pressure gradient of ðdP/dxÞ g is much less than that of ðdP/dxÞ w due to the difference of density.Therefore, the value range of 1/Φ g is also very small (within 0-0.1) and the nonlinearity cannot be well indicated in Figures 11 and 12.In Figures 13 and 14, the nonlinearity of both 1/Φ w and 1/Φ g can be well indicated.
As mentioned in the previous section, in order to predict the pressure field of the two-phase flow, an appropriate model between Φ w , Φ g , and χ is necessary.Delhaye et al. [38] proposed an empirical relationship, as shown in equations ( 5) and ( 6): in which C is a parameter that indicates the flow regime of each phase.When C is 5, both water and gas flow are laminar; when C is 10, water is turbulent and gas is laminar; when C is 12, liquid is laminar and gas is turbulent; when C is 20, both water and gas flow are turbulent.However, this method is initially proposed for describing the two-phase flow in pipes.Fourar and Bories [37] used this method to 9 Geofluids fit two-phase flow data in single fractures, and they found that this model can work well for approximate calculations.Since the two-phase flow in intersecting fractures has similarity to that in the rough single fracture in the aspect that the flow has nonlinearity induced by the inertial force, it is believed that this model may be also effective for the intersecting fractures.Therefore, the experiment data are fit with these two equations.
Figures 13 and 14 show the evolution of Φ w and Φ g with respect to Sqrt ðχÞ in case 1 and the corresponding fitting values of C. The curves of equations ( 5) and ( 6) are also plotted for reference.The five experiment curves approximately fall on the same curve, and they are close to C = 5, namely, the water laminar-gas laminar flow state.The flow rate has some weak impact on the values of C. To be accurate, when the water flow rates are 300 mL/min, 400 mL/min, 500 mL/  14 are also close to C = 5, which is consistent with the fitting values of C with equation ( 5) to some extent.12 Geofluids 500 mL/min, 600 mL/min, and 700 mL/min, the fitted values of C with equation ( 5) are 1.98, 2.30, 2.57, 2.59, and 2.61; the fitted values of C with equation ( 6) are 1.76, 1.81, 1.70, 1.55, and 1.05.Both groups of fitted values of C indicated that the turbulence in case 2 is smaller than that in case 1.This is reasonable since case 1 is composed of three fractures with an aperture of 0.3 mm-1 mm-0.3 mm, while case 2 is composed of three fractures with an aperture of 0.6 mm-1 mm-0.3 mm.The average fracture aperture of case 2 is larger than that of case 1, and in case 2, the superficial velocity is smaller, and consequently, the turbulence is small.In addition, the aperture change at the first intersection of case 2 (from 0.6 mm to 1 mm) is smaller than that of case 1 (from 0.3 mm to 1 mm), so in case 1, the turbulence induced by the aperture change at the fracture intersection is also serious.
To summarize, further studies should be conducted in the following aspects: (1) the values of C, which indicate different flow regimes, should be modified for the two-phase flow in intersecting fractures, to be different from those of the two-phase flow in pipes.(2) The relationship between C and the fracture configurations, such as the intersecting angle and the aperture variation, should be established.(3) The dependency of C on the fluid flow rate should be clarified.(4) In Delhaye's study, the values of C in equation (5) and C in equation ( 6) should be identical, but in our test, there is a little difference between the values of C calculated by these two equations.This means that equations ( 5) and ( 6) should be modified such as by multiplying a factor.

Summary
In this paper, a two-phase flow experiment is conducted in a 3D fracture model with smooth intersecting fractures to investigate its hydraulic characteristics.The following conclusions are supported by the results: (3) Delhaye's model can be used for approximate calculations of the two-phase pressure drop in the intersecting fractures, but not completely suitable.It should be modified for better fit.The phase multipliers Φ w and Φ g are believed to be related to the configuration of the intersecting fractures, such as the intersecting angle and the aperture change of the fractures.The range of values for C in Delhaye's model also needs to be reset.These aspects should be further investigated in future studies

5 F r a c t u r e 2 F r a c t u r e 1 F r a c t u r e 2 F 5 F r a c t u r e 4 F r a c t u r e 3 F r a c t u r e 4 FFigure 3 :
Figure 3: Schematic of the 3D fracture model (vertical view).

Figure 4 :
Figure 4: Test results of the single-phase flow.

( 1 )Figure 5 :
Figure 5: Evolution of the two-phase pressure drop with respect to the water flow rate and gas flow rate in case 1.

Figure 8 :
Figure 8: Saturations in fracture 5 at different moments when W r = 500 mL/min and G r = 500 mL/min in case 1.

Figure 9 :
Figure 9: Saturations in fracture 5 at different moments when Wr = 500 mL/min and Gr = 1000 mL/min in case 1.

Figure 10 :
Figure 10: The evolution of water saturation.

Figure 13 :
Figure 13: The relationship between Φ w and Sqrt ðχÞ in case 1.

Figure 15 :
Figure 15: The relationship between Φ w and Sqrt ðχÞ in case 2.

( 1 ) 2 )
The flow structures in the intersecting fractures show more similarity to those of the stratified wavy flow in pipes.The saturation kept fluctuating with respect to time, and it cannot reflect the fraction of two fluids very well.In addition, the nonlinearity induced by inertial force and turbulence in the intersecting fractures cannot be neglected in the two-phase flow.Due to these reasons, the generalized Darcy's law cannot well describe the hydraulic characteristics of the intersecting fractures; but the Lockhart-Martinelli model can fit the experiment data (For an approximate calculation of the two-phase flow pressure drop in the intersecting fractures, Delhaye's model is effective since the experimental data approximately tend to fall on the same curve; but C increases with respect to the water flow rate to some extent, indicating an increase of flow turbulence.By comparing different cases, the degrees of turbulence in different cases can be well indicated by the average values of the parameter C. Due to the significant aperture change in case 1, it showed more significant turbulence by the value of C Figure 12: The relationship between Φ w , Φ g , and χð1 + χÞ in case 2.Figure 11: The relationship between Φ w , Φ g , and χð1 + χÞ in case 1. 10 Geofluids min, 600 mL/min, and 700 mL/min, the corresponding C values are 4.05, 4.81, 5.30, 5.32, and 5.66, respectively.With the increase of the water flow rate, the value of C increases, which indicates an increase of turbulence of the flow.By fitting the test data to equation (6), the corresponding C values are 3.93, 3.62, 3.34, 2.96, and 2.78.The fitted values of C in Figure