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

Authors: Tianran Ma; Hao Xu; Chaobin Guo; Xuehai Fu; Weiqun Liu; Rui Yang
Journal: Geofluids
 (2020)

---

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


---

## Body

## 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 [3–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–14]. Hydraulic fractures are artificial fractures induced by industrial injection activities called hydrofracturing [15–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–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–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 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) phaseSmg and wetting (water) phase Smw 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 pmg and water saturation Smw, where the velocity of each phase vα is described by Darcy’s law. The governing equations for water and methane are described as follows.
(1)ϕmρmwSmwCmw∂pmg∂t+ϕmρmw−ϕmρmwSmwCmw∂pmc∂Smw∂Smw∂t+∇·−ρmwKmkmrwμmw∇pmg−∂pmc∂Smw∇Smw+ρmwg=0,(2)ϕmρmgSmgCmg+ρmcVmLpmLpmg+pmL2∂pmg∂t−ϕmρmg∂Smw∂t+∇·−ρmgKmkmrgμmg∇pmg+ρmgg=0,where Km is the absolute permeability of the matrix,kmrw and kmrgare the relative permeabilities of the water and gas phase, respectively, pmw 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/Cmdρmα/dpmα [14], Cmα is the fluid compressibility, VmL is the Langmuir volume constant, pmL is the Langmuir pressure, ρmc is the density of the coal matrix, and ϕm is the porosity of the matrix system. The capillary pressure pmc is the pressure difference between these two immiscible fluids, given as
(3)pmc=pmg−pmw.


### 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
(4)Ω=Ωm+Ωf=Ωm+∑i=0ndfiΩfi,where Ωm and Ωf represent the matrix and fracture domains, respectively, dfi is the fracture aperture of ith fracture subdomain Ωfi, and n is the total number of fractures.Figure 1
Schematic of the matrixΩm and fractures Ωf1~Ωf4 subdomains.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(5)dfϕfρfwSfwCfw∂pfg∂t+dfϕfρfw−ϕfρfwSfwCfw∂pfc∂Sfw∂Sfw∂t+∇τ·−dfρfwKkfrwμfw∇τpfg−∂pfc∂Sfw∇τSfw+ρfwg=0,(6)dfϕfρfgSfgCfg+ρfcVfLpfLpfg+pfL2∂pfg∂t−dfϕρfg∂Sfw∂t+∇τ·−dfρfgKfkfrgμfg∇τpfg+ρfgg=0,where df 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 Sfw~, integrating over the simulation domain Ωf:
(7)0=−∫ΩfdfϕfρfwSfwCfw∂pfg∂t+dfϕfρfw−ϕfρfwSfwCfw∂pfc∂Sfw∂Sfw∂tSfw~dΩf−∫Ω∇τ·−dfρfwKfkfrwμfw∇τpfg−∂pfc∂Sfw∇τSfw+ρfwgSfw~dΩf.According to Green’s first identity and divergence theorem, the third part of the right side in equation (7) is then expressed as
(8)∫Ωf∇·−dfρfwKfkfrwμfw∇τpfg−∂pfc∂Sfw∇τSfw+ρfwgSfw~dΩf=∫n−dfρfwKfkfrwμfw∇τpfg−∂pfc∂Sfw∇τSfw+ρfwg·nSfw~dn−∫Ωf∇τSfw~·−dfρfwKfkfrwμfw∇τpfg−∂pfc∂Sfw∇τSfw+ρfwgdΩf.Finally, equation (8) is rearranged and the governing equation is obtained as
(9)0=−dfϕfρfwSfwCfw∂pfg∂t+dfϕfρfw−ϕfρfwSfwCfw∂pfc∂Sfw∂Sfw∂tSfw~dΩf−∫n−dfρfwKfkfrwμfw∇τpfg−∂pfc∂Sfw∇τSfw+ρfwg·nSfw~dn+∫Ωf∇τSfw~·−dfρfwKfkfrwμfw∇τpfg−∂pfc∂Sfw∇τSfw+ρfwgdΩf.Similarly, we can obtain the weak expression for the water flow equation by multiplying equation (6) by the test function for water pressure pfg~, integrating over the simulation domain Ωf, and applying Green’s first identity and divergence theorem as
(10)0=∫Ωf−dfϕfρfgSfgCfg+ρfcVfLpfLpfg+pfL2∂pfg∂tpfg~+dfϕfρfg∂Sfw∂tpfg~dΩf−∫n−dfρfgKfkfrgμfg∇τpfg+ρfgg·npfg~dn+∫Ωf∇τpfg~·−dfρfgKfkfrgμfg∇τpfg+ρfggdΩf.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 follows5:
(11)pβc=pβcSβe−1/λβkβrg=1−Sβe21−Sβe2,kβrw=Sβe1−1−Sβe−1/mβmβ2,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
(12)Sβe=Sβw−Sβwr1−Sβwr−Sβgr,where Sβwr and Sβgr represent the residual saturations of the water and gas, respectively.The initial condition for the gas pressure and water saturation is(13)pβgt=0=pβg0,inΩ,Sβwt=0=Sβw0,inΩ.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(14)pβgΓ=pβgD,SβwΓ=SβWD.The flux conditions, called as natural boundary conditions, which are included in the weak form of two-phase equations(15)ρβgVg·nΓ=qg,ρβwVw·nΓ=qw.


## 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) phaseSmg and wetting (water) phase Smw 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 pmg and water saturation Smw, where the velocity of each phase vα is described by Darcy’s law. The governing equations for water and methane are described as follows.
(1)ϕmρmwSmwCmw∂pmg∂t+ϕmρmw−ϕmρmwSmwCmw∂pmc∂Smw∂Smw∂t+∇·−ρmwKmkmrwμmw∇pmg−∂pmc∂Smw∇Smw+ρmwg=0,(2)ϕmρmgSmgCmg+ρmcVmLpmLpmg+pmL2∂pmg∂t−ϕmρmg∂Smw∂t+∇·−ρmgKmkmrgμmg∇pmg+ρmgg=0,where Km is the absolute permeability of the matrix,kmrw and kmrgare the relative permeabilities of the water and gas phase, respectively, pmw 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/Cmdρmα/dpmα [14], Cmα is the fluid compressibility, VmL is the Langmuir volume constant, pmL is the Langmuir pressure, ρmc is the density of the coal matrix, and ϕm is the porosity of the matrix system. The capillary pressure pmc is the pressure difference between these two immiscible fluids, given as
(3)pmc=pmg−pmw.


## 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
(4)Ω=Ωm+Ωf=Ωm+∑i=0ndfiΩfi,where Ωm and Ωf represent the matrix and fracture domains, respectively, dfi is the fracture aperture of ith fracture subdomain Ωfi, and n is the total number of fractures.Figure 1
Schematic of the matrixΩm and fractures Ωf1~Ωf4 subdomains.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(5)dfϕfρfwSfwCfw∂pfg∂t+dfϕfρfw−ϕfρfwSfwCfw∂pfc∂Sfw∂Sfw∂t+∇τ·−dfρfwKkfrwμfw∇τpfg−∂pfc∂Sfw∇τSfw+ρfwg=0,(6)dfϕfρfgSfgCfg+ρfcVfLpfLpfg+pfL2∂pfg∂t−dfϕρfg∂Sfw∂t+∇τ·−dfρfgKfkfrgμfg∇τpfg+ρfgg=0,where df 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 Sfw~, integrating over the simulation domain Ωf:
(7)0=−∫ΩfdfϕfρfwSfwCfw∂pfg∂t+dfϕfρfw−ϕfρfwSfwCfw∂pfc∂Sfw∂Sfw∂tSfw~dΩf−∫Ω∇τ·−dfρfwKfkfrwμfw∇τpfg−∂pfc∂Sfw∇τSfw+ρfwgSfw~dΩf.According to Green’s first identity and divergence theorem, the third part of the right side in equation (7) is then expressed as
(8)∫Ωf∇·−dfρfwKfkfrwμfw∇τpfg−∂pfc∂Sfw∇τSfw+ρfwgSfw~dΩf=∫n−dfρfwKfkfrwμfw∇τpfg−∂pfc∂Sfw∇τSfw+ρfwg·nSfw~dn−∫Ωf∇τSfw~·−dfρfwKfkfrwμfw∇τpfg−∂pfc∂Sfw∇τSfw+ρfwgdΩf.Finally, equation (8) is rearranged and the governing equation is obtained as
(9)0=−dfϕfρfwSfwCfw∂pfg∂t+dfϕfρfw−ϕfρfwSfwCfw∂pfc∂Sfw∂Sfw∂tSfw~dΩf−∫n−dfρfwKfkfrwμfw∇τpfg−∂pfc∂Sfw∇τSfw+ρfwg·nSfw~dn+∫Ωf∇τSfw~·−dfρfwKfkfrwμfw∇τpfg−∂pfc∂Sfw∇τSfw+ρfwgdΩf.Similarly, we can obtain the weak expression for the water flow equation by multiplying equation (6) by the test function for water pressure pfg~, integrating over the simulation domain Ωf, and applying Green’s first identity and divergence theorem as
(10)0=∫Ωf−dfϕfρfgSfgCfg+ρfcVfLpfLpfg+pfL2∂pfg∂tpfg~+dfϕfρfg∂Sfw∂tpfg~dΩf−∫n−dfρfgKfkfrgμfg∇τpfg+ρfgg·npfg~dn+∫Ωf∇τpfg~·−dfρfgKfkfrgμfg∇τpfg+ρfggdΩf.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 follows5:
(11)pβc=pβcSβe−1/λβkβrg=1−Sβe21−Sβe2,kβrw=Sβe1−1−Sβe−1/mβmβ2,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
(12)Sβe=Sβw−Sβwr1−Sβwr−Sβgr,where Sβwr and Sβgr represent the residual saturations of the water and gas, respectively.The initial condition for the gas pressure and water saturation is(13)pβgt=0=pβg0,inΩ,Sβwt=0=Sβw0,inΩ.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(14)pβgΓ=pβgD,SβwΓ=SβWD.The flux conditions, called as natural boundary conditions, which are included in the weak form of two-phase equations(15)ρβgVg·nΓ=qg,ρβwVw·nΓ=qw.


## 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].Figure 2
Scheme depicting geometry, mesh boundary, and initial conditions of the simulation model.(2) Single-Fracture Case. In the simulation region, a single fracture with an arbitrary angle is modeled. Simulations were performed with three fracture orientations θ=0°,45°, and 135° to investigate the influence of fracture angle on flow behavior. The simulation time was 50 days.The model regions in the two configurations are1m×1m. The domain is initially nearly saturated with oil. The porosity and permeability of matrix are 0.20 and 9.87×10−16m2 (1 millidarcy), respectively. The fracture porosity is 1. All fractures in the domain are assumed to have the same aperture of 1.00×10−4m. Based on the cubic law, the corresponding permeability of the fractures is 8.33×10−10m2. The density and viscosity of the wetting phase are 1000 kg/m3 and 1×10−3Pa·s, respectively. The density and viscosity of the nonwetting phase are 600 kg/m3 and 0.45×10−3Pa·s. 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 2.32×10−8m3/s. 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 pressurepc and relative permeability of water krw and oil krn are described as a function of water saturation Sw, shown as [26, 43]
(16)pc=−BlogSw,krw=Sw2,krn=1−Sw2,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 Figure3. A good match is achieved between reference models and our simulation results, which demonstrates the accuracy of the proposed model.Figure 3
Distribution of water saturation after 25 days of water injection into a reservoir initially saturated with oil. Previous modeling results22 and our simulation results are displayed on (a) and (b), respectively.Figure4 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 θ=45° (blue lines in Figure 5). The simulation results illustrate that the fluid flow behavior is largely controlled by the angles of the fractures.Figure 4
Spatial distribution of (a–c) saturation and (d–f) pressure after 50 days of water injection with different fracture orientations.
(a)(b)(c)(d)(e)(f)Figure 5
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 L-DFM and E-DFM.
(a)(b)The curves for oil cumulative production with different fracture orientations are shown in Figure6. 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.Figure 6
Cumulative production with different fracture orientations.


## 4. Simulation Cases

### 4.1. CBM Production from Discrete Fractured Reservoirs

In this section, we introduce the four cases with different patterns (Figure7) 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°and45°, respectively. In the third and fourth cases, there are two sets of orthogonal fractures with angles θ=45°and−45° and θ=0°and90°, respectively. The simulation domain is 50×50m, 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 7
Model domain, fracture configuration, monitoring lines, and points for four cases. In cases 1 and 2, 45 fractures uniformly distributed through the entire coal reservoir are aligned at orientations ofθ=0°and45°, respectively, and cases 3 and 4 contain two sets of orthogonal fractures with angles θ=45°and−45° and θ=0°and90°, respectively.Table 1
Simulation parameters for water-flooding cases.
ParametersValueUnitSimulation area50×50m2Permeability of fracture,Kf1.00×10−10m2Porosity of fracture,ϕf1.00Fracture aperture,df1.00×10−4mPermeability of matrix,Km1.00×10−16m2Porosity of matrix,ϕm0.20Initial reservoir pressure,pfg and pmg1000PaInitial water saturation,Sfg and Smg0.70Entry capillary pressure,pfe and pme0.10MPaInitial residual water saturation0.00Initial residual nonwetting saturation0.00Viscosity of water,μfw and μmw10-3Pa·sViscosity of nonwetting phase,μfg and μmg1.84×10−5Pa·sDensity of water,ρfw and ρmw1000kg/m3Density of nonwetting phase,ρfg and ρmg0.864kg/m3Figure8 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.Figure 8
Spatial distribution of water saturation after 100 days of production for different fracture patterns.
(a)(b)(c)(d)Figure 9
(a, b) Gas pressure and (c, d) water saturation along the horizontal and vertical lines, respectively, after 250 days of production in the four cases.
(a)(b)(c)(d)The temporal evolution of gas pressure and water saturation at points A and B is shown in Figure10. 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.Figure 10
Temporal evolution of (a) gas pressure and (b) water saturation at points A and B for the four cases. Solid and dashed lines represent the results at points A and B, respectively.
(a)(b)


### 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θ=0°and90°. 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 11
Scheme depicting model domain, fracture configuration, monitoring surface, lines, and points for the three-dimensional cases.Figure12 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 x- and 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 12
Water saturation after 50 days of production from the three-dimensional model with randomly distributed fractures. The results for case 2 with consideration of gravity are shown in (b), (d), and (f).Figure 13
Gas saturation along the lines in the (a)x- and (b) y-directions after 50 days of production.
(a)(b)Figure14 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 14
Temporal evolution of gas saturation at point b in case 1 and points a, b, and c in case 2.


#### 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 10m×10m 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 2.0×104s and other parameters are listed in Table 1.Figure 15
The original fracture network (Case a) and connected fracture skeleton (Case b).
(a)(b)The spatial distribution of methane saturation after2.0×104s with different fracture permeabilities of kf=5.0×10−17, 1.0×10−16, and 5.0×10−16m2 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 16
The spatial distribution of methane saturation after2.0×10s in Case a and Case b with different fracture permeabilities of kf=5.0×10−17, 1.0×10−16, and 5.0×10−16m2.
(a)(b)Figure17 shows the evolution of average gas-phase pressure (pg¯=∫lpgdl/l) along the vertical line x=4.5 for two different cases with different fracture permeabilities of 5.0×10−17, 1.0×10−16, and 5.0×10−16m2. 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−16m2. Several reasons, involving boundary effect and fracture numbers according to the vertical line, lead to the phenomenon.Figure 17
The temporal evolution of average gas pressure along the vertical linex=0.45 with different fracture permeabilities of kf=5.0×10−17, 1.0×10−16, and 5.0×10−16m2.


## 4.1. CBM Production from Discrete Fractured Reservoirs

In this section, we introduce the four cases with different patterns (Figure7) 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°and45°, respectively. In the third and fourth cases, there are two sets of orthogonal fractures with angles θ=45°and−45° and θ=0°and90°, respectively. The simulation domain is 50×50m, 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 7
Model domain, fracture configuration, monitoring lines, and points for four cases. In cases 1 and 2, 45 fractures uniformly distributed through the entire coal reservoir are aligned at orientations ofθ=0°and45°, respectively, and cases 3 and 4 contain two sets of orthogonal fractures with angles θ=45°and−45° and θ=0°and90°, respectively.Table 1
Simulation parameters for water-flooding cases.
ParametersValueUnitSimulation area50×50m2Permeability of fracture,Kf1.00×10−10m2Porosity of fracture,ϕf1.00Fracture aperture,df1.00×10−4mPermeability of matrix,Km1.00×10−16m2Porosity of matrix,ϕm0.20Initial reservoir pressure,pfg and pmg1000PaInitial water saturation,Sfg and Smg0.70Entry capillary pressure,pfe and pme0.10MPaInitial residual water saturation0.00Initial residual nonwetting saturation0.00Viscosity of water,μfw and μmw10-3Pa·sViscosity of nonwetting phase,μfg and μmg1.84×10−5Pa·sDensity of water,ρfw and ρmw1000kg/m3Density of nonwetting phase,ρfg and ρmg0.864kg/m3Figure8 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.Figure 8
Spatial distribution of water saturation after 100 days of production for different fracture patterns.
(a)(b)(c)(d)Figure 9
(a, b) Gas pressure and (c, d) water saturation along the horizontal and vertical lines, respectively, after 250 days of production in the four cases.
(a)(b)(c)(d)The temporal evolution of gas pressure and water saturation at points A and B is shown in Figure10. 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.Figure 10
Temporal evolution of (a) gas pressure and (b) water saturation at points A and B for the four cases. Solid and dashed lines represent the results at points A and B, respectively.
(a)(b)


## 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θ=0°and90°. 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 11
Scheme depicting model domain, fracture configuration, monitoring surface, lines, and points for the three-dimensional cases.Figure12 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 x- and 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 12
Water saturation after 50 days of production from the three-dimensional model with randomly distributed fractures. The results for case 2 with consideration of gravity are shown in (b), (d), and (f).Figure 13
Gas saturation along the lines in the (a)x- and (b) y-directions after 50 days of production.
(a)(b)Figure14 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 14
Temporal evolution of gas saturation at point b in case 1 and points a, b, and c in case 2.


### 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 10m×10m 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 2.0×104s and other parameters are listed in Table 1.Figure 15
The original fracture network (Case a) and connected fracture skeleton (Case b).
(a)(b)The spatial distribution of methane saturation after2.0×104s with different fracture permeabilities of kf=5.0×10−17, 1.0×10−16, and 5.0×10−16m2 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 16
The spatial distribution of methane saturation after2.0×10s in Case a and Case b with different fracture permeabilities of kf=5.0×10−17, 1.0×10−16, and 5.0×10−16m2.
(a)(b)Figure17 shows the evolution of average gas-phase pressure (pg¯=∫lpgdl/l) along the vertical line x=4.5 for two different cases with different fracture permeabilities of 5.0×10−17, 1.0×10−16, and 5.0×10−16m2. 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−16m2. Several reasons, involving boundary effect and fracture numbers according to the vertical line, lead to the phenomenon.Figure 17
The temporal evolution of average gas pressure along the vertical linex=0.45 with different fracture permeabilities of kf=5.0×10−17, 1.0×10−16, and 5.0×10−16m2.


## 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θ=0°and90°. 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 11
Scheme depicting model domain, fracture configuration, monitoring surface, lines, and points for the three-dimensional cases.Figure12 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 x- and 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 12
Water saturation after 50 days of production from the three-dimensional model with randomly distributed fractures. The results for case 2 with consideration of gravity are shown in (b), (d), and (f).Figure 13
Gas saturation along the lines in the (a)x- and (b) y-directions after 50 days of production.
(a)(b)Figure14 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 14
Temporal evolution of gas saturation at point b in case 1 and points a, b, and c in case 2.


## 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 10m×10m 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 2.0×104s and other parameters are listed in Table 1.Figure 15
The original fracture network (Case a) and connected fracture skeleton (Case b).
(a)(b)The spatial distribution of methane saturation after2.0×104s with different fracture permeabilities of kf=5.0×10−17, 1.0×10−16, and 5.0×10−16m2 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 16
The spatial distribution of methane saturation after2.0×10s in Case a and Case b with different fracture permeabilities of kf=5.0×10−17, 1.0×10−16, and 5.0×10−16m2.
(a)(b)Figure17 shows the evolution of average gas-phase pressure (pg¯=∫lpgdl/l) along the vertical line x=4.5 for two different cases with different fracture permeabilities of 5.0×10−17, 1.0×10−16, and 5.0×10−16m2. 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−16m2. Several reasons, involving boundary effect and fracture numbers according to the vertical line, lead to the phenomenon.Figure 17
The temporal evolution of average gas pressure along the vertical linex=0.45 with different fracture permeabilities of kf=5.0×10−17, 1.0×10−16, and 5.0×10−16m2.


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



---

*Source: 8845348-2020-06-25.xml*