Pore-Scale Simulation of Gas and Water Two-Phase Flow in Rough-Walled Fractures Using the Volume of Fluid Method

: The gas and water ﬂow behavior in rough-walled hydrophilic fractures at the pore scale is crucial for understanding the gas production characteristics of naturally fractured formations. This paper presents a systematic analysis of the gas and water ﬂow characteristics in both the single-fracture and Y-shaped junction fracture models using the volume of ﬂuid (VOF) method. Numerical simulations showed that the gas/water rate ratio is the most signiﬁcant factor inﬂuencing gas bubble/slug geometry, phase distribution, and saturation. The effect of fracture roughness and tortuosity is less signiﬁcant than the gas/water ratio, whereas the total ﬂuid rate has a negligible effect. For Y-shaped junction models, the phase distribution and referential pathways are predominantly controlled only by the channel aperture ratio, whereas the effect of the intersecting angle and ﬂuid ﬂow rate can be neglected.


Introduction
Natural gas from naturally fractured gas formations such as coal, shale, and carbonate is an important addition to global energy supplies [1][2][3].Gas and water two-phase flows occur frequently in discrete fractures due to the presence of either inherent water or injected water-based fracturing fluids in these formations [4][5][6].It is well recognized that the interactions of gas and water at pore scales exert profound effects on the spatial distribution and flow capacities of the fluids and thus the production characteristics [7][8][9].
To date, two categories of investigation tools, namely, laboratory microfluidics and numerical simulation, have been proposed to investigate the gas and water interaction behaviors at the pore scale.Microfluidics provides a platform to directly visualize multiphase flows at micron-scale geometries [10].A number of experiments have been conducted to investigate the gas and water flow patterns in a single fracture or complex fracture networks.Mansour et al. [11] investigated the gas and water flow pattern in a singlefracture micro-chip model.Wang et al. [12] performed gas and water flow experiments on a microchip model representing carbonate reservoirs and proposed several gas trapping patterns.Chen et al. [13] investigated the relative permeability of gas and water in a Yjunction fracture model subject to varying fluid flow rates.Gerami et al. [14] performed gas and water flow tests on both fork-shaped fracture networks and coal cleat chip models.The previous experimental work gave direct observations of how gas and water interact at the pore scale.However, most of these previous works used microchip models with specific wall roughness and thus did not provide insights into how the fracture roughness affects the flow pattern.Another issue with the microfluidic method is the interpretation accuracy due to the limitation of image resolution.The image resolution of current commercial microscopes is generally on the magnitude of micrometers.However, if the wall roughness Energies 2022, 15, 9382 2 of 15 is on the magnitude of nanometers, the recorded images are incapable of capturing the accurate fluid distribution on the rough surface.Besides, microfluidic experiments are usually quite expensive and require the integration of specially designed microfluidic chips, microscope(s), and high-speed digital camera(s).
Compared with the microfluidic experiment, numerical simulations are generally less expensive, more versatile, and provide rapid predictions for a large range of flow conditions [15,16].To date, various methods have been proposed for simulating multiphase flow at the pore scale (or mesoscale), among which, pore network modeling (PNM) and computational fluid dynamics (CFD) are the most widely used.The PNM performs simulations on conceptualized models consisting of simplified geometries of pores and throats rather than on the true geometry of the flow path [17][18][19].Thus, the accuracy of the PNM depends heavily on the accuracy of the reconstructed discrete physical pore-throat model.It is recognized that the PNM is incapable of producing reliable predictions of fluid-solid interactions and flow characteristics for geometric models with high irregularity [20].Moreover, the PNM requires a predefinition of certain mechanisms such as the snap-off and pore-filling [21,22] for modeling multiphase flows, which inevitably brings in additional errors.
Compared with the PNM method, the CFD method is capable of simulating fluid flows directly on true physical models with complex geometries and thus is considered to be more accurate than the PNM [23].Moreover, the CFD method is capable of modeling a more complex multiphase flow that involves mix-wettability [24], heat transfer [25], and solid particle flow [26].Among the various CFD methods, the LBM and VOF are the most commonly used for modeling two-phase flows in porous media at the mesoscale.The LBM method is superior to the VOF method in terms of convergence stability and conservation accuracy for simulating incompressible fluid flow [27,28].However, conventional LBM is associated with the issue of numerical instability when simulating multiphase flow with a high density and viscosity ratio [29] (e.g., gas/water or gas/oil two-phase flow).Moreover, the LBM method has inherent difficulties in handling compressibility problems when compressible fluid(s) is(are) considered.The VOF method has an inherent advantage of mass conservation [30] and does not require a complex phase interface tracking algorithm that is important for calculating two-phase flow in complex geometric shapes [26].Thus, the VOF method is more mature and reliable than the LBM for simulating multiphase flow with compressible fluid(s) and a large density/viscosity ratio in pathways with complex geometries.The underlying philosophy of VOF is to couple the Navier-Stokes equations with the volume fraction model and then solve the discretized equations with numerical computation techniques [31][32][33][34][35].
Although the VOF method has been proven to be capable of simulating compressible fluid flows in channels in meso-and macroscale sand packs or pipelines, very limited research has reported the use of the VOF method to investigate gas and water interactions in pore-scale fracture networks.To the best knowledge of the authors, the only paper that reported gas and water flow simulation in microscale fractures was by Huang et al. [36], who showed the accuracy of the CFD simulation for gas and water flow patterns in Yjunction models.However, all the simulations were based on smooth fractures, which do not conform to real rough fractures.Moreover, the model used by Huang et al. neglects the gas compressibility effects.This paper used the compressible VOF method to simulate two-phase gas and water flows in rough-walled single-fracture and Y-shaped junction models.The interactions between the two-phase fluids in different fracture geometries will be analyzed, and new insights will be discussed.

Numerical Simulation Method for Two-Phase Flow
The CFD methods used for simulating fluid flow have been extensively presented in previous papers (e.g., [37][38][39]).For the completeness of this paper, the basic mathematical models will be briefly addressed as follows.The fluid flow considering fluid compressibility is written as: where ρ, U, and P are phase density, velocity, and pressure, respectively.τ and F represent the sheer stress and surface tension, respectively.For a two-phase flow problem, it is required that the fluid interface be captured with certain techniques such as the VOF and the Level-set (LS) methods.It is commonly recognized that the VOF is superior to the LS method in terms of mass conservation and the simulation accuracy in complex structures [40,41].The basic philosophy behind the VOF method is to solve two single-phase flow problems independently outside the interfacial region and then couple the two-phase fluids in the interfacial region with specific assumptions and simplifications.The concept of volume fraction (α) is used to define the interfacial and non-interfacial zones, which is given as: where V i is the volume of the phase i, and V f is the control volume.For a two-phase gas and water flow problem, we can assume that α = 0 and α = 1 represent the gaseous and water regions, respectively.For the interface region, the volume fraction falls in the range of (0, 1), and the fluid properties such as the density and viscosity can be approximated as [42]: The VOF model describing the volume fraction is given as: Equations ( 1), ( 2) and ( 6) are coupled and iteratively solved using the finite volume method.
In this study, all simulations were conducted using the compressibleInterFoam solver implemented in the open-source CFD toolbox OpenFoam (https://openfoam.org/,accessed on 24 November 2022).The compressibleInterFoam solver is capable of simulating compressible and immiscible two-phase flow with a non-isothermal effect.The VOF method is implemented in the solver to capture the fluid interface.For models considered in this study, we used the SnappyHexMesh toolbox to construct the meshing.The strategy that combines coarse and fine meshing was used in order to ensure accuracy while maintaining high computational accuracy.The coarse meshing and fine meshing were applied in the center and boundary regions, respectively.The maximum total number of grids was set to be 2 × 10 6 considering Cha et al.'s simulation results [19] as a reference.

Validity Test of the VOF Method
To validate the accuracy of the VOF for simulating gas and water flow at the microscale, simulations were conducted and then compared with experimental results by Gerami et al. [14] and Mansour et al. [11] in this paper.
Gerami et al. [14] conducted the experiment of water charging into hydrophobic fork-shaped flow channels manufactured on polydimethylsiloxane (PDMS).The geometry of the flow channels is illustrated in Figure 1a.The flow channels were initially saturated with air.Deionized water was subsequently injected into the micro-channel at a constant flow rate of 1 µL/min.The distribution of water and air was recorded as the water invaded gradually into the flow channels.The experimental process was simulated in this paper, and the simulation results on gas/water distribution were compared with the experimental observations.
To validate the accuracy of the VOF for simulating gas and water flow at the microscale, simulations were conducted and then compared with experimental results by Gerami et al. [14] and Mansour et al. [11] in this paper.
Gerami et al. [14] conducted the experiment of water charging into hydrophobic forkshaped flow channels manufactured on polydimethylsiloxane (PDMS).The geometry of the flow channels is illustrated in Figure 1a.The flow channels were initially saturated with air.Deionized water was subsequently injected into the micro-channel at a constant flow rate of 1 μL/min.The distribution of water and air was recorded as the water invaded gradually into the flow channels.The experimental process was simulated in this paper, and the simulation results on gas/water distribution were compared with the experimental observations.Mansour et al. [11] performed experiments on the simultaneous flow of water and gas (nitrogen) through a T-junction channel.The flow channel was manufactured using the UV lithography technique on PDMS material.The wall surface of the flow channel was modified using the oxygen plasma treatment to alter the wettability from hydrophobicity to hydrophilicity.The geometry of the T-junction channel used for experiments is illustrated in Figure 1b.Simulation scenarios with different gas/water ratios are to be conducted to identify the flow pattern, which will then be compared with the experimental results.

Simulation of Gas/Water Flow in Rough-Walled Fracture Models
Previous studies [43][44][45] and our preliminary tests showed that natural fractures are generally non-smooth and exhibit a certain degree of tortuosity (Figure 2a), which may affect fluid flow therein.For naturally fractured reservoirs, the geometry of the fracture system is complex and typically composed of a number of inter-connected single fractures with relatively large ranges of apertures [46,47].Our preliminary tests confirmed that the basic joints of the intersected fractures are generally Y-shaped or T-shaped junctions (Figure 2b).As a matter of fact, the T-shape junctions can be considered as specific cases of Yshaped junctions with a right intersecting angle.Mansour et al. [11] performed experiments on the simultaneous flow of water and gas (nitrogen) through a T-junction channel.The flow channel was manufactured using the UV lithography technique on PDMS material.The wall surface of the flow channel was modified using the oxygen plasma treatment to alter the wettability from hydrophobicity to hydrophilicity.The geometry of the T-junction channel used for experiments is illustrated in Figure 1b.Simulation scenarios with different gas/water ratios are to be conducted to identify the flow pattern, which will then be compared with the experimental results.

Simulation of Gas/Water Flow in Rough-Walled Fracture Models
Previous studies [43][44][45] and our preliminary tests showed that natural fractures are generally non-smooth and exhibit a certain degree of tortuosity (Figure 2a), which may affect fluid flow therein.For naturally fractured reservoirs, the geometry of the fracture system is complex and typically composed of a number of inter-connected single fractures with relatively large ranges of apertures [46,47].Our preliminary tests confirmed that the basic joints of the intersected fractures are generally Y-shaped or T-shaped junctions (Figure 2b).As a matter of fact, the T-shape junctions can be considered as specific cases of Y-shaped junctions with a right intersecting angle.
To validate the accuracy of the VOF for simulating gas and water flow at the microscale, simulations were conducted and then compared with experimental results by Gerami et al. [14] and Mansour et al. [11] in this paper.
Gerami et al. [14] conducted the experiment of water charging into hydrophobic forkshaped flow channels manufactured on polydimethylsiloxane (PDMS).The geometry of the flow channels is illustrated in Figure 1a.The flow channels were initially saturated with air.Deionized water was subsequently injected into the micro-channel at a constant flow rate of 1 μL/min.The distribution of water and air was recorded as the water invaded gradually into the flow channels.The experimental process was simulated in this paper, and the simulation results on gas/water distribution were compared with the experimental observations.Mansour et al. [11] performed experiments on the simultaneous flow of water and gas (nitrogen) through a T-junction channel.The flow channel was manufactured using the UV lithography technique on PDMS material.The wall surface of the flow channel was modified using the oxygen plasma treatment to alter the wettability from hydrophobicity to hydrophilicity.The geometry of the T-junction channel used for experiments is illustrated in Figure 1b.Simulation scenarios with different gas/water ratios are to be conducted to identify the flow pattern, which will then be compared with the experimental results.

Simulation of Gas/Water Flow in Rough-Walled Fracture Models
Previous studies [43][44][45] and our preliminary tests showed that natural fractures are generally non-smooth and exhibit a certain degree of tortuosity (Figure 2a), which may affect fluid flow therein.For naturally fractured reservoirs, the geometry of the fracture system is complex and typically composed of a number of inter-connected single fractures with relatively large ranges of apertures [46,47].Our preliminary tests confirmed that the basic joints of the intersected fractures are generally Y-shaped or T-shaped junctions (Figure 2b).As a matter of fact, the T-shape junctions can be considered as specific cases of Yshaped junctions with a right intersecting angle.

Single-Fracture Model
The geometry of a basic single-fracture model is shown in Figure 3.The fracture length and mean aperture were set to be 4500 and 300 µm, respectively.The roughness of a fracture has been demonstrated previously to exert profound effects on fluid flow in microchannels [48,49].It is well-recognized that almost all natural and hydraulic fractures are with a certain degree of roughness.In this study, the roughness is set by adding a random noise on the aperture along the fracture, which is quantitatively represented as the Energies 2022, 15, 9382 5 of 15 averaged deviations of the fracture wall inner surface from the average fracture aperture (Figure 3).The mathematical expression for calculating the roughness is: where L is the length of a fracture; Z(x) is the distance between the inner surface and the averaged center of the fracture at point x; Z 0 is the half of the average fracture aperture, which is quantitatively calculated as the average distance between the inner surface and the center of the of the fracture.

Single-Fracture Model
The geometry of a basic single-fracture model is shown in Figure 3.The fracture length and mean aperture were set to be 4500 and 300 μm, respectively.The roughness of a fracture has been demonstrated previously to exert profound effects on fluid flow in microchannels [48,49].It is well-recognized that almost all natural and hydraulic fractures are with a certain degree of roughness.In this study, the roughness is set by adding a random noise on the aperture along the fracture, which is quantitatively represented as the averaged deviations of the fracture wall inner surface from the average fracture aperture (Figure 3).The mathematical expression for calculating the roughness is: where L is the length of a fracture; Z(x) is the distance between the inner surface and the averaged center of the fracture at point x; Z0 is the half of the average fracture aperture, which is quantitatively calculated as the average distance between the inner surface and the center of the of the fracture.The tortuosity is defined as the ratio of the length of the true flow path to the geometrical length of the flow channel (Figure 3).Generally, a higher tortuosity results in an increased flow path length and thus increases the flow resistance [50,51].To investigate the effect of the tortuosity on the two-phase flow, single-fracture models with varying tortuosity were synthetized.The synthesis of the torturous rough fractures includes three steps.First, two sine functions with different frequencies were used to generate the base trend lines of the fracture's upper and lower wall boundaries, respectively.Second, the trend lines were shifted with a random frequency to form the wall boundary of the fracture.Third, both the upper and lower boundaries were added randomly distributed noises to represent rough walls.The true flow path was approximated as the middle points between the upper and lower boundaries of the fracture surface.To calculate the tortuosity, the true flow and the geometrical length were obtained using computer-aided design software.
To investigate the gas and water flow characteristics in single-fracture models, simulations were conducted considering the effect of roughness, tortuosity, flow rate, and gas/water ratio.The simulation scenarios are shown in Table 1.Since this paper is primarily concerned with low-velocity fluid flow in fractured gas reservoirs, relatively small flow rate values were considered to mimic the underground conditions.Previous studies [52,53] showed that most naturally fractured formations such as coal and carbonate are water-wet.However, the contact angle may vary significantly with mineral content, roughness, pressure, etc.Therefore, the inner wall surface was assumed to have a contact angle of 45°, which is the median between the complete water-wetting condition (with a contact angle of 0°) and the cutting-off value of the gas-wetting condition (with a contact angle of 90°).To simulate the "simultaneous" injection of gas and water, a T-junction representing two inlets was used (Figure 3).It is noted that the T-junction was assumed to The tortuosity is defined as the ratio of the length of the true flow path to the geometrical length of the flow channel (Figure 3).Generally, a higher tortuosity results in an increased flow path length and thus increases the flow resistance [50,51].To investigate the effect of the tortuosity on the two-phase flow, single-fracture models with varying tortuosity were synthetized.The synthesis of the torturous rough fractures includes three steps.First, two sine functions with different frequencies were used to generate the base trend lines of the fracture's upper and lower wall boundaries, respectively.Second, the trend lines were shifted with a random frequency to form the wall boundary of the fracture.Third, both the upper and lower boundaries were added randomly distributed noises to represent rough walls.The true flow path was approximated as the middle points between the upper and lower boundaries of the fracture surface.To calculate the tortuosity, the true flow and the geometrical length were obtained using computer-aided design software.
To investigate the gas and water flow characteristics in single-fracture models, simulations were conducted considering the effect of roughness, tortuosity, flow rate, and gas/water ratio.The simulation scenarios are shown in Table 1.Since this paper is primarily concerned with low-velocity fluid flow in fractured gas reservoirs, relatively small flow rate values were considered to mimic the underground conditions.Previous studies [52,53] showed that most naturally fractured formations such as coal and carbonate are water-wet.However, the contact angle may vary significantly with mineral content, roughness, pressure, etc.Therefore, the inner wall surface was assumed to have a contact angle of 45 • , which is the median between the complete water-wetting condition (with a contact angle of 0 • ) and the cutting-off value of the gas-wetting condition (with a contact angle of 90 • ).To simulate the "simultaneous" injection of gas and water, a T-junction representing two inlets was used (Figure 3).It is noted that the T-junction was assumed to have a smooth inner surface in order to ensure a complete mix of the two phases once injected into the fracture model.The water and gas flow into the fracture model at their respective constant flow rates through the horizontal and vertical inlets, respectively.The outlet was constrained with a constant pressure of 0.1 MPa.Different scenarios of Y-shaped fracture junction models were constructed to investigate the effect of Y-shaped junction geometry, flow rate, and gas/water ratio on the phase distribution.The basic Y-shaped junction model is shown in Figure 4.The roughness and gas/water ratio for the base model were assigned to be 0.4 and 1, respectively.Varying intersecting angles (30 • , 60 • , and 90 • ) and fracture aperture ratios were considered to investigate the fracture geometry on the phase distribution.The apertures of the inlet and the upper outlet remained fixed at 300 µm and 200 µm, respectively.The aperture of the lower outlet was varied at 200, 300, and 400 µm to examine the effect of the fracture aperture ratio.have a smooth inner surface in order to ensure a complete mix of the two phases once injected into the fracture model.The water and gas flow into the fracture model at their respective constant flow rates through the horizontal and vertical inlets, respectively.The outlet was constrained with a constant pressure of 0.1 MPa.Different scenarios of Y-shaped fracture junction models were constructed to investigate the effect of Y-shaped junction geometry, flow rate, and gas/water ratio on the phase distribution.The basic Y-shaped junction model is shown in Figure 4.The roughness and gas/water ratio for the base model were assigned to be 0.4 and 1, respectively.Varying intersecting angles (30°, 60°, and 90°) and fracture aperture ratios were considered to investigate the fracture geometry on the phase distribution.The apertures of the inlet and the upper outlet remained fixed at 300 μm and 200 μm, respectively.The aperture of the lower outlet was varied at 200, 300, and 400 μm to examine the effect of the fracture aperture ratio.

Validation of the VOF Method
Figure 5 depicts the comparison between the VOF method and Gerami et al.'s experimental results of the fluid distribution in the fork-shaped gas-wet channel model during the water invasion process.It can be seen that the VOF method is capable of accurately reproducing the experiment results.Both the experiment and simulation showed that the injected water first invaded into the largest channel and displaced the gas phase-out.As the injection continued, the water phase gradually invaded into the middle channel with the moderate aperture, whereas the narrowest channel was never invaded throughout the test.The underlying mechanism is quite straightforward because the resistance force to water flow decreases with the channel width due to the hydrophobic nature of the model; therefore, it is easier for the water phase to invade the larger channel with less resistance.

Validation of the VOF Method
Figure 5 depicts the comparison between the VOF method and Gerami et al.'s experimental results of the fluid distribution in the fork-shaped gas-wet channel model during the water invasion process.It can be seen that the VOF method is capable of accurately reproducing the experiment results.Both the experiment and simulation showed that the injected water first invaded into the largest channel and displaced the gas phase-out.As the injection continued, the water phase gradually invaded into the middle channel with the moderate aperture, whereas the narrowest channel was never invaded throughout the test.The underlying mechanism is quite straightforward because the resistance force to water flow decreases with the channel width due to the hydrophobic nature of the model; therefore, it is easier for the water phase to invade the larger channel with less resistance.Figure 6 compares the simulation results with the experimental observations of the simultaneous flow of gas and water by Mansour et al. [8].It is shown that the simulation results agree well with the experiment in terms of the bubble patterns (Figure 6a). Figure 6b compares the experimental and simulated bubble aspect ratio (BAR).The BAR is defined as the ratio of the bubble length along the major axis to that along the minor axis.In this study, the major and minor axes correspond to the directions along and perpendicular to the flow direction, respectively.The BAR is a basic parameter describing the shape of a bubble and reflects the underlying interactions between two fluids.It can be seen that the simulated and experimental BARs are also in fine agreement.The results shown in Figures 5 and 6

Effect of Roughness
Figure 7 depicts the gas and water distribution in the single-fracture models with different roughness.It can be seen that all the gas bubbles have generally identical sizes, and the bubbles are distributed evenly along the smooth fracture (with zero roughness).As the roughness increases, the shape of the gas bubble becomes more irregular and the shape of each bubble depends heavily on the local surface roughness.Moreover, the distance between any two neighboring bubbles is also affected by the surface roughness.Figure 6 compares the simulation results with the experimental observations of the simultaneous flow of gas and water by Mansour et al. [8].It is shown that the simulation results agree well with the experiment in terms of the bubble patterns (Figure 6a). Figure 6b compares the experimental and simulated bubble aspect ratio (BAR).The BAR is defined as the ratio of the bubble length along the major axis to that along the minor axis.In this study, the major and minor axes correspond to the directions along and perpendicular to the flow direction, respectively.The BAR is a basic parameter describing the shape of a bubble and reflects the underlying interactions between two fluids.It can be seen that the simulated and experimental BARs are also in fine agreement.The results shown in Figures 5 and 6 confirm the validity of the VOF method for compressible two-phase flow in microscale channels.Figure 6 compares the simulation results with the experimental observations of the simultaneous flow of gas and water by Mansour et al. [8].It is shown that the simulation results agree well with the experiment in terms of the bubble patterns (Figure 6a). Figure 6b compares the experimental and simulated bubble aspect ratio (BAR).The BAR is defined as the ratio of the bubble length along the major axis to that along the minor axis.In this study, the major and minor axes correspond to the directions along and perpendicular to the flow direction, respectively.The BAR is a basic parameter describing the shape of a bubble and reflects the underlying interactions between two fluids.It can be seen that the simulated and experimental BARs are also in fine agreement.The results shown in Figures 5 and 6 confirm the validity of the VOF method for compressible two-phase flow in microscale channels.

Effect of Roughness
Figure 7 depicts the gas and water distribution in the single-fracture models with different roughness.It can be seen that all the gas bubbles have generally identical sizes, and the bubbles are distributed evenly along the smooth fracture (with zero roughness).As the roughness increases, the shape of the gas bubble becomes more irregular and the shape of each bubble depends heavily on the local surface roughness.Moreover, the distance between any two neighboring bubbles is also affected by the surface roughness.3.2.Gas/Water Flow Simulations in Rough-Walled Single-Fracture Models 3.2.1.Effect of Roughness Figure 7 depicts the gas and water distribution in the single-fracture models with different roughness.It can be seen that all the gas bubbles have generally identical sizes, and the bubbles are distributed evenly along the smooth fracture (with zero roughness).As the roughness increases, the shape of the gas bubble becomes more irregular and the shape of each bubble depends heavily on the local surface roughness.Moreover, the distance between any two neighboring bubbles is also affected by the surface roughness.Figure 8 shows the evolution of water saturation as simultaneous gas and water flow occurs along the rough fracture.As depicted, the water saturation exhibits a general linear decrease trend for each roughness value at the initial injection stage when gas is introduced into the fracture.As gas and water injection continues, the water saturation fluctuates around an average value for each roughness value.The fluctuation is due to the alternating introduction of water and gas into the fractures.Generally, the gas and water flow can be considered as "stable" after a duration of approximately 10 s for all the roughness cases.The residual water saturation is considered to be the average water saturation at the stable stage.Figure 9 shows the effect of roughness on BAR and residual water saturation.As shown in Figure 9a, the average BAR exhibits a slightly increasing trend with increasing roughness.Overall, the effect of roughness on the BAR can be neglected.As shown in Figure 9b, the residual water saturation climbs as the roughness increases.The increase in the residual water saturation is due to the increase in the end-corner structures due to increasing roughness, which provides increasing accommodation space for the residual water (see Figure 7d).Such observation is consistent with [54].   Figure 8 shows the evolution of water saturation as simultaneous gas and water flow occurs along the rough fracture.As depicted, the water saturation exhibits a general linear decrease trend for each roughness value at the initial injection stage when gas is introduced into the fracture.As gas and water injection continues, the water saturation fluctuates around an average value for each roughness value.The fluctuation is due to the alternating introduction of water and gas into the fractures.Generally, the gas and water flow can be considered as "stable" after a duration of approximately 10 s for all the roughness cases.The residual water saturation is considered to be the average water saturation at the stable stage.Figure 8 shows the evolution of water saturation as simultaneous gas and water flow occurs along the rough fracture.As depicted, the water saturation exhibits a general linear decrease trend for each roughness value at the initial injection stage when gas is introduced into the fracture.As gas and water injection continues, the water saturation fluctuates around an average value for each roughness value.The fluctuation is due to the alternating introduction of water and gas into the fractures.Generally, the gas and water flow can be considered as "stable" after a duration of approximately 10 s for all the roughness cases.The residual water saturation is considered to be the average water saturation at the stable stage.Figure 9 shows the effect of roughness on BAR and residual water saturation.As shown in Figure 9a, the average BAR exhibits a slightly increasing trend with increasing roughness.Overall, the effect of roughness on the BAR can be neglected.As shown in Figure 9b, the residual water saturation climbs as the roughness increases.The increase in the residual water saturation is due to the increase in the end-corner structures due to increasing roughness, which provides increasing accommodation space for the residual water (see Figure 7d).Such observation is consistent with [54].   Figure 9 shows the effect of roughness on BAR and residual water saturation.As shown in Figure 9a, the average BAR exhibits a slightly increasing trend with increasing roughness.Overall, the effect of roughness on the BAR can be neglected.As shown in Figure 9b, the residual water saturation climbs as the roughness increases.The increase in the residual water saturation is due to the increase in the end-corner structures due to increasing roughness, which provides increasing accommodation space for the residual water (see Figure 7d).Such observation is consistent with [54]. Figure 8 shows the evolution of water saturation as simultaneous gas and water flow occurs along the rough fracture.As depicted, the water saturation exhibits a general linear decrease trend for each roughness value at the initial injection stage when gas is introduced into the fracture.As gas and water injection continues, the water saturation fluctuates around an average value for each roughness value.The fluctuation is due to the alternating introduction of water and gas into the fractures.Generally, the gas and water flow can be considered as "stable" after a duration of approximately 10 s for all the roughness cases.The residual water saturation is considered to be the average water saturation at the stable stage.Figure 9 shows the effect of roughness on BAR and residual water saturation.As shown in Figure 9a, the average BAR exhibits a slightly increasing trend with increasing roughness.Overall, the effect of roughness on the BAR can be neglected.As shown in Figure 9b, the residual water saturation climbs as the roughness increases.The increase in the residual water saturation is due to the increase in the end-corner structures due to increasing roughness, which provides increasing accommodation space for the residual water (see Figure 7d).Such observation is consistent with [54].

Effect of the Tortuosity
Figure 10 depicts the effect of tortuosity on gas and water distribution.It can be seen that the gas bubble tends to be 'stretched' along the flow direction as the tortuosity increases.The stretch of the gas bubble can be quantitatively described with the BAR as shown in Figure 11a.It can be seen that the stretch effect becomes stable at high tortuosity.It is presumed that the reason why the BAR remains almost constant at a critical tortuosity value of 1.3 is due to the snap-off effect constrained by the tortuous fracture.However, further investigations are needed to validate this presumption and to examine the predominant factors influencing the critical tortuosity value.As depicted in Figure 11b, the residual water saturation tends to increase with the tortuosity, which is due to the increased possibility of water trapping occurring in the end-corners of the tortuous flow paths (Figure 10d).However, the percentage of the increase in the water saturation is limited (within 10%).

Effect of the Tortuosity
Figure 10 depicts the effect of tortuosity on gas and water distribution.It can be seen that the gas bubble tends to be 'stretched' along the flow direction as the tortuosity increases.The stretch of the gas bubble can be quantitatively described with the BAR as shown in Figure 11a.It can be seen that the stretch effect becomes stable at high tortuosity.It is presumed that the reason why the BAR remains almost constant at a critical tortuosity value of 1.3 is due to the snap-off effect constrained by the tortuous fracture.However, further investigations are needed to validate this presumption and to examine the predominant factors influencing the critical tortuosity value.As depicted in Figure 11b, the residual water saturation tends to increase with the tortuosity, which is due to the increased possibility of water trapping occurring in the end-corners of the tortuous flow paths (Figure 10d).However, the percentage of the increase in the water saturation is limited (within 10%).

Effect of Gas/Water Ratio
Figure 12 shows the gas and water phase distributions with different gas/water ratios.It is depicted that with the increase in the gas/water ratio, the gas phase expands gradually from a bubble pattern to a slug pattern.Moreover, the distance between two neighboring bubbles or slugs is decreased as the gas/water ratio increases.For the scenario of a high gas/water ratio of 9:1, it can be seen clearly that residual water is trapped in the wall corners due to the roughness nature.Comparing Figures 6, 8, and 10 leads to the conclusion that the gas/water ratio exerts a more significant effect on the gas and water distributions in the rough-walled fracture than roughness and tortuosity.

Effect of the Tortuosity
Figure 10 depicts the effect of tortuosity on gas and water distribution.It can be seen that the gas bubble tends to be 'stretched' along the flow direction as the tortuosity increases.The stretch of the gas bubble can be quantitatively described with the BAR as shown in Figure 11a.It can be seen that the stretch effect becomes stable at high tortuosity.It is presumed that the reason why the BAR remains almost constant at a critical tortuosity value of 1.3 is due to the snap-off effect constrained by the tortuous fracture.However, further investigations are needed to validate this presumption and to examine the predominant factors influencing the critical tortuosity value.As depicted in Figure 11b, the residual water saturation tends to increase with the tortuosity, which is due to the increased possibility of water trapping occurring in the end-corners of the tortuous flow paths (Figure 10d).However, the percentage of the increase in the water saturation is limited (within 10%).

Effect of Gas/Water Ratio
Figure 12 shows the gas and water phase distributions with different gas/water ratios.It is depicted that with the increase in the gas/water ratio, the gas phase expands gradually from a bubble pattern to a slug pattern.Moreover, the distance between two neighboring bubbles or slugs is decreased as the gas/water ratio increases.For the scenario of a high gas/water ratio of 9:1, it can be seen clearly that residual water is trapped in the wall corners due to the roughness nature.Comparing Figures 6, 8, and 10 leads to the conclusion that the gas/water ratio exerts a more significant effect on the gas and water distributions in the rough-walled fracture than roughness and tortuosity.

Effect of Gas/Water Ratio
Figure 12 shows the gas and water phase distributions with different gas/water ratios.It is depicted that with the increase in the gas/water ratio, the gas phase expands gradually from a bubble pattern to a slug pattern.Moreover, the distance between two neighboring bubbles or slugs is decreased as the gas/water ratio increases.For the scenario of a high gas/water ratio of 9:1, it can be seen clearly that residual water is trapped in the wall corners due to the roughness nature.Comparing Figures 6, 8 and 10 leads to the conclusion that the gas/water ratio exerts a more significant effect on the gas and water distributions in the rough-walled fracture than roughness and tortuosity.The qualitative analyses above are consistent with the quantitative correlations shown in Figure 13, which depicts significant changes in the BAR and residual water saturation with respect to the varying gas/water ratio.Overall, the BAR increases linearly with the gas/water ratio, whereas the residual water saturation exhibits a logarithmic decreasing trend with the gas/water ratio.

Effect of Total Flow Rate
Figure 14 depicts that the gas bubble geometry and phase distribution pattern exhibit minor changes when the total flow rate was varied from 0.2 to 8 × 10 −3 m/s.Our statistics show that the BARs and residual water saturations for these flow rate conditions are in relatively small ranges of 1.45 to 1.5, and 0.54 to 0.56, respectively.Chen et al.'s [10] microfluidic experiments concluded that gas and water flow patterns in the fracture junction can be affected by the flow rate, which is inconsistent with our observations.The inconsistency may be possibly associated with the inherent experimental errors due to the fluctuations of alternating gas and water flows.These statistical values suggest that the total flow rate exerts a negligible effect on the gas and water phase distribution and flow pattern.In other words, the so-called waterblockage effect can be neglected for a single fracture.These simulation results indicate that relatively large dewatering rates or large pressure differences can be applied in order to accelerate the gas recovery factor for reservoirs where single fractures are developed.The qualitative analyses above are consistent with the quantitative correlations shown in Figure 13, which depicts significant changes in the BAR and residual water saturation with respect to the varying gas/water ratio.Overall, the BAR increases linearly with the gas/water ratio, whereas the residual water saturation exhibits a logarithmic decreasing trend with the gas/water ratio.The qualitative analyses above are consistent with the quantitative correlations shown in Figure 13, which depicts significant changes in the BAR and residual water saturation with respect to the varying gas/water ratio.Overall, the BAR increases linearly with the gas/water ratio, whereas the residual water saturation exhibits a logarithmic decreasing trend with the gas/water ratio.

Effect of Total Flow Rate
Figure 14 depicts that the gas bubble geometry and phase distribution pattern exhibit minor changes when the total flow rate was varied from 0.2 to 8 × 10 −3 m/s.Our statistics show that the BARs and residual water saturations for these flow rate conditions are in relatively small ranges of 1.45 to 1.5, and 0.54 to 0.56, respectively.Chen et al.'s [10] microfluidic experiments concluded that gas and water flow patterns in the fracture junction can be affected by the flow rate, which is inconsistent with our observations.The inconsistency may be possibly associated with the inherent experimental errors due to the fluctuations of alternating gas and water flows.These statistical values suggest that the total flow rate exerts a negligible effect on the gas and water phase distribution and flow pattern.In other words, the so-called waterblockage effect can be neglected for a single fracture.These simulation results indicate that relatively large dewatering rates or large pressure differences can be applied in order to accelerate the gas recovery factor for reservoirs where single fractures are developed.

Effect of Total Flow Rate
Figure 14 depicts that the gas bubble geometry and phase distribution pattern exhibit minor changes when the total flow rate was varied from 0.2 to 8 × 10 −3 m/s.Our statistics show that the BARs and residual water saturations for these flow rate conditions are in relatively small ranges of 1.45 to 1.5, and 0.54 to 0.56, respectively.Chen et al.'s [10] microfluidic experiments concluded that gas and water flow patterns in the fracture junction can be affected by the flow rate, which is inconsistent with our observations.The inconsistency may be possibly associated with the inherent experimental errors due to the fluctuations of alternating gas and water flows.The qualitative analyses above are consistent with the quantitative correlations shown in Figure 13, which depicts significant changes in the BAR and residual water saturation with respect to the varying gas/water ratio.Overall, the BAR increases linearly with the gas/water ratio, whereas the residual water saturation exhibits a logarithmic decreasing trend with the gas/water ratio.

Effect of Total Flow Rate
Figure 14 depicts that the gas bubble geometry and phase distribution pattern exhibit minor changes when the total flow rate was varied from 0.2 to 8 × 10 −3 m/s.Our statistics show that the BARs and residual water saturations for these flow rate conditions are in relatively small ranges of 1.45 to 1.5, and 0.54 to 0.56, respectively.Chen et al.'s [10] microfluidic experiments concluded that gas and water flow patterns in the fracture junction can be affected by the flow rate, which is inconsistent with our observations.The inconsistency may be possibly associated with the inherent experimental errors due to the fluctuations of alternating gas and water flows.These statistical values suggest that the total flow rate exerts a negligible effect on the gas and water phase distribution and flow pattern.In other words, the so-called waterblockage effect can be neglected for a single fracture.These simulation results indicate that relatively large dewatering rates or large pressure differences can be applied in order to accelerate the gas recovery factor for reservoirs where single fractures are developed.These statistical values suggest that the total flow rate exerts a negligible effect on the gas and water phase distribution and flow pattern.In other words, the so-called waterblockage effect can be neglected for a single fracture.These simulation results indicate that relatively large dewatering rates or large pressure differences can be applied in order to accelerate the gas recovery factor for reservoirs where single fractures are developed.Figure 15 shows the water and gas distributions for Y-shaped junction models with varying aperture ratios.As can be seen from Figure 15b,c, the gas and water phases are generally distributed evenly along the flow channels where the downstream outlets have identical apertures of 200 µm.In this case, the flow paths for the water and gas phases are identical, viz., the two phases do not have their own preferential flow path.Figure 15 shows the water and gas distributions for Y-shaped junction models with varying aperture ratios.As can be seen from Figure 15b,c, the gas and water phases are generally distributed evenly along the flow channels where the downstream outlets have identical apertures of 200 μm.In this case, the flow paths for the water and gas phases are identical, viz., the two phases do not have their own preferential flow path.However, as the lower outlet has increased apertures of 300 and 400 μm (corresponding to an aperture ratio of 1.5 and 2.0), the gas phase flows out of the model only through the larger outlet channel, and no gas exists within the smaller outlet channel (Figure 15f,i).For these two models, the water phase flows through both the larger and the smaller outlet channels (Figure 16).In other words, preferential flow through the larger channel occurs for the gas phase.This observation is within expectations, considering the Young-Laplace equation that results in reduced resistance for the gas phase in a hydrophilic channel with a larger aperture.However, as the lower outlet has increased apertures of 300 and 400 µm (corresponding to an aperture ratio of 1.5 and 2.0), the gas phase flows out of the model only through the larger outlet channel, and no gas exists within the smaller outlet channel (Figure 15f,i).For these two models, the water phase flows through both the larger and the smaller outlet channels (Figure 16).In other words, preferential flow through the larger channel occurs for the gas phase.This observation is within expectations, considering the Young-Laplace equation that results in reduced resistance for the gas phase in a hydrophilic channel with a larger aperture.

Figure 1 .
Figure 1.Illustration of the geometries of the fracture models used for validating the simulation method.(a) Gerami model; (b) Mansour model.

Figure 2 .
Figure 2. Examples of SEM images of rough-walled tortuous single-fracture (a) and Y-shaped fracture junction (b) in coals.

Figure 1 .
Figure 1.Illustration of the geometries of the fracture models used for validating the simulation method.(a) Gerami model; (b) Mansour model.

Figure 1 .
Figure 1.Illustration of the geometries of the fracture models used for validating the simulation method.(a) Gerami model; (b) Mansour model.

Figure 2 .
Figure 2. Examples of SEM images of rough-walled tortuous single-fracture (a) and Y-shaped fracture junction (b) in coals.

Figure 2 .
Figure 2. Examples of SEM images of rough-walled tortuous single-fracture (a) and Y-shaped fracture junction (b) in coals.

Figure 3 .
Figure 3.The basic geometry of a rough-walled single-fracture model.

Figure 3 .
Figure 3.The basic geometry of a rough-walled single-fracture model.

Figure 4 .
Figure 4. Illustration of the geometry of the Y-shaped fracture model.

Figure 4 .
Figure 4. Illustration of the geometry of the Y-shaped fracture model.Similar to the constraint conditions of the single-fracture models, for each Y-shaped junction model, gas and water were simultaneously injected through a T-junction into the fracture model at given flow rates.The outlets of the model were assumed to operate at a constant pressure of 0.1 MPa.

Figure 5 .
Figure 5.Comparison of the simulation with experimental results for water invasion into forkshaped adjacent channel model at time steps of (a) 10 s; (b) 30 s; (c) 40 s.The upper and lower rows are experimental and simulation results, respectively.The water and gas phases are in orange and white, respectively, for the experiment results in the upper row.

Figure 6 .
Figure6compares the simulation results with the experimental observations of the simultaneous flow of gas and water by Mansour et al.[8].It is shown that the simulation results agree well with the experiment in terms of the bubble patterns (Figure6a).Figure6bcompares the experimental and simulated bubble aspect ratio (BAR).The BAR is defined as the ratio of the bubble length along the major axis to that along the minor axis.In this study, the major and minor axes correspond to the directions along and perpendicular to the flow direction, respectively.The BAR is a basic parameter describing the shape of a bubble and reflects the underlying interactions between two fluids.It can be seen that the simulated and experimental BARs are also in fine agreement.The results shown in Figures5 and 6confirm the validity of the VOF method for compressible two-phase flow in microscale channels.

Figure 5 .
Figure 5.Comparison of the simulation with experimental results for water invasion into fork-shaped adjacent channel model at time steps of (a) 10 s; (b) 30 s; (c) 40 s.The upper and lower rows are experimental and simulation results, respectively.The water and gas phases are in orange and white, respectively, for the experiment results in the upper row.

Figure 5 .
Figure 5.Comparison of the simulation with experimental results for water invasion into forkshaped adjacent channel model at time steps of (a) 10 s; (b) 30 s; (c) 40 s.The upper and lower rows are experimental and simulation results, respectively.The water and gas phases are in orange and white, respectively, for the experiment results in the upper row.

Figure 6 .
Figure 6.Comparison of the simulation with experimental results for the single-fracture model.(a) Comparison of water and gas distribution patterns, with the color bar representing water phase saturation; (b) comparison of the BAR.

Figure 6 .
Figure 6.Comparison of the simulation with experimental results for the single-fracture model.(a) Comparison of water and gas distribution patterns, with the color bar representing water phase saturation; (b) comparison of the BAR.

Figure 8 .
Figure 8. Variation of average water saturation with time for single-fracture models.

Figure 9 .
Figure 9.Effect of roughness on (a) bubble aspect ratio and (b) residual water saturation.

Figure 8 .
Figure 8. Variation of average water saturation with time for single-fracture models.

Figure 9 .
Figure 9.Effect of roughness on (a) bubble aspect ratio and (b) residual water saturation.

Figure 8 .
Figure 8. Variation of average water saturation with time for single-fracture models.

Figure 8 .
Figure 8. Variation of average water saturation with time for single-fracture models.

Figure 9 .Figure 9 .
Figure 9.Effect of roughness on (a) bubble aspect ratio and (b) residual water saturation.

Figure 11 .
Figure 11.Effect of tortuosity on (a) bubble aspect and (b) residual water saturation.

Figure 13 .
Figure 13.Effect of the gas/water ratio on (a) BAR and (b) residual water saturation.

Figure 14 .
Figure 14.Snapshots of the gas and water distribution in single-fracture models with flow rates of (a) 0.0002; (b) 0.002; (c) 0.004; and (d) 0.008 m/s.

Figure 13 .
Figure 13.Effect of the gas/water ratio on (a) BAR and (b) residual water saturation.

Figure 14 .
Figure 14.Snapshots of the gas and water distribution in single-fracture models with flow rates of (a) 0.0002; (b) 0.002; (c) 0.004; and (d) 0.008 m/s.

Figure 13 .
Figure 13.Effect of the gas/water ratio on (a) BAR and (b) residual water saturation.

Figure 13 .
Figure 13.Effect of the gas/water ratio on (a) BAR and (b) residual water saturation.

Figure 14 .
Figure 14.Snapshots of the gas and water distribution in single-fracture models with flow rates of (a) 0.0002; (b) 0.002; (c) 0.004; and (d) 0.008 m/s.

Figure 14 .
Figure 14.Snapshots of the gas and water distribution in single-fracture models with flow rates of (a) 0.0002; (b) 0.002; (c) 0.004; and (d) 0.008 m/s.

3 .
Gas/Water Flow Simulations in Rough-Walled Y-Shaped Junction Models 3.3.1.Effect of Aperture Ratio on the Phase Distribution

Figure 15 .
Figure 15.Snapshots of the gas and water distribution in Y-shaped junction models.(a-c) aperture ratio of 1.0; (d-f) aperture ratio of 1.5; (g-i) aperture ratio of 2.0.The left, middle, and right columns represent the snapshots at simulation steps of 0 s, 1.5 s, and 8.5 s, respectively.

Figure 15 .
Figure 15.Snapshots of the gas and water distribution in Y-shaped junction models.(a-c) aperture ratio of 1.0; (d-f) aperture ratio of 1.5; (g-i) aperture ratio of 2.0.The left, middle, and right columns represent the snapshots at simulation steps of 0 s, 1.5 s, and 8.5 s, respectively.

Figure 16 .
Figure 16.Flow velocity at outlets of the Y-shaped junction.(a) Aperture ratio = 1.5, upper outlet; (b) aperture ratio = 1.5, lower outlet; (c) aperture ratio = 2, upper outlet; (d) aperture ratio = 2, lower outlet.3.3.2.Effect of Intersecting Angle and Flow Rate on the Phase Distribution Figures 17 and 18 depict the phase distributions considering the intersecting angle and flow rate, respectively.As shown in the two figures, the intersecting angle and flow rate exert a negligible effect on the phase distribution.For all the cases considered in Figures 17 and 18, the gas phase flows out of the model only from the larger outlet channel, whereas the water phase flows through both the larger and smaller channels.

Figure 18 .
Figure 18.Snapshots of the phase distribution for the Y-shaped junction models with flow rates of (a) 0.0002 m/s; (b) 0.004 m/s.

Table 1 .
Simulation scenario setup for the single-fracture models.

Table 1 .
Simulation scenario setup for the single-fracture models.