Abstract

As a complex two-phase flow in naturally fractured coal formations, the prediction and analysis of CBM production remain challenging. This study presents a discrete fracture approach to modeling coalbed methane (CBM) and water flow in fractured coal reservoirs, particularly the influence of fracture orientation, fracture density, gravity, and fracture skeleton on fluid transport. The discrete fracture model is first verified by two water-flooding cases with multi- and single-fracture configurations. The verified model is then used to simulate CBM production from a discrete fractured reservoir using four different fracture patterns. The results indicate that fluid behavior is significantly affected by orientation, density, and fracture connectivity. Finally, several cases are performed to investigate the influence of gravity and fracture skeleton. The simulation results show that gas migrates upwards to the top reservoir during fluid extraction owing to buoyancy and the connected fracture skeleton plays a dominant role in fluid transport and methane production efficiency. Overall, the developed discrete fracture model provides a powerful tool to study two-phase flow in fractured coal reservoirs.

1. Introduction

Coalbed methane (CBM) as a form of high-quality clean energy has attracted considerable interest for sustainable development in both industrial and academic realms [1, 2]. Although several countries (e.g., China, Australia, USA, and Canada) have achieved commercial CBM production from coal reservoirs, the prediction and analysis of CBM production remain challenging because of complex two-phase flow in naturally fractured coal formations [37].

Coal is typically composed of matrix and fractures [8, 9]. The matrix refers to the collection of pores of different scales including micropores, mesopores, and macropores [7, 8, 10, 11]. The fracture system comprises four types of fractures: cleats, fracture swarms (or cracks/tertiary cleats), structure fractures (or faults), and hydraulic fractures [7]. Cleats refer to two sets of perpendicular fractures, called face and butt cleats. Fracture swarms and structure fractures refer to randomly distributed microfractures and large-scale fractures, respectively [7, 1214]. Hydraulic fractures are artificial fractures induced by industrial injection activities called hydrofracturing [1517]. The unique coal structure results in complex coupled fluid transport between matrix and fractures in a coal reservoir [18]. The coal matrix acts as a primary reservoir for CBM storage, although it has relatively low permeability and may even be impermeable. Large amounts of CBM are adsorbed on the inner surface of the coal matrix [19]. The fracture system provides a primary pathway for the migration of CBM and water through underground coal reservoirs. Storage and transfer of fluid are the essential properties of fractured coal reservoirs and are described by adsorption/absorption and diffusion and seepage models, respectively.

The complexity of pores and fracture networks complicates in situ analysis of water and methane transport, which introduces errors into the evaluation of methane production performance. Hence, the analysis of flow in fractured porous media is of great importance in fluid flow and thus methane production [20, 21]. A series of conceptual and numerical models across multiple scales has been developed and proposed to clarify transport mechanisms [2229]. In general, fluid in fractured porous media is approached with two conceptual models, which are continuum models and discrete fracture-matrix (DFM) models. In continuum models, fractures are represented implicitly in fractured porous media. The equivalent properties are calculated with crack tensor theory [30, 31], in which orientation, size, and aperture of fractures are considered. This type of upscaling technique is widely used for large-scale simulation, especially for the reservoirs with dense fracture networks. Generally, coal is treated as a structure with a single porosity and single permeability (SPSP), dual porosity and single permeability (DPSP), dual porosity and dual permeability (DPDP), or even triple porosity and dual permeability (TPDP) [29, 3236], in which matrix and fractures overlap. The representative elementary volume (REV) inside the reservoir is assumed to simultaneously satisfy the flow mass balance equations of matrix and fractures. The aforementioned methods use a continuum model or equivalent porous media for modeling fractured rocks or coal rock. Several well-known reservoir simulators, including ECLIPSE [37], CMG [38], COMET2 [39], and TOUGH2 [40], utilize the continuum models. Considering the dominant role of fractures in fluid transport, an alternative conceptual discrete fracture model has been proposed where the matrix is assumed to be impermeable and fluid processes are controlled by the fracture network [41]. In the discrete fracture model, the fractures are described explicitly by lower-dimensional lines or interface, which has the advantage of mesh generation and thus reduces the computational time greatly. In light of the mass exchange between matrix and fracture, a single-phase discrete fracture-matrix model has been developed to investigate the influence of adsorbed and free gas and fracture networks on gas production [42]. In this model, flow behavior of fluid occurs in both fracture and surrounding matrix system.

In this paper, we first develop a mathematical model to simulate water and methane flow through fractured coal reservoirs. Two water-flooding cases containing multiple fractures and a single fracture are then simulated to verify the accuracy of the proposed model. We then test four cases with multifracture configurations to investigate the influence of fracture orientation and distribution pattern on fluid behavior and methane production performance. Finally, we carried out several simulation cases with discrete fracture networks to investigate the effects of gravity and connectivity on fluid transport.

2. Governing Equations

2.1. Water and CBM Flow in Porous Media

In the mathematical model presented here, the coal reservoir is assumed to be saturated with methane and water in gas and aqueous phases, respectively. Hence, the sum of the saturated gas (nonwetting) phase and wetting (water) phase is equal to 1. Moreover, the model assumes that methane adsorbed on the coal grain surface diffuses instantaneously into the pores. Thus, the methane mass in the matrix system consists of free and adsorbed phases. The general mass balance equations for immiscible-phase (water and gas) flow in the coal reservoir matrix are given by gas phase pressure and water saturation , where the velocity of each phase is described by Darcy’s law. The governing equations for water and methane are described as follows. where is the absolute permeability of the matrix, and are the relative permeabilities of the water and gas phase, respectively, is the water pressure, is the density of each phase ( and refer to water and methane, resp.), which is calculated by [14], is the fluid compressibility, is the Langmuir volume constant, is the Langmuir pressure, is the density of the coal matrix, and is the porosity of the matrix system. The capillary pressure is the pressure difference between these two immiscible fluids, given as

2.2. Water and CBM Flow in Fractures

In this study, we represent fractures in the coal reservoir as low-dimensional grid cells [26]. Fractures are described as two-dimensional interfaces and one-dimensional lines in a three-dimensional or two-dimensional domain, respectively. A two-dimensional domain contains discontinuous fractures, as shown in Figure 1. The total simulation space is decomposed into where and represent the matrix and fracture domains, respectively, is the fracture aperture of fracture subdomain , and is the total number of fractures.

We assume pressure continuity across the whole model space, which means the gas pressure, water pressure, and capillary pressure are the same for the matrix and fracture grids, which means that the jump of pressure and saturation on the interface of matrix and fracture are not considered in this work. The equations for methane and water flow through fractures are expressed as where is the fracture aperture or thickness. The variables in equations (3)–(6) have the same physical characteristics and subscripted and represent these variables inside the matrix and fracture system, respectively. We demonstrate the detailed process of the weak form of equation (5). All items in equation (5) are first moved to the right-hand side, and both sides are multiplied by the test function for wetting saturation , integrating over the simulation domain :

According to Green’s first identity and divergence theorem, the third part of the right side in equation (7) is then expressed as

Finally, equation (8) is rearranged and the governing equation is obtained as

Similarly, we can obtain the weak expression for the water flow equation by multiplying equation (6) by the test function for water pressure , integrating over the simulation domain , and applying Green’s first identity and divergence theorem as

Equations (9) and (10) are referred to as the weak forms of water and CBM mass balance equations.

To solve equations (1), (2), (9), and (10), the auxiliary equations of capillary pressure, , gas and water relative permeability of the nonwetting and wetting phases are adopted as follows5: where and refer to variables inside the fracture and matrix systems, respectively, is the entry pressure, and and are coefficients determined by laboratory experiments. The effective saturation is defined as where and represent the residual saturations of the water and gas, respectively.

The initial condition for the gas pressure and water saturation is

As boundary conditions, the two-phase flow can have the following.

The Dirichlet boundary conditions for the gas pressure and water saturation are given as

The flux conditions, called as natural boundary conditions, which are included in the weak form of two-phase equations

3. Model Verification

We solve the above equations with the finite element software COMSOL. The equations of flow mass balance in the matrix are implemented with equations (1) and (2) using the partial differential equation (PDE) interface. The two immiscible phase flow in fractures with equations (9) and (10) are implemented with a weak contribution module. We then test two configurations of water flooding in an oil reservoir to investigate the accuracy of the model and numerical solution proposed in the paper. Two cases with different fracture configurations are adopted as follows.

(1) Multifracture Case. Figure 2 depicts the model geometry and mesh scheme. In this case, water is injected into a fractured porous medium with six fractures for 25 days. Detailed information of these fractures is provided in the references [26].

(2) Single-Fracture Case. In the simulation region, a single fracture with an arbitrary angle is modeled. Simulations were performed with three fracture orientations and to investigate the influence of fracture angle on flow behavior. The simulation time was 50 days.

The model regions in the two configurations are . The domain is initially nearly saturated with oil. The porosity and permeability of matrix are 0.20 and (1 millidarcy), respectively. The fracture porosity is 1. All fractures in the domain are assumed to have the same aperture of . Based on the cubic law, the corresponding permeability of the fractures is . The density and viscosity of the wetting phase are 1000 kg/m3 and , respectively. The density and viscosity of the nonwetting phase are 600 kg/m3 and . Fluid compressibility is neglected for both phases, which is justifiable because flow velocities are very small. The injection and production wells are located in the lower left and upper right corners. Water is injected into the fractured porous media at a constant rate of . The initial pressure and nonwetting phase saturation are set to 3.99 MPa and 0.99, respectively.

Pointwise constraints are applied at the production well with constant pressure and saturation. All boundaries are impermeable. The capillary pressure and relative permeability of water and oil are described as a function of water saturation , shown as [26, 43] where the parameter in the matrix and fracture is equal to 1 atm.

The spatial distribution of water saturation after 25 days of water injection into a nearly saturated oil reservoir is shown in Figure 3. A good match is achieved between reference models and our simulation results, which demonstrates the accuracy of the proposed model.

Figure 4 shows the spatial distribution of water saturation and gas pressure after 50 days of water injection with different fracture angles. The results obtained from the numerical simulations are in good agreement with reference studies22,41. Comparisons of water saturation and pressure along the diagonal line from the injection well to the production well after 50 days of water injection with low-dimensional discrete fracture model (L-DFM) and equidimensional discrete fracture model (E-DFM) are shown in Figure 5. It can be seen that the result of L-DFM is in good agreement with the result of E-DFM, which indicates that the L-DFM proposed in the paper can accurately simulate the two-phase flow in fractured porous media. For the large absolute permeability of a fracture, the fluid preferentially propagates into the reservoir through the fracture and causes significant pressure changes. A steady-state flow along the fracture is observed in the case with fracture angle (blue lines in Figure 5). The simulation results illustrate that the fluid flow behavior is largely controlled by the angles of the fractures.

The curves for oil cumulative production with different fracture orientations are shown in Figure 6. The cumulative production is the lowest in the case with fracture angle , which can be explained by the fact that injected water prefers to migrate aligned with fracture orientation, and thus, less oil is pushed out of the reservoir by water injection.

4. Simulation Cases

4.1. CBM Production from Discrete Fractured Reservoirs

In this section, we introduce the four cases with different patterns (Figure 7) that were tested to investigate the influence of fractures on flow fluid behavior and methane production. In the first two cases, 45 parallel fractures are uniformly distributed through the entire coal reservoir at orientations of , respectively. In the third and fourth cases, there are two sets of orthogonal fractures with angles and , respectively. The simulation domain is , in which the aperture and length of all the fractures are assumed to be 10−4 m and 5 m. The production well is located at the center of the simulation model with a constant gas pressure of 1 MPa and water saturation of 0.2. The initial gas pressure and water saturation are 6 MPa and 0.7, respectively. The total simulation time is 100 days. The surrounding boundaries are set to no flow. Other simulation parameters are listed in Table 1.

Figure 8 shows the spatial distributions of water saturation after 100 days of production for the four cases with different fracture configurations. The simulation results show that fracture geometry has a critical influence on flow path. Figure 9 shows the gas pressure and water saturation along the vertical and horizontal lines. During production, the water saturation and gas pressure decrease from the outer lateral boundaries to the production well. The speed of saturation and pressure front extraction from the coal seam differ for the four cases. In case 4, the pressure and saturation along the lines are lower than the initial conditions, which signify that drainage has approached the surrounding boundaries.

The temporal evolution of gas pressure and water saturation at points A and B is shown in Figure 10. A decrease in pore pressure and saturation is observed in the early stage of all cases because of the pressure and saturation drawdown at the production well. The water saturation and gas pressure in cases 3 and 4 are lower than those in cases 1 and 2 likely owing to the increased density (or number) of fractures, which enhances the overall reservoir permeability and fluid velocity. In case 4, the fractures in the vertical direction coincidentally connect to form a long fracture with a length of 45 m. Fluid migration is the fastest in case 4, which demonstrates that fracture connectivity dominantly impacts fluid transport and production efficiency.

4.2. Sensitive Analysis
4.2.1. Effect of Gravity

In this section, we set up two simulation cases with three-dimensional models that consider two sets of orthogonal fractures with angles . The model geometry is shown in Figure 11. The apertures of the whole fractures is assumed to be 10−4 m. The length and position of the two sets of fractures are randomly distributed in the simulation domain. Two simulation cases are performed to investigate the effect of gravity on fluid migration. Gravity is neglected in case 1 and considered in case 2. In these simulations, the coal height is 10 m and the simulation time is 50 days. Other settings and parameters are the same as those in Section 4.1.

Figure 12 shows the spatial distribution of water saturation of the whole domain (Figures 12(a) and 12(b)), inner surfaces of fractures (Figures 12(c) and 12(d)), and surface monitoring after 50 days of production in the three-dimensional reservoir with randomly distributed fractures (Figure 12(d)). The water saturation in the reservoir decreases extensively during drainage gas production. Without considering gravity, the water saturation is uniformly distributed in the vertical direction, whereas a nonuniform saturation distribution is observed along the fractures and monitoring surface in case 2. Gas saturation after 50 days of production along lines in the - and -directions is shown in Figure 13. The gas saturation exhibits a “wave-type” reduction from the wellbore to the lateral boundaries. Gas saturation along lines d, e, and f in the -direction is smoother than that along lines a, b, and c in the -direction because of fewer fractures cross the -direction lines. The gas saturation is the largest along the upper lines (a, d) and lowest along lower lines (c, f) as a result of the buoyancy effect.

Figure 14 shows the temporal evolution of gas saturation at point b in case 1 and points a, b, and c in case 2. The gas saturation increases rapidly in a relatively short time as a result of continuous dewatering. Gas continues to migrate upwards, which causes more gas to gather at the top reservoir surface during production (Figures 13 and 14).

4.2.2. Effect of Fracture Skeleton

This section carries our several cases to study the influence of fracture skeleton on the migration path of fluid. The discrete fracture network is generated by the open source tool DFNE [44]. Two sets of fractures are oriented of 45 and 135 degrees in a coal reservoir. Each set of fractures has 40 individual fractures. The minimum and maximum lengths of fracture lines are 1 and 7 m, respectively. Then, disconnected and isolated fractures are removed to investigate the effect of fracture skeleton. The original fracture network (Case a) and a connected fracture skeleton after processing (Case b) are shown in Figure 15. The total simulation is and other parameters are listed in Table 1.

The spatial distribution of methane saturation after with different fracture permeabilities of , , and in two different simulation domains is shown in Figure 16. Generally, the simulation results show that the distributions of gas saturation in two domains are similar. The distribution of saturation in Case b is smoother than that in Case a. The fact can be explained by that a large saturation gradient appears at the end of the disconnected fracture due to large permeability between fracture and matrix.

Figure 17 shows the evolution of average gas-phase pressure () along the vertical line for two different cases with different fracture permeabilities of , , and . The pressure decreases as the methane is continuously produced out of the fractured coal reservoir. Pressure drops significantly in the case with larger fracture permeability. In Figure 17, solid lines are simulation results in Case a while the dashed lines are the results in Case b. It can be seen that a good agreement has been achieved between those two cases, which demonstrates that the skeleton of fracture networks has an influential contribution to methane production. A larger difference is observed for the two cases with a higher fracture permeability of . Several reasons, involving boundary effect and fracture numbers according to the vertical line, lead to the phenomenon.

5. Conclusions

In this study, we developed and applied a discrete fracture model to simulate two-phase (coalbed methane and water) flow through fractured coal reservoirs. The proposed two-phase model in fractured porous media was verified by two oil-reservoir water-flooding cases with single and multiple fractures. The simulation results are in good agreement, which confirms the model feasibility and accuracy. We simulated CBM production from discrete fractured coal reservoirs with four types of fracture configurations. The simulation results clearly show that the patterns of fluid flow and production performance are significantly affected by fracture orientation, density, and connectivity. The fluid prefers to migrate aligned with the fracture orientation. Increasing fracture density enhances production efficiency. Moreover, fracture connectivity seems to contribute significantly to fluid transport and methane production efficiency. Later, two three-dimensional cases were studied to investigate the influence of gravity. The results show that gas continues to migrate upwards to the top reservoir surface during fluid extraction as a result of the buoyancy of methane, which provides the possibility of methane leakage. Finally, we performed two cases of original discrete fracture network and a connected fracture network to study the effect of fracture skeleton. Simulation results demonstrate that the connected fracture skeleton is of great importance to fluid migration and methane production. Overall, the developed model provides a powerful approach to study coalbed methane and water flow in fractured coal reservoirs.

Data Availability

The data used to support the findings of this study are available from the first author upon request.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

This study was supported by the Special Subject Grant of National Basic Research Program of China (973 program) (No. 2015CB251602), the National Science and Technology Major (2016ZX05043), the Jiangsu Natural Science Foundation (BK20180636), the Independent Innovation Project for Double First-Level Construction of CUMT (2018ZZCX04), the Advance Research Program (LTKY201803), and the China and Jiangsu Postdoctoral Science Foundation (2019M65201). We also thank Curtis Oldenburg (Lawrence Berkeley National Laboratory) for comments on an earlier draft.