3D reconstruction of coal pore network and its application in CO₂-ECBM process simulation at laboratory scale

Abstract

Three-dimensional (3D) reconstruction of the equivalent pore network model (PNM) using X-ray computed tomography (CT) data are of significance for studying the CO₂-enhanced coalbed methane recovery (CO₂-ECBM). The docking among X-ray CT technology, MATLAB, with COMSOL software not only can realize the 3D reconstruction of PNM, but also the CO₂-ECBM process simulation. The results show that the Median filtering algorithm enabled the de-noising of the original 2D CT slices, the image segmentation of all slices was realized based on the selected threshold, and the PNM can be constructed based on the Maximum Sphere algorithm. The mathematical model of CO₂-ECBM process fully coupled the expanded Langmuir equation. At the same time for CO₂ injection, CH₄ pressure tends to decrease with the increase of CO₂ pressure, but its difference is not obvious. The CH₄ pressure in the slice center changed a lot, while at the edge it changed a little under different CO₂ pressures. The injected CO₂ was transported to matrix along the macro and micro-fractures with continuous flow. The injected CO₂ first replaced the adsorbed CH₄ by covering the inner surface of macro-pores and meso-pores to form the single molecular layer adsorption of CO₂. Then they migrated to micro-pores by Fick’s diffusion, sliding flow, and surface diffusion. Furthermore, the CO₂ replaced CH₄ adsorbed by volumetric filling in micro-pores, and formed the multi-molecular layer adsorption of CO₂. The gas pressure and migration path between CO₂ and CH₄ are opposite. This study can provide a theoretical basis for studying digital rock physics technology and enrich the development of CO₂-ECBM technology.

Graphical abstract

Keywords

CO₂-ECBM / 3D reconstruction / numerical simulation / X-ray CT / COMSOL / Qinshui Basin

1 Introduction

The CO₂-enhanced coalbed methane recovery (CO₂-ECBM) cannot only increase CH₄ production but also reduce CO₂ emissions (Fang et al., 2019a; Zheng et al., 2020). As a geological carrier for CO₂ injection, the pore system is the storage carrier of CO₂/CH₄, and the fracture system is the diffusion and seepage carrier of CO₂/CH₄ (Wang et al., 2018a, 2018b; Cheng et al., 2020; Wang et al., 2020a; Li et al., 2020a). Therefore, the three-dimensional (3D) reconstruction of coal pore and fracture can provide not only a framework for studying the seepage characteristic of fluid in porous media, but also theoretical support for the practical application of CO₂-ECBM technology.

As for pore and fracture, the characterization method can be summarized as follows. 1) Pore size distribution: low temperature CO₂/N₂ adsorption and mercury intrusion porosimetry (Connell et al., 2016; Li et al., 2017; Ni et al., 2020a; Nie et al., 2020). 2) Morphological observation: naked eye observation, optical microscope analysis, and scanning electron microscope analysis (Cai et al., 2018; Zhang et al., 2018; Li et al., 2020b). 3) 3D imaging: X-ray computed tomography (CT) and focused ion beam scanning electron microscopy (Fang et al., 2019b; Wang et al., 2020b). Through the comparison of the above experimental methods, it can be concluded that X-ray CT, as a non-destructive and high-resolution scanning techno-logy (Fang et al., 2020; Zhu et al., 2020), is the most direct and accurate method to reconstruct the pore and fracture structure.

Previous studies have shown that the X-ray CT technology has been used for the 3D reconstruction of coal structure (Fan et al., 2020), and the finite element software has been used for the CO₂-ECBM process simulation (Fang et al., 2019c, 2019d), but few researchers have performed the CO₂-ECBM process simulation in finite element software with pore and fracture structure as the geological support. Therefore, based on the data extracted by X-ray CT technology and the application of MATLAB software, the interconnected pores and throats can be extracted from the equivalent pore and fracture network model (PNM), and the STL file output can be imported into the finite element software for the CO₂-ECBM process simulation, which is a new research hotspot focus in the CO₂-ECBM direction (Fig. 1).

In this study, first, based on the data obtained by X-ray CT and the application of MATLAB software, the geological carrier needed for numerical simulation can be extracted from PNM. Secondly, the STL file of the interconnected PNM extracted from the MATLAB software can be imported into the COMSOL software (available at COMSOL website) for grid partitioning and debugging. Then, the CO₂-ECBM process can be numerically simulated in the COMSOL software. Finally, the dynamic characteristics of the CO₂-ECBM process and the constraining effect of the pore and fracture structure on this process are discussed (Fig. 1). the significance and innovation of this study can be reflected as follows: 1) the mathematical model for simulating the CO₂-ECBM process at laboratory scale was deduced. 2) The pore network model was reconstructed and visualized in 3D. 3) The comprehensive application of X-ray CT technology, MATLAB, and COMSOL software for the CO₂-ECBM process simulation was established. 4) The dynamic characteristics of the CO₂-ECBM process and the confining effect of the pore and fracture structure on such a process were analyzed and visualized. This study can provide a theoretical basis for the study of digital rock physics technology and enrich the development of CO₂-ECBM technology.

2 Geological background and sampling point distribution

2.1 Geological background of the study area

The research area is located in the southern part of the Qinshui Basin, which is bounded by the Jinhuo Fault Zone to the east and adjacent to the Taihang Mountain Uplift, the Huoshan Uplift, and the Wutai Mountain Uplift to the west, south, and north direction, respectively (Figs. 2(a)–2(b); Cai et al., 2011; Wang et al., 2019). The stratigraphic distribution has the typical characteristics of a synclinal basin, the outcrop layer at the margin is old and the outcrop layer in the basin is new (Wang et al., 2021; Zhang et al., 2021). The coal-bearing strata are mainly Carboniferous-Permian, and the 3# coal seam of the Shanxi Formation and the 15# coal seam of the Taiyuan Formation are the main mining coal seams in the study area, and their transverse distribution is relatively stable (Fang et al., 2017; Huang et al., 2019; Ni et al., 2020b), which provides favorable conditions for the implementation of the CO₂-ECBM project.

2.2 Sampling point

In this study, coal samples collected from Bofang (BF) colliery in the Qinshui Basin were taken as the research object (Fig. 2(b); Liu et al., 2017; Fang et al., 2017), and the packing, transportation, storage and basic testing of samples comply with relevant international standards (Liu et al., 2017; Fang et al., 2019e). The anti-oxidation treatment of the fresh samples collected can prevent the fresh samples from being oxidized (Fig. 2(c)), which will affect the accuracy of the physicochemical properties of sample. After testing, the key parameters of the coal samples are shown in Table 1.

3 Methods

3.1 Reconstruction method of pore and fracture structure

3.1.1 X-ray CT scan imaging

In this study, the X-ray CT scanning imaging system is composed of the X-ray source, precision sample table, high-resolution detector, data processing system and controller system (Fig. 3). The sample used for scanning is a small coal column with a diameter of 2 mm and a height of about 6 mm, which was drilled with a prototype mechanical drill (Fig. 3(a)). The scanning area is a small coal column with a diameter of 2 mm and a height of 1 mm, the total number of scans is 3600, the pixel resolution is 0.5 μm, the spatial resolution is 200 nm, and the typical 2D CT slices are shown in Fig. 3(c). The gray, black and white areas within the perimeter can represent the distribution of organic matter, pores and minerals in coal, respectively (Fig. 3(c)).

3.1.2 3D reconstruction method of pore and fracture structure

To extract the geological carrier needed for the CO₂-ECBM process simulation, i.e., to extract the interconnected PNM in coal, the scanned 2D CT slices must be visually reconstructed, and the main tasks are as follows: preprocessing of the 2D slices, threshold selection and image segmentation, analysis of the representative elementary volume (REV), and construction of the pore network model.

Preprocessing of 2D slices: Due to uncontrollable factors, such as the X-ray CT scanner, the original CT slices extracted under noise lead to certain errors in the reconstruction of the coal structure at the later stage. Previous studies have shown that the noise reduction, such as the Median Filtering algorithm, for the original 2D slices can well protect the integrity of pores and smooth the transition between pore and matrix (Li and Zhang 2019; Fang et al., 2020), which provides a good basic for 3D reconstruction of coal structure.

Threshold selection and image segmentation: Converting 2D slices into a 3D image is the purpose of threshold segmentation, so that the pores can be separated from the matrix. The image segmentation method based on the selected threshold value is widely used in image segmentation. Its core idea is to determine the threshold value for image segmentation based on the gray histogram of the image. In the current study, the gray histogram of the slices has single peak mode, while only a few samples have double peak mode. When the gray histogram of the image is bimodal, the local gray minimum can be used as the threshold for image segmentation.

Analysis of representative elementary volume: The representative elementary volume (REV) has the smallest size, but it contains all the information that can represent the physical property of the coal reservoir, and it is statistically significant (Vik et al., 2013; Yuan et al., 2016; Harpreet, 2017; Fang et al., 2020). The petrophysical property obtained with a size smaller than REV varies significantly, while the size larger than REV tends to be stable. By analyzing the law of variation of porosity and REV size, the size of REV can be determined.

Construction of pore network model: The Maximum Sphere algorithm is good at capturing the topological and geometrical structures of the pore and fracture, and is widely used to construct the PNM (Silin and Patzek, 2006; Lei et al., 2018). First, taking any point in the pore space as a basing point, the maximum inscribed sphere with this point as the circle center and tangent to the skeleton boundary is constantly searched. Secondly, when all the inscribed spheres are found, the other inscribed spheres contained in the inscribed spheres are removed, and the remaining inscribed spheres form the sphere set. Then, the clustering algorithm can be used to classify and summarize the maximum sphere and identify the pore and the throat. Finally, the pores can be represented as larger spheres and the throats as a set of smaller spheres (Fig. 4).

3.2 Numerical simulation method of CO₂-ECBM process

3.2.1 Mathematical model

To realize the numerical analysis of CO₂-ECBM process on laboratory scale, the mathematical model considering the parameter of pore and fracture structure should be established. In this study, the expanded Langmuir equation of competitive adsorption and the theoretical equation of adsorption, desorption and diffusion of CO₂ and CH₄ should be fully coupled in the mathematical model.

The gas diffusion in matrix is mainly controlled by its own concentration and follows Fick’s law (Clarkson and Bustin, 1999; Sun et al., 2018). Based on Fick’s first law, the continuity equation of gas adsorption and diffusion in matrix is shown in Eq. (1), namely Fick’s second law with containing gas source term S (Fick, 1855; Sun et al., 2018):

(1){\displaystyle \frac{\partial C}{\partial t}}-D\nabla C=S,

where C is the gas concentration, mol/L, which is closely related to the spatial position (x, y, z) of gas and the time (t) analyzed; and D is the gas diffusion coefficient, m²/s.

For fully coupled theoretical equations of gas adsorption and diffusion, the S in Eq. (1) can be characterized by the change of gas concentration adsorbed by the matrix over time:

(2)S=\mathrm{-}{\displaystyle \frac{\partial C_{\mathrm{a}\mathrm{d}}}{\partial t}},

where C_ad is the concentration of adsorbed gas in matrix, mol/L, which can be expressed as follows:

(3)C_{\mathrm{ad}}={\displaystyle \frac{n_{\mathrm{ad}}}{V_{\mathrm{e}}}}={\displaystyle \frac{sv/V_{\mathrm{m}}}{V_{\mathrm{e}}}}={\displaystyle \frac{sv}{V_{\mathrm{m}}{V}_{\mathrm{e}}}},

where n_ad is the gas amount adsorbed by matrix, mol; V_e is the volume of grid elements, m³; s is the pore surface area within the grid element, nm²; v is the gas amount adsorbed within the pore surface area per unit, mL; and V_m is the gas molar volume, 22.4 L/mol.

The v in Eq. (3) can be calculated by the following equation (Sun et al., 2018):

(4)v={\displaystyle \frac{V{N}_{\mathrm{solid}}{V}_{\mathrm{voxel}}{\rho }_{\mathrm{true}}}{f}},

where V is the gas volume adsorbed in coal per unit mass, m³; f is the total pore surface area, m²; N_solid is the total number of solid voxel; V_voxel is the volume per unit voxel, m³; and ρ_true is the true coal density, kg/m³.

For CO₂-ECBM process simulation on laboratory scale, it should be assumed that the gas is adsorbed only on the inner surface of the pore. When the gray value of voxel with g(x, y, z) of the position (x, y, z) meets the following conditions (Eq. (5)), the inner surface of the pore can be marked:

(5)\left\{\begin{array}{l}g(x,y,z)-g(x+\mathrm{1},y,z)=-\mathrm{1},\mathrm{o}\mathrm{r}\hfill \\ g(x,y,z)-g(x,y+\mathrm{1},z)=-\mathrm{1},\mathrm{o}\mathrm{r}\hfill \\ g(x,y,z)-g(x,y,z+\mathrm{1})=-\mathrm{1},\mathrm{o}\mathrm{r}\hfill \\ g(x-\mathrm{1},y,z)-g(x,y,z)=-\mathrm{1},\mathrm{o}\mathrm{r}\hfill \\ g(x,y-\mathrm{1},z)-g(x,y,z)=-\mathrm{1},\mathrm{o}\mathrm{r}\hfill \\ g(x,y,z-\mathrm{1})-g(x,y,z)=-\mathrm{1}\hfill \end{array}\right\},

where the gray value of pore voxels is 1, and the gray value of other voxels is 0.

For the bi-component gas of CH₄/CO₂, the C_ad in Eq. (2) can be characterized by the following extended Langmuir equation, and the diffusion theory of each component independently of other components has also been verified in the CO₂-ECBM pilot test (Shi and Durucan, 2005; Shi et al., 2008; Sun et al., 2018):

(6)V_{{\mathrm{CH}}_{\mathrm{4}}}={\displaystyle \frac{a_{{\mathrm{CH}}_{\mathrm{4}}}{b}_{{\mathrm{CH}}_{\mathrm{4}}}{P}_{{\mathrm{CH}}_{\mathrm{4}}}}{\mathrm{1}+b_{{\mathrm{CH}}_{\mathrm{4}}}{P}_{{\mathrm{CH}}_{\mathrm{4}}}+b_{{\mathrm{CO}}_{\mathrm{2}}}{P}_{{\mathrm{CO}}_{\mathrm{2}}}}},

(7)V_{{\mathrm{CO}}_{\mathrm{2}}}={\displaystyle \frac{a_{{\mathrm{CO}}_{\mathrm{2}}}{b}_{{\mathrm{CO}}_{\mathrm{2}}}{P}_{{\mathrm{CO}}_{\mathrm{2}}}}{\mathrm{1}+b_{{\mathrm{CH}}_{\mathrm{4}}}{P}_{{\mathrm{CH}}_{\mathrm{4}}}+b_{{\mathrm{CO}}_{\mathrm{2}}}{P}_{{\mathrm{CO}}_{\mathrm{2}}}}},

where V_{{\mathrm{CH}}_{\mathrm{4}}} and V_{{\mathrm{CO}}_{\mathrm{2}}} are the adsorption volume of CH₄ and CO₂, respectively, m³/kg; P_{{\mathrm{CH}}_{\mathrm{4}}} and P_{{\mathrm{CO}}_{\mathrm{2}}}are gas pressure of CH₄ and CO₂, respectively, Pa; a_{{\mathrm{CH}}_{\mathrm{4}}}and a_{{\mathrm{CO}}_{\mathrm{2}}}are the Langmuir volume of CH₄ and CO₂, respectively, m³/kg; and b_{{\mathrm{CH}}_{\mathrm{4}}}and b_{{\mathrm{CO}}_{\mathrm{2}}}are Langmuir pressure of CH₄ and CO₂, respectively, 1/Pa.

Considering the gas compressibility, the gas state equation can be characterized as follows:

(8)C={\displaystyle \frac{n}{V}}={\displaystyle \frac{P}{ZRT}}\mathrm{=}{\displaystyle \frac{P}{f(T,P)RT}},

where Z is the compression factor, which is related to temperature and pressure; T is temperature, K; and P is pressure, MPa.

Based on the above analysis, the fully coupled mathematical model of CO₂-ECBM process on laboratory scale can be expressed as follows:

(9){\displaystyle \frac{\partial P_{{\mathrm{CH}}_{\mathrm{4}}}}{\partial t}}-D_{{\mathrm{CH}}_{\mathrm{4}}}\mathrm{\Delta }{P}_{{\mathrm{CH}}_{\mathrm{4}}}=-Z_{{\mathrm{CH}}_{\mathrm{4}}}RT{\displaystyle \frac{\partial }{\partial t}}\left[{\displaystyle \frac{N_{\mathrm{solid}}{V}_{\mathrm{voxel}}{\rho }_{\mathrm{true}}s{a}_{{\mathrm{CH}}_{\mathrm{4}}}{b}_{{\mathrm{CH}}_{\mathrm{4}}}{P}_{{\mathrm{CH}}_{\mathrm{4}}}}{f{V}_{\mathrm{m}}{V}_{\mathrm{e}}(\mathrm{1}+b_{{\mathrm{CH}}_{\mathrm{4}}}{P}_{{\mathrm{CH}}_{\mathrm{4}}}+b_{{\mathrm{CO}}_{\mathrm{2}}}{P}_{{\mathrm{CO}}_{\mathrm{2}}})}}\right],

(10){\displaystyle \frac{\partial P_{{\mathrm{CO}}_{\mathrm{2}}}}{\partial t}}-D_{{\mathrm{CO}}_{\mathrm{2}}}\mathrm{\Delta }{P}_{{\mathrm{CO}}_{\mathrm{2}}}=-Z_{{\mathrm{CO}}_{\mathrm{2}}}RT{\displaystyle \frac{\partial }{\partial t}}\left[{\displaystyle \frac{N_{\mathrm{solid}}{V}_{\mathrm{voxel}}{\rho }_{\mathrm{true}}s{a}_{{\mathrm{CO}}_{\mathrm{2}}}{b}_{{\mathrm{CO}}_{\mathrm{2}}}{P}_{{\mathrm{CO}}_{\mathrm{2}}}}{f{V}_{\mathrm{m}}{V}_{\mathrm{e}}(\mathrm{1}+b_{{\mathrm{CH}}_{\mathrm{4}}}{P}_{{\mathrm{CH}}_{\mathrm{4}}}+b_{{\mathrm{CO}}_{\mathrm{2}}}{P}_{{\mathrm{CO}}_{\mathrm{2}}})}}\right].

3.2.2 Numerical schemes

For CO₂-ECBM process simulation on laboratory scale, the numerical schemes are shown in Table 2. Scheme 1 mainly represents the visual results of the CO₂-ECBM process, and scheme 2 mainly discusses the effect of the CO₂ pressure injected on the CO₂-ECBM process.

3.2.3 Loading of boundary conditions

On the laboratory scale, the loading of boundary conditions for this simulation is shown in Fig. 5. The pores in the reservoir are saturated with CH₄ (1 × 10⁻² Pa), and the CO₂ pressure is set to 0 Pa under the initial conditions. During the CO₂-ECBM process, the CH₄ pressure outside the pores, i.e., on the outer surface of the cube, remains saturated condition at 1 × 10⁻² Pa, and the CO₂ pressure is set according to the numerical schemes in Table 2 (Fig. 5).

3.2.4 Numerical parameters

The simulation parameters required for this CO₂-ECBM process simulation were all derived from correlation analysis experiment, and the gas attribute parameters required for this simulation were set according to the attribute parameters of CH₄ and CO₂ (Table 3).

3.2.5 Development of numerical software

The COMSOL, an advanced multiphysical field finite element software, has a wide simulation capability and strong post-processing capability to analyze and solve the fully coupled mathematical equations, but no corresponding processing interface can be added, and the data optimization capability of post-processing is low. The MATLAB software can be used to extract the geological model required for numerical simulation, i.e., the PNM, and the simulation results can be optimized in MATLAB software after COMSOL processing (Fig. 6). Therefore, the comprehensive application of COMSOL and MATLAB software is very necessary.

For CO₂-ECBM process simulation, the COMSOL software can create a graphical user interface (GUI) to realize the construction of COMSOL and MATLAB simulation system, and the GUI can form an independent software package based on typical calculation methods. First, the MATLAB script can be called in the GUI of the COMSOL software to build the geological model, and the grid partitioning and optimization of the geological model can be realized. Secondly, based on the MATLAB script, the geological model grid can be recognized in the GUI of COMSOL software. Then, based on the built-in PDE function of COMSOL software, parameters, variables and boundary conditions can be set in the GUI, and the numerical simulation can be successfully completed. Finally, the script function and language in MATLAB wrote by ourselves can be called to perform 3D visualization and data optimization after COMSOL processing (Fig. 6). Based on the MATLAB script, the interaction and sharing of data between COMSOL and MATLAB software can be realized.

4 Results

4.1 3D reconstruction of pore and fracture structure

The 3D reconstruction of pore and fracture structure for typical 2D CT slices can extract the geological carrier, i.e., the PNM, needed for this numerical simulation. The preprocessing of 2D CT slices, threshold selection and image segmentation, analysis of REV, and construction of PNM are the key works, which are as follows:

4.1.1 Preprocessing of 2D CT slices

From Fig. 7(a), it can be seen that many noise points appeared in the original 2D CT slices. The Median Filtering algorithm was used to process the 2D CT slices, and the results can be seen in Fig. 7(b). After filtering, the transition between skeleton and pore is smooth and natural, and the unrealistic outlier is no longer present in the image (Fig. 7). During the actual filtering, the filtered slices are often compared with the original slices to see whether some of the pores in coal have been erased.

4.1.2 Threshold selection and image segmentation

First, the frequency distribution histogram of gray can be calculated according to the gray distribution of the image. The gray range (H) is between 159 and 175 for BF sample (Fig. 8(b)). Second, a certain threshold (T_V) is selected based on the gray range (H) to divide the gray into two parts:H_{\mathrm{0}\to T_V}, H_{T_V}{}_{+\mathrm{1}\to \mathrm{255}}, to calculate the variance a\left(T_V\right)=Var\left(H_{\mathrm{159}}{}_{\to T_V}\right), b\left(T_V\right)=Var\left(H_T{{}_{+\mathrm{1}}}_{\to \mathrm{175}}\right), and to calculate the difference of variance D(T_V) = |a(T_V)-b(T_V)|, where T_V∈H. Then, the threshold value T_V = find[max(D)] is defined to find the value of T_V in order to maximize D. Finally, after selecting T_V, the slice processed by T_V is often compared with the original 2D CT slice to detect whether some of the pores are deleted, and T_V can be partially fine-tuned (Fig. 8(b)). Based on this, the pores, organic matter and inorganic minerals of the BF sample can be extracted and reconstructed in 3D (Fig. 8(c)).

4.1.3 Extraction of pore network model

Analysis of the relationship between the porosity and the size of REV shows that when REV is larger than 500 × 500 × 500 voxel, porosity is relatively stable with the change of the REV size. Therefore, the 500 × 500 × 500 voxel can be selected to represent the size of REV (Fig. 9(a)). Based on the application of the Maximum Sphere algorithm (Fig. 9(b)), the pore and throat can be extracted from the PNM, separately (Fig. 9(c)).

4.2 Numerical result of the CO₂-ECBM process

4.2.1 Pre-treatment of the geological model

Based on the MATLAB software, the interconnected pores and throats can be extracted from the PNM, and the STL file can be imported into COMSOL software for the CO₂-ECBM process simulation, so as to bridge from the geometric model to the numerical simulation (Fig. 10). Given the computer memory capacity requirements of COMSOL software, when the sample size analyzed are larger than 60 × 60 × 60 voxel, the simulation process in COMSOL software will overflow due to the lack of computer memory. Therefore, a mesh size of 60 × 60 × 60 voxel was chosen for the CO₂-ECBM numerical analysis (Fig. 10(a)). Due to the complexity of pore and fracture structure, errors may occur in the meshing of geological model. Therefore, manual repair and debugging should be performed for the faulty geological model (Fig. 10(b)). Through continuous repair and debugging, the COMSOL software can generate the error-free tetrahedral mesh required for this numerical simulation (Fig. 10(c)).

4.2.2 Visualization of numerical results

Based on the COMSOL software, the post-processing and analysis of numerical results in laboratory scale were carried out, and the distribution of the pressure field of CO₂ and CH₄ in 3D, 2D and 1D can be obtained (Figs. 11–13).

Figure 11 shows the 3D distribution of CO₂ and CH₄ pressure during the CO₂-ECBM process. For CH₄, with the increase of CO₂ displacement time for CH₄, the CH₄ pressure gradually decreases at both the edge and the center of the simulation model (Fig. 11(a)). For CO₂, the CO₂ pressure gradually increases from the edge to the center of the simulation model with increasing CO₂ displacement time for CH₄ (Fig. 11(b)). During the whole period of CO₂ displacement for CH₄, the pressure drop was the same (ΔP= 0.01 Pa), but the 3D distribution of gas pressure at different times and positions was greatly different (Fig. 11), which lies in the differences of radius, shape and connectivity of pore and throat. When CH₄ is displaced by CO₂, the area where CH₄ pressure gradually decreases are also the area where CO₂ pressure gradually increases (Fig. 11).

The pressure field in different slices was quantitatively analyzed. The CO₂ gradually disinfects the center from the edge of the model. Therefore, analysis should be performed with the positions of 85 μm (5th slice), 90 μm (10th slice), 95 μm (15th slice), 100 μm (20th slice), 105 μm (25th slice), and 110 μm (30th slice) and with the analysis time of 15th second, 30th second, 45th second and 60th second, respectively (Fig. 12). In the original 2D CT slices, the black area represents the matrix and the white area represents the pore and fracture, which are the geological carrier for CO₂-ECBM process simulation.

Figure 12 shows the 2D distribution of CO₂ and CH₄ pressure during the CO₂-ECBM process. For CH₄, as the time of CO₂ displacement increases, the gas pressure decreases from the center to the edge in all slices, and the total CH₄ pressure in reservoir also gradually decreases with time. At the same time, the CH₄ pressure in the central slices, such as the 10th, 15th, 20th and 25th slice, was relatively high, while that in the edge slices, such as the 5th and 30th slice, was relatively low (Fig. 12). For CO₂, the pressure distribution is opposite to that of CH₄. As the time of CO₂ displacement CH₄ increase, the CO₂ pressure increases from the edge to the center in all slices, and the total CO₂ pressure in reservoir also gradually increases with time. At the same time, the CO₂ pressure in the central slices, such as the 10th, 15th, 20th and 25th slice, was relatively low, while that in the edge slices, such as the 5th and 30th slice, was relatively high (Fig. 12).

In this study, three monitoring points of A(82, 128, 109), B(87, 120, 107) and C(96, 115, 104) were selected to quantitatively investigate the distribution law of CH₄ and CO₂ pressure at different positions over time (Fig. 13), and points A to C represent the position from the edge of the reservoir to the center.

Figure 13 shows the 1D distribution of CO₂ and CH₄ pressure during the CO₂-ECBM process. For CH₄, the CH₄ pressure at different positions gradually decreases with increasing CO₂ displacement time of CH₄. Here, the closer to the center of the model, the higher the CH₄ pressure, and the closer to the edge, the lower the CH₄ pressure (Fig. 13(a)). In terms of CO₂, the CO₂ pressure at different positions gradually increases with the increase of CO₂ displacement CH₄ time. Here, the closer it is to the center of the model, the lower the CO₂ pressure is, while the closer it is to the edge of the model, the higher the CO₂ pressure is (Fig. 13(b)). The gas pressure changes rapidly in the early stage (0‒20 s) of CO₂ displacement CH₄ and slowly in the later stage (>20 s) (Fig. 13).

4.2.3 Effect of CO₂ pressure injected on the CO₂-ECBM process

The saturation pressure of CH₄ was kept at 1 × 10⁻² Pa, and the effect of gas pressure injected on CO₂-ECBM process was analyzed by changing the CO₂ pressure injected according to the scheme 2. This part focuses on the analysis of the distribution law of CO₂.

Figure 14 shows the 3D distribution of the CO₂ pressure field when the CO₂ pressure injected is changed during the CO₂-ECBM process. At the same CO₂ injection time, the CO₂ pressure shows a trend of gradual increase with the increase of CO₂ pressure injected, and the differences of CO₂ pressure in reservoir are relatively clear at each CO₂ pressure injected. At different CO₂ pressures injected, the CO₂ pressure changes in reservoir are all large, and the CO₂ pressure changes in the center of the slice are relatively small, while those at the edge of the slice are relatively large (Fig. 14).

To further analyze the distribution of CO₂ pressure field in different slices under different pressures injected, the position of 95 μm on the X-axis was taken as an example for analysis, and the analysis time was 1st second, 10th second, 20th second, 40th second and 60th second (Fig. 15).

Figure 15 shows the 2D distribution of the CO₂ pressure field when the CO₂ pressure injected is changed during the CO₂-ECBM process. At the same gas injection time, the CO₂ pressure shows a trend of gradual increase with the increase of CO₂ pressure injected. For the same gas pressure injected, the CO₂ pressure gradually increases from the slice edge to the center as the increase of CO₂ injection time (Fig. 15).

Monitoring point B(87, 120, 107) was selected to quantitatively investigate the gas pressure changes in reservoir under different CO₂ pressures injected in the CO₂-ECBM process (Fig. 16).

Figure 16 shows the distribution of CO₂ pressure at point B at different CO₂ pressures injected. At different CO₂ pressures injected, the CO₂ pressure in pore gradually increases with increasing CO₂ injection time. At the same time, the higher the gas pressure injected, the faster the increasing rate of CO₂ pressure, which indicates that the diffusion and adsorption rates of CO₂ increase with the increase of gas pressure injected (Fig. 16). The higher the CO₂ pressure, the sooner the CO₂ pressure in pore reaches the steady-state, i.e., T_S5>T_S4>T_S3>T_S2>T_S1 (Fig. 16).

5 Discussion

5.1 Dynamic characteristics of the CO₂-ECBM continuous process

Figure 17 is a simple diagram of the continuous CO₂-ECBM process. Based on this figure and the previous discussion on CO₂-ECBM theory, the continuous process with adsorption, desorption, diffusion and percolation of CH₄ and CO₂ during the CO₂-ECBM process can be analyzed.

As for the CO₂, the CO₂ injected mainly migrates along the macro- and micro-fractures with continuous flow into the matrix. The CO₂ injected first replaces the CH₄ adsorbed by covering the inner surface of the macro- and the meso-pores, so that a single molecular layer of CO₂ is adsorbed. Then it migrates to the micro-pores by Fick’s diffusion, slip flow and surface diffusion. Moreover, CO₂ replaces CH₄ adsorbed by volumetric filling of the micro-pores and forms a multi-molecular layer adsorption of CO₂ (Fig. 17). In the case of CH₄, the CH₄ in reservoir is saturated in the initial state, and CH₄ molecules maintain a dynamic equilibrium state of adsorption and desorption behavior in matrix. The injection of CO₂ disrupts this equilibrium of CH₄ in matrix. Since the matrix has a higher adsorption capacity for CO₂ than for CH₄, the matrix has priority in adsorption of CO₂ over desorption of CH₄ in the competition between CO₂ and CH₄, which can complete the process of CO₂ replacement of CH₄ (Fig. 17).

The CH₄ and CO₂ adsorbed on the inner surface of the matrix complete in the adsorption position. In general, the adsorption capacity of matrix is twice as high for CO₂ as for CH₄ (Fig. 18). The CH₄ desorbed diffuses from the surface of matrix into the micro-pores under the action of the concentration gradient, and the process follows Fick’s law. Then, it moves from the pore to the fracture and then to the wellbore through seepage under the action of the pressure gradient, which follows Darcy’s Law (Fig. 19). Based on the analysis of Dalton’s law, in the CO₂-ECBM process, although the total gas pressure in reservoir remains unchanged, the partial pressure of CO₂ in reservoir increases with time, while the partial pressure of CH₄ decreases with time. The CO₂ injected in reservoir gradually replaces the CH₄ in reservoir through competitive adsorption, which improves the CH₄ recovery while storing CO₂ (Fig. 19).

5.2 Influential effect of pore and fracture structure on continuous CO₂-ECBM process

The coal reservoir can be abstracted as a dual-pore medium consisting of pores and fractures. The term “Double-Pore” refers to the pore system and the fracture system (Fig. 20). Pores are mainly evolved by physical and chemical processes during coal formation. Fractures are mainly formed by the late stress effect. In reservoir, the multilevel pore and fracture structure network formed by interconnected pores and fractures will affect the adsorption, desorption, diffusion and seepage of fluid in reservoir. The fluid is in dynamic equilibrium in reservoir, the pore system in matrix is mainly the adsorption space of fluid, and the migration channel of the fluid is mainly the fracture system.

The pore and fracture systems in reservoir form the multilevel network structure, and this is the main occurrence space and migration channel of CH₄ and CO₂ in the CO₂-ECBM process. The influential effect on the continuous process during the CO₂-ECBM process can be summarized as follows: (1) micro-pores and meso-pores in reservoir are the main occurrence place of CH₄ and CO₂ during the CO₂-ECBM process. (2) The migration and production path of CH₄ is from micro-pores, meso-pores, macro-pores, micro-fractures, endogenous fractures, macro-fractures, to fracturing fractures during CO₂-ECBM process. (3) The migration path of CO₂ is exactly opposite to that of CH₄ during CO₂-ECBM process. (4) The output of CH₄ passes through three flow levels, i.e., from pores, natural fractures, fracturing fractures, to wellbore. (5) At the macro level, CO₂ injection also passes through three flow levels, which is from wellbores, fracturing fractures, natural fractures, to pores.

6 Conclusions

1) The perfect coupling among the X-ray CT technology, the MATLAB, with the COMSOL software cannot only realize the 3D reconstruction of pore network model, but also complete the CO₂-ECBM process simulation, and further investigate the influential effect of PNM on such process. The Median filtering algorithm enables the de-noising of the original 2D CT slices, the threshold selection and image segmentation are realized based on the image segmentation method, and the PNM can be constructed based on the Maximum Sphere algorithm.

2) The CO₂-ECBM mathematical model considering PNM parameter completely couples the expanded Langmuir equation of competitive adsorption and the theoretical equation of adsorption, desorption and diffusion of CO₂/CH₄. At the same injection time, CH₄ pressure gradually decrease with the increase of CO₂ pressure injected, but the difference of CH₄ pressure at each CO₂ pressure is not obvious. At different CO₂ pressures, the CH₄ pressure changed a lot at the slice center, while it changed a little at the slice edge. The law of change of CO₂ pressure is opposite to that of CH₄ pressure.

3) The CO₂ injected is transported into the matrix along the macro- and micro-fractures with continuous flow. The CO₂ injected first replaces the CH₄ adsorbed by covering the inner surface of the macro-pores and the meso-pores to form the adsorption of CO₂ in a single molecular layer. Then they migrate to the micro-pores by Fick’s diffusion, slip flow and surface diffusion. Moreover, the CO₂ replaces the CH₄ adsorbed in the micro-pores by volumetric filling and forms the multi-molecular layer adsorption of CO₂. The migration path of CH₄ is opposite to that of CO₂.

References

[1] Cai Y D, Liu D M, Yao Y B, Li J G, Qiu Y K (2011). Geological controls on prediction of coalbed methane of No. 3 coal seam in Southern Qinshui Basin, North China. Int J Coal Geol, 88(2–3): 101–112

[2] Cai T T, Feng Z C, Zhou D (2018). Multi-scale characteristics of coal structure by X-ray computed tomography (X-ray CT), scanning electron microscope (SEM) and mercury intrusion porosimetry (MIP). AIP Adv, 8(2): 025324

[3] Cervik J (1967). Behavior of coal-gas reservoirs. In: SPE Eastern Regional Meeting

[4] Cheng G X, Jiang B, Li M, Li F L, Song Y (2020). Effects of pore-fracture structure of ductile tectonically deformed coals on their permeability: an experimental study based on raw coal cores. J Petrol Sci Eng, 193: 107371

[5] Clarkson C R, Bustin R M (1999). The effect of pore structure and gas pressure upon the transport properties of coal: a laboratory and modeling study. Fuel, 78(11): 1345–1362

[6] Connell L D, Mazumder S, Sander R, Camilleri M, Pan Z J, Heryanto D (2016). Laboratory characterization of coal matrix shrinkage, cleat compressibility and the geo-mechanical properties determining reservoir permeability. Fuel, 165: 499–512

[7] Fan N, Wang J R, Deng C B, Fan Y P, Wang T T, Guo X Y (2020). Quantitative characterization of coal microstructure and visualization seepage of macropores using CT-based 3D reconstruction. J Nat Gas Sci Eng, 81: 103384

[8] Fang H H, Sang S X, Wang J L, Liu S Q, Ju W (2017). Simulation of paleotectonic stress fields and distribution prediction of tectonic fractures at the Hudi Coal Mine, Qinshui Basin. Acta Geol Sin-Engl, 91(6): 2007–2023

[9] Fang H H, Sang S X, Liu S Q (2019a). Establishment of dynamic permeability model of coal reservoir and its numerical simulation during the CO₂-ECBM process. J Petrol Sci Eng, 179: 885–898

[10] Fang H H, Sang S X, Liu S Q, Du Y (2019b). Methodology of three-dimensional visualization and quantitative characterization of nanopores in coal by using FIB-SEM and its application with anthracite in Qinshui Basin. J Petrol Sci Eng, 182: 106285

[11] Fang H H, Sang S X, Liu S Q (2019c). The coupling mechanism of the thermal-hydraulic-mechanical fields in CH₄-bearing coal and its application in the CO₂-enhanced coalbed methane recovery. J Petrol Sci Eng, 181: 106177

[12] Fang H H, Sang S X, Liu S Q (2019d). Numerical simulation of enhancing coalbed methane recovery by injecting CO₂ with heat injection. Petrol Sci, 16(1): 32–43

[13] Fang H H, Sang S X, Liu S Q, Liu S P (2019e). Experimental simulation of replacing and displacing CH₄ by injecting supercritical CO₂ and its geological significance. Int J Greenh Gas Control, 81: 115–125

[14] Fang H H, Sang S X, Liu S Q (2020). Three-dimensional spatial structure of the macro-pores and flow simulation in anthracite coal based on X-ray mu-CT scanning data. Petrol Sci, 17(5): 1221–1236

[15] Fick A (1855). On liquid diffusion. Philos Mag, 10(63): 30–39

[16] Huang S P, Liu D M, Cai Y D, Gan Q (2019). In situ stress distribution and its impact on CBM reservoir properties in the Zhengzhuang area, southern Qinshui Basin, North China. J Nat Gas Sci Eng, 61: 83–96

[17] Lei J, Pan B Z, Zhang L H (2018). Advance of microscopic flow simulation based on digital cores and pore network. Prog Geophys, 33(2): 653–660 (in Chinese)

[18] Li Y, Zhang C, Tang D, Gan Q, Niu X, Wang K, Shen R (2017). Coal pore size distributions controlled by the coalification process: an experimental study of coals from the Junggar, Ordos, and Qinshui basins in China. Fuel, 206: 352–363

[19] Li Z W, Zhang G Y (2019). Fracture segmentation method based on contour evolution and gradient direction consistency in sequence of coal rock CT images. Math Probl Eng, 1: 1–82980747

[20] Li Y, Wang Y, Wang J, Pan Z (2020a). Variation in permeability during CO₂-CH₄ displacement in coal seams: part I-Experimental insights. Fuel, 263: 116666

[21] Li Y, Yang J, Pan Z, Tong W (2020b). Nanoscale pore structure and mechanical property analysis of coal: an insight combining AFM and SEM images. Fuel, 260: 116352

[22] Liu S Q, Sang S X, Wang G, Ma J S, Wang T, Wang W, Du Y, Wang T (2017). FIB-SEM and X-ray CT characterization of interconnected pores in high-rank coal formed from regional metamorphism. J Petrol Sci Eng, 148: 21–31

[23] Ni G H, Li S, Rahman S, Xun M, Wang H, Xu Y H, Xie H C (2020a). Effect of nitric acid on the pore structure and fractal characteristics of coal based on the low-temperature nitrogen adsorption method. Powder Technol, 367: 506–516

[24] Ni X M, Zhao Z, Wang Y B, Wang L (2020b). Optimisation and application of well types for ground development of coalbed methane from No. 3 coal seam in shizhuang south block in Qinshui Basin, Shanxi Province, China. J Petrol Sci Eng, 193: 107453

[25] Nie B S, Lun J Y, Wang K D, Shen J S (2020). Three-dimensional characterization of open and closed coal nanopores based on a multi-scale analysis including CO₂ adsorption, mercury intrusion, low-temperature nitrogen adsorption, and small-angle X-ray scattering. Energy Sci Eng, 8(6): 2086–2099

[26] Shi J Q, Durucan S (2005). A model for changes in coalbed permeability during primary and enhanced methane recovery. SPE Reservoir Eval Eng, 8(04): 291–299

[27] Shi J Q, Mazumder S, Wolf K H, Durucan S (2008). Competitive methane desorption by supercritical CO₂ injection in coal. Transp Porous Media, 75(1): 35–54

[28] Silin D, Patzek T (2006). Pore space morphology analysis using maximal inscribed spheres. Physica A, 371(2): 336–360

[29] Singh H (2017). Representative Elementary Volume (REV) in spatio-temporal domain: a method to find REV for dynamic pores. J Earth Sci, 28(2): 391–403

[30] Sun Y F, Zhao Y X, Yuan L (2018). CO₂-ECBM in coal nanostructure: modelling and simulation. J Nat Gas Sci Eng, 54: 202–215

[31] Vik B, Bastesen E, Skauge A (2013). Evaluation of representative elementary volume for a vuggy carbonate rock-Part: porosity, permeability, and dispersivity. J Petrol Sci Eng, 112(3): 36–47

[32] Wang G, Wang K, Jiang Y J, Wang S Q (2018a). Reservoir permeability evolution during the process of CO₂-enhanced coalbed methane recovery. Energies, 11(11): 2996

[33] Wang G, Wang K, Wang S G, Elsworth D, Jiang Y J (2018b). An improved permeability evolution model and its application in fractured sorbing media. J Nat Gas Sci Eng, 56: 222–232

[34] Wang X M, Wang X M, Pan Z J, Yin X B, Chai P C, Pan S D, Yang Q (2019). Abundance and distribution pattern of rare earth elements and yttrium in vitrain band of high-rank coal from the Qinshui Basin, northern China. Fuel, 248: 93–103

[35] Wang G, Han D Y, Jiang C H, Zhang Z Y (2020a). Seepage characteristics of fracture and dead-end pore structure in coal at micro- and meso-scales. Fuel, 266: 117058

[36] Wang M F, Wang J J, Tao S, Tang D Z, Wang C C, Yi J (2020b). Quantitative characterization of void and demineralization effect in coal based on dual-resolution X-ray computed tomography. Fuel, 267: 116836

[37] Wang H, Yao Y B, Huang C C, Liu D M, Cai Y D (2021). Fault development characteristics and their effects on current gas content and productivity of No.3 coal seam in the Zhengzhuang Field, southern Qinshui Basin, north China. Energ Fuels, 35(3): 2268–2281

[38] Yuan C, Chareyre B, Darve F (2016). Pore-scale simulations of drainage in granular materials: finite size effects and the representative elementary volume. Adv Water Resour, 95: 109–124

[39] Zhang Y X, Yuan Q, Huang D, Kong S F, Zhang J, Wang X F, Lu C Y, Shi Z B, Zhang X Y, Sun Y L, Wang Z F, Shao L Y, Zhu J H, Li W J (2018). Direct observations of fine primary particles from residential coal burning: insights into their morphology, composition, and hygroscopicity. J Geophys Res Atmos, 123(22): 12964–12979

[40] Zhang K, Sang S X, Zhou X Z, Liu C J, Ma M Y, Niu Q H (2021). Influence of supercritical CO₂-H₂O-caprock interactions on the sealing capability of deep coal seam caprocks related to CO₂ geological storage: a case study of the silty mudstone caprock of coal seam No. 3 in the Qinshui Basin, China. Int J Greenh Gas Control, 106: 103282

[41] Zheng S J, Yao Y B, Elsworth D, Liu D M, Cai Y D (2020). Dynamic fluid interactions during CO₂-ECBM and CO₂ sequestration in coal seams. Part 2: CO₂-H₂O wettability. Fuel, 279: 118560

[42] Zhu Q Z, Wang X L, Zuo Y Q, Pan J N, Ju Y W, Su X F, Yu K (2020). Numerical simulation of matrix swelling and its effects on fracture structure and permeability for a high-rank coal based on X-ray micro-CT image processing techniques. Energ Fuels, 34(9): 10801–10809

RIGHTS & PERMISSIONS

Higher Education Press