A Discrete Fracture Modeling Approach for Analysis of Coalbed Methane and Water Flow in a Fractured Coal Reservoir

Key Laboratory of CBM Resources and Dynamic Accumulation Process, Ministry of Education, China University of Mining and Technology, Xuzhou, Jiangsu 221116, China School of Mechanics and Civil Engineering, China University of Mining and Technology, Xuzhou, Jiangsu 221116, China Lawrence Berkeley National Laboratory, Energy Geosciences Division, Berkeley, CA 94720, USA Chinese Academy of Geological Sciences, Beijing 100037, China State Key Laboratory Geomechanics and Deep Underground Engineering, China University of Mining and Technology, Xuzhou, Jiangsu 221116, China


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 [3][4][5][6][7].
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,[12][13][14]. Hydraulic fractures are artificial fractures induced by industrial injection activities called hydrofracturing [15][16][17]. 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 [22][23][24][25][26][27][28][29]. 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,[32][33][34][35][36], 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 wellknown 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 singlephase 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.

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 S mg and wetting (water) phase S mw 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 p mg and water saturation S mw , where the velocity of each phase v α is described by Darcy's law. The governing equations for water and methane are described as follows.
where K m is the absolute permeability of the matrix, k mrw and k mrg are the relative permeabilities of the water and gas phase, respectively, p mw is the water pressure, ρ mα is the density of each phase (α = w and g refer to water and methane, resp.), which is calculated by 1/C m ðdρ mα /dp mα Þ [14], C mα is the fluid compressibility, V mL is the Langmuir volume constant, p mL is the Langmuir pressure, ρ mc is the density of the coal matrix, and ϕ m is the porosity of the matrix system. The capillary pressure p mc 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 Ω m and Ω f represent the matrix and fracture domains, respectively, d f i is the fracture aperture of i th fracture subdomain Ω f i , and n 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 d f is the fracture aperture or thickness. The variables in equations (3)-(6) have the same physical characteristics and subscripted m and f 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 f S f w , integrating over the simulation domain Ω f : 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 f p f g , integrating over the simulation domain Ω f , and applying Green's first identity and divergence theorem as

Geofluids
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, p βc , gas and water relative permeability of the nonwetting k βrg and wetting k βrw phases are adopted as follows 5 : where β = f and m refer to variables inside the fracture and matrix systems, respectively, p βc is the entry pressure, and λ β and m β are coefficients determined by laboratory experiments. The effective saturation S βe is defined as where S βwr and S βgr represent the residual saturations of the water and gas, respectively.

Geofluids
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

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 numer-ical 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].    7 Geofluids permeability of water k rw and oil k rn are described as a function of water saturation S w , shown as [26,43] where the parameter B 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 differ-ent fracture angles. The results obtained from the numerical simulations are in good agreement with reference studies 22,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 twophase 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 θ = 45°(blue lines in Figure 5). The simulation results illustrate that the fluid flow behavior is largely controlled by the angles of the fractures.

Geofluids
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 θ = 45°, 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.

CBM Production from Discrete Fractured Reservoirs.
In this section, we introduce the four cases with different pat-terns ( 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 θ = 0°and 45°, respectively. In the third and fourth cases, there are two sets of orthogonal fractures with angles θ = 45°and − 45°and θ = 0°and 90°, respectively. The simulation domain is 50 × 50 m, 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.

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

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 θ = 0°and 90°. The model geometry is shown in Figure 11. The apertures of the 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 threedimensional 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 xand y-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 y-direction is smoother than that along lines a, b, and c in the x-direction because of fewer fractures cross the y-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).  Figure 15. The total simulation is 2:0 × 10 4 s and other parameters are listed in Table 1.

Geofluids
The spatial distribution of methane saturation after 2:0 × 10 4 s with different fracture permeabilities of k f = 5:0 × 10 −17 , 1:0 × 10 −16 , and 5:0 × 10 −16 m 2 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 (p g = Ð l p g dl/l) along the vertical line x = 4:5 for two 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 5:0 × 10 −16 m 2 . Several reasons, involving boundary effect and fracture numbers according to the vertical line, lead to the phenomenon.

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 twophase 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 investi-gate 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.