Three-dimensional spatial structure of the macro-pores and flow simulation in anthracite coal based on X-ray μ-CT scanning data

The three-dimensional (3D) structures of pores directly affect the CH4 flow. Therefore, it is very important to analyze the 3D spatial structure of pores and to simulate the CH4 flow with the connected pores as the carrier. The result shows that the equivalent radius of pores and throats are 1–16 μm and 1.03–8.9 μm, respectively, and the throat length is 3.28–231.25 μm. The coordination number of pores concentrates around three, and the intersection point between the connectivity function and the X-axis is 3–4 μm, which indicate the macro-pores have good connectivity. During the single-channel flow, the pressure decreases along the direction of CH4 flow, and the flow velocity of CH4 decreases from the pore center to the wall. Under the dual-channel and the multi-channel flows, the pressure also decreases along the CH4 flow direction, while the velocity increases. The mean flow pressure gradually decreases with the increase of the distance from the inlet slice. The change of mean flow pressure is relatively stable in the direction horizontal to the bedding plane, while it is relatively large in the direction perpendicular to the bedding plane. The mean flow velocity in the direction horizontal to the bedding plane (Y-axis) is the largest, followed by that in the direction horizontal to the bedding plane (X-axis), and the mean flow velocity in the direction perpendicular to the bedding plane is the smallest.


Introduction
The coal reservoir is a complex porous medium, which contains matrixes, pores, and minerals (Ni et al. 2017;Liu et al. 2019). The pore structure directly affects the physical property of the coal reservoir (Fu et al. 2017;Su et al. 2019) and further affects the CH 4 flow (Naveen et al. 2018;Zhou et al. 2018;Liu et al. 2019). Among the multi-scale pores developed in coal reservoir, most of the macro-pores are connected with nanometer-scale pores and millimeter-scale fractures (Ni et al. 2017), which play important roles in the CH 4 flow. The analysis of the 3D spatial structure of the macro-pores and the simulation of CH 4 flow with the interconnected pores as the carrier will lay foundations for the study of CH 4 flow on the microscopic scales.
In the study of pore structure, traditional methods such as the mercury intrusion porosimetry and the low temperature liquid nitrogen adsorption analysis (Guo et al. 2019;Wang et al. 2018) often destroy the original pore structure during the sample preparation process, so the real 3D pore structures cannot to be characterized. Compared with the Edited by Jie Hao Hui-Huang Fang and Shu-Xun Sang have contributed equally to this work. traditional methods, the focused ion beam scanning electron microscopy (FIB-SEM) and the X-ray computed tomography (i.e., Nano-CT and μ-CT) are new methods to study the pore structure (Yao et al. 2009a;Ma et al. 2014;Shaina et al. 2016;Xiong et al. 2016;Wu et al. 2019a). The FIB-SEM and Nano-CT can only characterize the nanopore structure, and their application and promotion are further limited by the high resolution, long experimental period and strict requirements of sample preparation (Ni et al. 2017). Due to the limited sample size observed, the analysis of the pore structure of coal reservoir with strong heterogeneity is not very representative (Ni et al. 2017).
Compared with the FIB-SEM, the X-ray micro-focus computed tomography (μ-CT) is a non-destructive scanning technology (Sleutel et al. 2008;Hamamoto et al. 2016), which can not only identify the pores without destroying the original structure, but also solve the contradiction between sample size and scanning resolution (Ni et al. 2017). Based on the X-ray μ-CT, the evolution of coal structure, permeability and heterogeneity (Wolf et al. 2008;Heriawan and Koike 2015), the analysis of coal damage (Viljoen et al. 2015) and the spatial distribution of pores, fractures and minerals (Yao et al. 2009a, b) have been extensively studied. The X-ray μ-CT technology can also be combined with the specific algorithms and the visualization technology to achieve the 3D visualization of coal structure (Yao et al. 2009a, b;Ni et al. 2017). The application of the X-ray μ-CT technology is shown to be effective to study the 3D pore structure (Sun et al. 2016;Ni et al. 2017).
The visualization analysis of the gas flow with the interconnected pores as the carrier and the discussion of the pressure field and velocity field are of important engineering value for the development of low-permeability reservoirs (Zhou et al. 2016;Ni et al. 2017). Based on the 2D CT slices of the sandstone, Ju et al. (2014) have simulated the gas flow with the interconnected pores in sandstone as the carrier through using the Lattice Boltzmann method, but this method will take a long time and consume a lot of manpower Yang et al. 2018). With the development of the 3D visualization technology and the multi-field coupling finite element method, Bird et al. (2014) have used the AVIZO software and the COMSOL software to simulate the gas flow with the interconnected pores in carbonate as the carrier. At present, the visualization technology of gas flow at micro-scale is mainly applied to sandstone and carbonate reservoirs Yang et al. 2018). Moreover, the simulation and visualization of the gas flow within the pore scale in coal reservoir are less, especially the CH 4 flow. Therefore, it is necessary to explore new ways to simulate and visualize the CH 4 flow. Recent studies have shown that the Pore Network Model extracted from the 2D CT slices is composed of the pores, throats and their coordination number (Youssef et al. 2007;Zhao et al. 2018;Wu et al. 2019b), which not only can characterize the complexity of the pore space, but also can be used as the carrier for the CH 4 flow simulation on the macro-pore scale (Ni et al. 2017;Blunt 2001;Wu et al. 2019b).
The objective of this study is to realize the docking technology of the programming software and the simulation software through using the STL files and then to realize the analysis of velocity field and pressure field of CH 4 flow under various conditions. In this study, the coal sample collected from the southern Qinshui basin in China was firstly scanned by the X-ray μ-CT to obtain the 2D CT slices. Secondly, the original 2D CT slices were reconstructed and divided to obtain the equivalent pore network model. Then, a mathematical model of CH 4 flow simulation with macropores as the carrier was established. Finally, the simulation and its visualization of CH 4 flow were realized. The innovations of this study are as follows: (1) a numerical model for coal is built and its equivalent pore network model is extracted based on the X-ray μ-CT data.
(2) The distribution of the pressure field and velocity field in signal-, dual-, and multi-channel flows along different directions are revealed. The analysis of the 3D spatial structure and flow simulation of the macro-pores based on the perfect docking of the X-ray μ-CT technology, the image processing technology and the COMSOL can provide new approaches for exploring the CH 4 flow in the macro-pore scale.

Sample
One anthracite coal from the Bofang Mine (BF) in the Jincheng area was collected from the No.3 coal seam of the Shanxi formation in the southern Qinshui basin of China ( Fig. 1; Liu et al. 2017;Fang et al. 2017). The collection, retention, and preparation of the coal sample were conducted in line with the relevant Chinese and international standards, such as the GB/T 19222-2003, the GB/T 16773-2008, and the ISO 7404-2:1985(Liu et al. 2017Fang et al. 2019). Sample was sealed with absorbent paper to avoid oxidation.
The key properties of the sample are shown in Table 1.

X-ray μ-CT scanning
In this study, the X-ray micro-focus CT scanner (i.e., Xradia 510 Versa, the Carl Zeiss Foundation Group) was used for sample scanning, which is a digital core analysis system with resolution up to 0.5 µm. The experimental coal sample was cylindrical in shape of Φ 2.0 × 2.0 mm, which is placed in the center of the CT system workbench. The energy of the X-ray is 180 kV/15 W, the scanning speed is 175 mm/s, and the image resolution is 0.5 µm. It is necessary to dry the sample before scanning. In this study, the coal sample can be cut into 1013 vertical slices with each thickness of 0.5 µm. The 1013 slices with 988 × 988 pixels can be obtained. Based on the scanning analysis of the tested sample, some typical 2D CT slices can be seen in Fig. 2 the matrixes, pores, and minerals, respectively, inside the circle (Fig. 2).

Image processing
The original 2D CT slices of the BF sample should be reconstructed and segmented to obtain the equivalent Pore Network model, which includes the pre-processing of 2D CT slices, the image segmentation, the selection of the representative elementary volume and the extraction of the equivalent pore network model.

Pre-processing of 2D CT slices
It is generally known that the higher the image processing accuracy is, the better the 3D reconstruction effect of the pore structure will be (Zhou et al. 2016). The original 2D CT slices obtained based on the X-ray μ-CT scanning are more or less affected by noise, so there are scatter dots noise in the CT slices ( Fig. 3a), which will have an adverse impact on image segmentation in the following operation; therefore, the original 2D CT slices need to be de-noised. The median filter can not only protect the integrity of the pores but also smooth the transition between the pores and the matrixes (Li and Zhang 2019), which can de-noise the original 2D CT slices. In this study, the original 2D CT slices are processed by the 3D Median filter with 18 neighborhoods, and the results are shown in Fig. 3b.

Image segmentation
The conversion from the 2D CT slices into the segmented images is the primary purpose of threshold segmentation, which can separate the pore and matrix from the coal, and further separate the organic matter and mineral from the matrix. The porosity (φ) of the scanned sample can be obtained through the laboratory experiment, and the segmentation threshold obtained based on the porosity of the digital core can be used to segment the image (Fang et al. 2018). First, a rough threshold range can be obtained by using the porosity method. Then, multiple pore extraction should be carried out to compare the distribution of pores in the original 2D CT Slices. Finally, the threshold range should be further adjusted to get the optimal threshold. The equation is as follows: where K is the segmentation threshold; φ is the porosity; I max and I min are the maximum and minimum grayscale values of the image; and P(i) is the voxel number with grayscale value of i. Voxels with grayscale values below the segmentation threshold represent pores, while those with grayscale values greater than the segmentation threshold represent matrixes and minerals. The 3D reconstruction of coal sample, and the pores, matrixes and minerals can be separated from the reconstructed coal based on the segmentation threshold ( Fig. 4).

Representative elementary volume
The representative elementary volume (REV) is the smallest unit, which can effectively characterize the macroscopic physical property of coal (Vik et al. 2013;Yuan et al. 2015). Selecting an appropriate REV is the key to analyzing the 3D spatial structure of pore, which can effectively reduce  preet 2017). After many calculations, the results demonstrate that the porosity of coal has nothing to do with the sample size when the sample size is over 500 × 500 × 500 voxels (Fig. 5). Therefore, the 500 × 500 × 500 voxels are chosen to characterize the sample size, and the actual sample size is 250 × 250 × 250 μm 3 .

Equivalent pore network model
The equivalent pore network model is the carrier of CH 4 flow. Understanding the number, radius and length of throats, as well as the relative position and connectivity of pores and throats are the prerequisite for analyzing the geometrical and topological structures of the pore network model, as well as the basis for further simulating and visualizing the CH 4 flow. The Maximum Ball method can be used to establish the pore network model, and the specific principles are as follows (Silin and Patzek 2006;Lei et al. 2018): the Maximum Ball method takes any point in pore space as the reference point (red part in Fig. 6) and constantly looks for the maximum inscribed sphere with this point as the circle center and tangent to the skeleton boundary (red circle in Fig. 6). After all the inscribed spheres are found, the spheres contained in other inscribed spheres will be removed, and the remaining spheres will constitute the largest ball set. The clustering algorithm is introduced to classify and merge the maximum sphere and identify pores and throats. Pores can be represented by the larger spheres  Figure 7 denotes the pore network model of the BF sample, and the cylinder and sphere represent the throats and the pores, respectively.

Flow simulation
In order to simulate and visualize the CH 4 flow after extracting the Pore Network model, the following steps such as the pre-processing of pore grid, the controlling equations and the basic setting of the boundary conditions need to be established. Although a computer with 24 GB of the random-access memory has been used, the 3D mesh is often interrupted by a memory overflow in COMSOL when the sample size exceeds 60 × 60 × 60 voxels. Therefore, the sample sizes have to be reduced to 60 × 60 × 60 voxels (i.e., 30 × 30 × 30 μm 3 ) for CH 4 flow simulation as derived from REV interception in this study (Fig. 8).

Pre-processing of pore grid
The interconnected pores and throats are the major channels for CH 4 flow. Therefore, it is necessary to extract the interconnected pores and throats in the pore network model. The interconnected pores and throats which can represent the real connected environment of the BF sample can be extracted based on the corresponding programming in the MATLAB software (Fig. 8a), and the output STL file can be imported into the COMSOL (www.comso l.com) to simulate the CH 4 flow. Due to the complex pore structure, the errors or warnings often occur when the Pore Network model is meshed in the COMSOL. Therefore, the Pore Network model should be manually debugged in the COMSOL, which includes the restoration of the notches and sharp corners, and the elimination of the intersect, aperture, small hole, chamfer, interface, and coincident edge (Zhou et al. 2016;Ni et al. 2017). Through continuous debugging, the flawless tetrahedral mesh for numerical simulation can be generated in COMSOL (Fig. 8b).

Controlling equations of flow simulation
In this study, it is assumed that the fluid is incompressible and the flow is laminar in the pore network model. The Navier-Stokes model and the continuity equation can be applied to the flow simulation, which are as follows: where ρ, v, P, and μ denote the density, velocity, pressure, and viscosity of fluid, respectively. ( The constant-pressure boundary is adopted for the inlet boundary (Eqs. 5, 6) and exit boundary (Eqs. 6, 7) in COM-SOL software, which are as follows: The boundary conditions of the remaining boundaries are shown below:

Basic settings of flow simulation
The parameters of the fluid (i.e., density and viscosity) were determined according to the CH 4 properties under standard conditions during the simulation processes, but the influence of temperature and pressure on viscosity and density was not taken into account in this study. The model boundary conditions are loaded as shown in Fig. 9. There are three directions in the model, and two relative boundaries are selected as the

Results and discussion
According to the pore diameter, the pores can be divided into micro-pores (< 10 nm), transition pores (10-10 2 nm), mesopores (10 2 -10 3 nm), and macro-pores (> 10 3 nm; Hodot 1966;Zhao et al. 2019). The pores extracted from the BF sample in this study are macro-pores based on the pore structure classification system mentioned above. In order to better understand the geometrical and topological structures of the pore network model, the following parameters are defined (Lindquist et al. 2000;Sok et al. 2002;Herring et al. 2013;Zhao et al. 2018;Liu et al. 2019): (1) pore radius is the radius of the equivalent sphere with the same volume; (2) throat radius and length are the radius and length of the channel connecting two interconnected pores; (3) coordination number is the connective number between one pore and other pores; (4) Euler's characteristic is also known as the connectivity function, which can be obtained by calculating the relative Euler characteristics of the minimum radius of the pores and throats.

Geometric structure
The number of macro-pores gradually decreases with the increase of pore radius (Fig. 10a). When the pore radius is greater than 7 μm, the pore number is less than 10 in each range of radius (Fig. 10a). When the pore radius is approximately 2 μm, the pore volume reaches its maximum value (Fig. 10b). The equivalent radius of throats is between 1.03 and 8.9 μm, and their number decreases with the increase of radius. The number of throats with radius less than 3 μm occupies 96.07% of the total throats (Fig. 10c). The length of throats is between 3.28 and 231.25 μm, and their number also decreases with the increase of radius (Fig. 10d).

Topological structure
The coordination number plays an important role in controlling the CH 4 flow (Rabbani et al. 2014;Song et al. 2017). The greater the value is, the better the connectivity of pores will be (Lindquist et al. 2000;Sok et al. 2002;Liu et al. 2019). The coordination number of pores concentrates around three in BF sample (Fig. 11a), which indicates that most of the pores can be interconnected with other three pores, and the flow  Fig. 9 Diagram of loading condition in boundaries of three directions. Notes: X and Y direction means the direction horizontal to the bedding plane, and Z direction means the direction perpendicular to the bedding plane 1 3 path is wide. The closer the intersection point between the connectivity function and the X-axis is to zero, the poorer the pore connectivity will be (Silin and Patzek 2006;Vogel and Roth 2001). The intersection point is 3-4 μm in BF sample (Fig. 11b), which indicates that macro-pores with the width of 3-4 μm play a major role in the connectivity of pores.

Single-channel flow characteristic
The distribution of the pressure field and velocity field on single-channel flow condition along the direction Euler's characteristic, mm -3

Fig. 11
Topological structure of pore network model. a Coordination number; b Euler's characteristic number perpendicular and horizontal to the bedding plane can be obtained by analyzing the post-processing module in COM-SOL software (Figs. 12, 13, 14, 15). During the single-channel flow, the pressure-dropping (ΔP) along the direction perpendicular and horizontal to the bedding plane are the same (0.01 Pa), but the 3D distribution of the pressure in different directions is quite different, and the mean flow pressure along the direction horizontal (X-axis and Y-axis) and perpendicular to the bedding plane is 5.64 × 10 −3 Pa, 7.84 × 10 −3 Pa, and 3.08 × 10 −3 Pa,  (Fig. 12). This lies in the difference of radius, shape, and connectivity of pores and throats. Specifically, the pressure decreases gradually along the direction of the CH 4 flow, and the maximum pressure appears in the narrow pores near the inlet boundary. The pressure rapidly changes at the narrow pores, which decreases first and then increases (Fig. 12). The velocity distribution of 3D view and different sections (i.e., x = 80, 87.5, 95, 102.5 and 110 μm; y = 110, 102.5, 95, 87.5 and 80 μm; z = 110, 102.5, 95, 87.5 and 80 μm) along the single-channel flow in the direction perpendicular and horizontal to the bedding plane can be seen in Figs. 13, 14 and 15.
Although the pressure dropping in three directions are the same, the velocity distribution in the direction perpendicular and horizontal to the bedding plane are quite different (Figs. 13a, 14, 15a). Some findings are as follows: (1) Although all the pores are interconnected in the largest poreconnected group (Fig. 8), the CH 4 flow can only be observed in some pores (Figs. 13, 14, 15). The CH 4 firstly flows along the path with large radius, then the short distance from the outlet.
(2) As for the CH 4 flow in a single pore, the velocity gradually decreases from the center of the pore to the wall, and until it becomes zero (Figs. 13,14,15). In terms of the whole pore structure, the maximum velocity appears in the central region where the CH 4 flow is fully developed.
(3) Although the mean flow velocity on single-channel along the direction horizontal (X-axis and Y-axis) and perpendicular to the bedding plane is about 2.97 × 10 −12 m/s, 5.23 × 10 −12 m/s, and 1.97 × 10 −12 m/s, respectively, the flow velocity in different sections are quite different (Figs. 13,14,15). (4) In different sections, the maximum velocity appears in the center of the pore channels. Due to the rapid bending of pore channels and the rapid decrease of pore diameter, the velocity of the CH 4 flow in the narrow pores is the fastest (Figs. 13, 14, 15).
In order to further analyze the distribution of pressure field and velocity field at different slices (i.e., the distance of slices from the flow inlet slice is 0, 7.5, 15, 22.5 and 30 μm, respectively) along a single-channel in the direction horizontal and perpendicular to the bedding plane, the average value of pressure and velocity on different slices is analyzed (Fig. 16). In different directions, the distribution of mean flow pressure from the inlet to the outlet has the same distribution law, and the mean flow pressure gradually decreases with the distance from the flow inlet  (Fig. 16a). The changes of mean flow pressure in the direction horizontal to the bedding plane are relatively stable, while the changes in the direction perpendicular to the bedding plane are relatively large (Fig. 16a). In the direction horizontal to the bedding plane, the changes of mean flow velocity of various slices are relatively stable with the value of (0.4-0.6) × 10 −11 m/s (Fig. 16b). In the direction perpendicular to the bedding plane, there is a significant difference in mean flow velocity of various slices with the value of (0.05-1.2) × 10 −11 m/s (Fig. 16b). The mean flow velocity in the direction horizontal (Y-axis) to the bedding plane was the largest, followed by that in the direction horizontal (X-axis) to the bedding plane, and the mean flow velocity in the direction perpendicular to the bedding plane was the smallest (Fig. 16b), which is consistent with the distribution of mean flow pressure in different directions (Fig. 16).

Dual-and multi-channel flow characteristics
For the dual-channel and multi-channel flow, the former means the CH 4 flows along the X + Y, X + Z and Y + Z directions at the same time, and the latter means the CH 4 flows along the X + Y + Z directions simultaneously, so the direction of the dual-channel flow is the X + Y, X + Z and Y + Z directions, and the direction of the multi-channel flow is the X + Y + Z directions (Figs. 17, 18, 19). The distribution of pressure field and velocity field during the dual-channel and multi-channel flow can be seen in Figs. 17, 18 and 19.
The mean flow pressure along the X + Y, X + Z, Y + Z directions in the dual-channel flow, and along the X + Y + Z directions in the multi-channel flow is 4.83 × 10 −3 Pa, 6.92 × 10 −3 Pa, 7.16 × 10 −3 Pa, and 5.71 × 10 −3 Pa, respectively. Similar to the pressure condition in the single-channel flow along the X, Y, and Z directions, the pressure decreases  Fig. 18 Distributions of velocity field along dual-channels X + Y, X + Z, and Y + Z directions. Notes: X and Y direction means the direction horizontal to the bedding plane, and Z direction means the direction perpendicular to the bedding plane gradually along the CH 4 flow direction in the dual-channel and the multi-channel flow, and the maximum pressure occurs at the inlet boundary. Different from the pressure condition in the single-channel flow along the X, Y, and Z directions, the distribution area of larger pressure in the dual-channel and the multi-channel flow is relatively wide (Figs. 17, 19a).
The mean flow velocity along the X + Y, X + Z, Y + Z directions in the dual-channel flow, and along the X + Y + Z directions in the multi-channel flow is 1.55 × 10 −11 m/s, 9.12 × 10 −11 m/s, 1.94 × 10 −11 m/s, and 10.75 × 10 −11 m/s, respectively. Compared with the velocity condition in the single-channel flow along the X, Y, and Z directions, the velocity increases obviously under the dual-channel and the multi-channel flow (Figs. 18,19b), and the maximum velocity is nearly 4.31 × 10 −9 m/s in the X + Z and X + Y + Z directions, which indicates that the development of interconnected pores with various directions is conducive to CH 4 migration and production.

Analysis of uncertainty
The 3D reconstruction and visualization results of the macro-pores are more or less affected by the man-made differences in image processing, and the small sample size simulated will affect the simulation result of CH 4 flow, but the real intention of this study is to analyze the 3D spatial structure of the macro-pores and CH 4 flow simulation in anthracite coal based on μ-CT scanning data.
In order to reduce the negative effect caused by the uncertainty, multiple methods should be used in conjunction during the 3D visualization and reconstruction processes of digital core. So, the optimal 3D reconstruction and visualization results can be obtained after image processing. Further improvement of the hardware equipment of the data processor can further expand the sample size simulated, so that the simulation effect of CH 4 flow will be better.

Conclusions
In this study, the BF coal sample was firstly scanned by the X-ray μ-CT technology to obtain the 2D CT slices. Secondly, the original 2D CT slices were reconstructed and divided to obtain the pore network model. Then, a mathematical model of CH 4 flow simulation with macro-pores as the carrier was established. Finally, the CH 4 flow simulation and its visualization were realized based on the perfect docking of the X-ray μ-CT technology and the COMSOL. The main conclusions are as follows: 1. The actual reconstructed volume of digital core is 250 × 250 × 250 μm 3 , and the number of pores and throats gradually decreases with the increase of radius. The equivalent radius of pores and throats are 1-16 μm and 1.03-8.9 μm, respectively, and the throat length is 3.28-231.25 μm. The coordination number of pores concentrates around three, and the intersection point Distributions of pressure field and velocity field in X + Y + Z directions. a Pressure filed, b velocity filed. Notes: X and Y direction means the direction horizontal to the bedding plane, and Z direction means the direction perpendicular to the bedding plane between the connectivity function and the X-axis is 3-4 μm, which indicate the macro-pores have good connectivity in 3D space. 2. Regarding the single-channel flow, the pressure decreases gradually along the direction of CH 4 flow. The CH 4 firstly flows along the path with large radius and the short distance from the outlet. The flow velocity of CH 4 gradually decreases from the center of the pore to the wall. Compared with the single-channel flow, the pressure also decreases gradually along the CH 4 flow direction, while the velocity increases obviously in the dual-channel and multi-channel flow. 3. The mean flow pressure gradually decreases with the increase of the distance from the inlet slice. The change of mean flow pressure is relatively stable in the direction horizontal to the bedding plane, while it is relatively large in the direction perpendicular to the bedding plane. The mean flow velocity value in the direction horizontal to the bedding plane (Y-axis) is the largest, followed by that in the direction horizontal to the bedding plane (X-axis), and the mean flow velocity value in the direction perpendicular to the bedding plane is the smallest, which is consistent with the distribution of mean flow pressure in different directions. 4. The perfect docking of the X-ray μ-CT technology, the image processing technology and the COMSOL software can not only realize the 3D visualization reconstruction of pores and throats, but also complete the CH 4 fluid simulation with the interconnected pores as the carrier. Multiple methods should be synthetically used during the 3D visualization and reconstruction processes of digital core, and the hardware equipment of the data processor should be further improved.