Microstructure-based multi-scale evaluation of fluid flow in an anthracite coal sample with partially-percolating voxels

Understanding fluid flow behavior in coal is of great significance for coal-bed methane exploration. X-ray CT and image segmentation have been widely used to extract pore network and generate flow field grids for flow simulation in coal samples. However, these techniques have fundamental limitations for the multi-scale characterization of coal samples, where the sub-voxel scale details could not be resolved for millimeter scale macroscopic samples. This makes it difficult to simulate the multi-scale flow behavior of fluid transport in coal sample with varying pore scales. The primary challenge is to make connection between simulation results of different scales. In the present work, multi-scale fluid flow in an anthracite coal sample was simulated by incorporating the data-constrained modeling (DCM), molecular dynamics (MD) method and partially-percolating lattice Boltzmann method (PP-LBM). In this multi-scale simulation method, three-dimensional (3D) flow field containing multi-scale structural information of the coal sample was generated by combining DCM with multi-energy synchrotron radiation CT. Multi-scale fluid flow was simulated by PP-LBM. In PP-LBM, an effective percolation fraction parameter which represents the effective volume fraction of the fluid that contributed to the flow for the voxel was used as a bridge to connect the fluid flow pattern of sub-voxel scales and voxel scales. The effective percolation fraction of a voxel versus its porosity was derived by MD simulations at the sub-voxel size level. The 3D distribution of fluid speed in the coal sample and its permeability were obtained by this multi-scale method. The numerical results are consistent with published laboratory measurements. Our proposed approach incorporated multi-scale effects and offered a more realistic fluid transport simulation method for a coal sample with varying pore size scales from the microscopic to macroscopic level. The method would be applicable for fluid transport simulations for other multi-scale porous materials.


Introduction
The understanding of fluid transport properties in coal is relevant for coal-bed methane exploration and the development of clean coal technology. Fluid transport in coal is affected by many factors including lithotype [1], pore size and porosity [2,3], interaction between fluid and coal [4], shrinkage of coal during fluid extraction process [5], fines production [6], coal surface roughness properties [7], relative permeability of different fluid phases [8]. Furthermore, coal is a dual porosity system which contains micro-pores and cleats. Cleats are main seepage channel of fluid and micro-pores are important control factor for fluid flow. X-ray CT as a nondestructive characterization technique can provide three-dimensional (3D) structural information of sample. Many researchers have used this technique for extraction of pore networks in coal [4,[9][10][11][12][13][14], studying gas adsorption and desorption in coal [15][16][17], and evaluation of fluid flow in coal [1,2,5,14,18]. Although extensive research in this area has been carried out, a severe problem still exists in the application of x-ray CT for coal which is an dual porosity medium, in that all pores are difficult to be captured using a single imaging technique. As has mentioned above, the pore sizes in coal ranging from nanometers to millimeters are all relevant. However the existing x-ray CT imaging technique has a limited spatially dynamical range. For a millimeter-sized sample, the imaging pixel size cannot be smaller than micron-or sub-micron scale. A higher imaging resolution is unpractical because that the data size and experiment time will increase exponentially with the improvement of imaging resolution. With the main-stream approach, the structural information below the imaging resolution is lost in the later image threshold segmentation process. Consequently, the coal microstructure as characterized with the conventional technique is inadequate for fluid transport simulation as fine length scale contributions could not be accounted for. For a long time, the difficulty of multiple scales characterization of sample 3D internal structure has restricted the development of multi-scale flow simulation.
Recently, a method of data-constrained modeling (DCM) combined with multi-energy synchrotron x-ray CT has gained a degree of popularity in multi-scale characterization of material microstructures [19]. One of the distinctive advantage of the DCM approach as compared with the conventional image segmentation technique is that it quantifies fractional material and pore distributions finer than x-ray CT imaging resolution. This makes it possible for incorporating fine length scale effect in fluid transport simulation. DCM has been applied successfully to multi-scale microstructure characterization of a broad range of materials [20][21][22][23][24]. In our previous study [10], this approach has been used for multi-scale characterization of an anthracite coal sample, where the 3D structure of the sample was generated, including the volume fractions of those compositions for each imaging voxel. This opens a new opportunity for accurate multi-scale simulation of fluid flow in coal.
As a numerical technique for modeling fluid transport in porous media, the lattice Boltzmann method (LBM) has been accepted for simulating fluid flow in porous media [25,26]. LBM simulates fluid particles propagation on a cubic lattice which consists of cubic voxels, which is equivalent to discretized solution of the Navier Stocks equation. The states of the fluid particles are described by a distribution function. The interactions between fluid particles are incorporated through collisions at lattice nodes (voxels) such that mass and momentum are conserved. Thus, LBM requires an explicit description of the pore-space geometry which is often obtained through x-ray CT and image segmentation.
In order to improve computational efficiency for simulating fluid flow through a porous media, a partial-bounce-back LBM was proposed by Dardis and McCloskey in 1998 [27, 28], where neighboring voxels were grouped to reduce the number of computational nodes. Unlike the traditional LBM, where a lattice is comprised of nodes that are either solid or fluid, this approach uses a probabilistic or partial-bounce back model, where lattice node properties are varied to reflect the local solid fraction of the material [29][30][31].
Recently, a partially-percolating lattice Boltzmann method (PP-LBM) was formulated by Li et al [32], which combines the partial-bounce-back LBM with the fluid-solid boundaries proposed by Walsh et al [29]. In PP-LBM an effective percolation fraction (p f ) is introduced to represent the effective fluid volume fraction that contributes to flow. For a discretely segmented 3D image, p f takes a value of either 0 or 1. For a partially porous voxel as obtained by DCM for a multi-scale materials such as coal, p f can take any value between 0 and 1. PP-LBM has been used for flow simulations for various tight oil and gas reservoir rocks [32].
In this paper, PP-LBM was used to simulate single phase fluid (water) flow in a realworld anthracite coal sample collected from Qinshui basin of China. In contrast to previous work [32], this study took the p f as a function of pore fractions in each voxels, which enables a more accurate permeability simulation. Then molecular dynamics (MD) simulation was employed to obtain the quantitative relationship between p f and porosity. This approach was able to incorporate rich sub-voxel information into the flow simulation. Therefore, it more accurately represents fluid flow in coal.

The coal sample
The anthracite coal sample used in this study is the same one as used by Wang et al [10]. It was collected from Yangquan mine, which belongs to Qinshui basin. It is one of the gas rich coal basins in China. The vitrinite reflectance and carbon content of the sample are 2.3% and 89.0%, respectively. The total porosity of coal sample is 4.43%, which was obtained according to the Chinese standard GB/T 23561. . The x-ray CT data of the sample was collected at the BL13W beam-line in the Shanghai Synchrotron Radiation Facility. The sample was imaged under the monochromatic x-ray beam energy of 14 keV, 18 keV, and 30 keV respectively with the imaging pixel size of 3.7 μm. The DCM method was applied to obtain the 3D multi-scale microstructure of the sample. DCM is based on the assumption that the internal distribution of compositions of a material sample follows the general principle of statistical physics, and such distribution is consistent with multi-energy x-ray CT imaging measurements. Numerically, voxel compositions can be solved approximately with a nonlinear optimal programming method using the DCM software. As the DCM calculation incorporate partial volume effect in the CT imaging process, the compositions information smaller than the imaging resolution size is preserved. In this study, a rectangular region of 100 × 100 × 20 voxels was selected for the following PP-LBM calculation. The physical dimension of each voxel is 3.7×3.7×3.7 μm 3 . Figure 1 is the spatial distribution of voxel porosity and mineral volume fractions in the coal sample. The DCM computed porosity of this region is 6.47%. Detailed DCM calculation process and parameter sets can be found in [10]. The x-ray CT data and more DCM computed results are available from CSIRO Data-Access-Portal [33]. It can be seen that in the sample, the pore and heavy minerals are mainly distributed in the dispersion pattern, while the light minerals is relatively concentrated. The display intensities of each composition in different areas in figure 1 are different and vary continuously, especially for the pores. This indicates that the most voxels are partially porous at imaging resolution of 3.7 μm. This is an indication that the pores and other compositions cannot be distinguished using image threshold segmentation.

The PP-LBM method
Fluid flow between voxels is simulated by PP-LBM. A brief account of background information for this study to be based on is given below. A more detailed description of this method can be found in [32]. In the PP-LBM, the states of the fluid particles in each voxel are represented by the distribution function a ( ) f x t , , where x is the spatial coordinate of a voxel on the lattice, t is the time, a represents 19 directions of D3Q19 velocity model. In the flow simulation, fluid particles migrate between neighboring voxels. The movement of each fluid particle at every time step is divided into three steps: streaming step, collision step, and porous media step. At the streaming step, fluid particles migrate between neighboring voxels. At the collision step, the fluid particles collide on individual voxel. At the porous media step, the fluid particles rebound partially. The distribution functions after these three steps are given by equations (1)-(3): , , 1 In equation (3), p f is the effective percolation fraction of the voxel at x as mentioned above [32]. In order to incorporate the influence of the fine structure smaller than the CT imaging resolution on fluid flow, p f was taken as a function of pore volume fraction in this paper. The value of p f will be derived by MD simulations in the next section.

Molecular simulation method
In the molecular simulations, the modeled coal molecules are flexible and movable. They are described using the GROMOS force field [34]. In the coal model, the bonded interactions include bond, bond angle, harmonic (improper) dihedral angle, and trigonometric (torsional) dihedral angle terms. The non-bonded interactions are the sum of van der Waals and electrostatic (Coulomb) interactions between all pairs of atoms. Detailed parameter sets can be found in [34]. In contrast to the parameterization of other biomolecular force fields, this parameterization of the GROMOS force field is based primarily on reproducing the free enthalpies of hydration and solvation for a range of compounds. The relative free enthalpy is a key property in many biomolecular processes of interest and is why this force field was selected [35].
Methane molecules are used to investigate the effect of porosity on fluid diffusivity and based on their relationship the effective percolation fraction as a function of porosity was extracted. In the MD simulation, methane molecules are treated as united atoms. The nonbonded interactions between atoms which are separated by more than three bonds, or belong to different molecules, are described by pair wise-additive Lennard-Jones (LJ) 12-6 potentials. Cross-interactions are calculated by the Lorentz-Berthelot mixing rules. The LJ parameters and energy terms and the parameters for coal are taken from [34]. The Coulomb interactions in the system are calculated by Ewald summation for periodic systems [36]. A truncated and shifted potential was used with a cutoff radius of 1.4 nm to reduce the computational effort. Interactions for separation of longer than 1.4 nm were omitted from the energy and force computations as beyond this distance the potential and force become negligible.
Simulation procedure consists of Monte Carlo simulations performed in the grand canonical ensemble (GCMC) in which the chemical potentials of the CH 4 , the total volume, and the temperature of the system were fixed and MD simulation in the constant volume and temperature ensemble. The former used the open source package RASPA 1.0 developed by Dubbeldam et al [37]. The latter was carried out using Gromacs software [38][39][40]. The combined GCMC and MD methods were applied to the study of methane adsorption on dry and moist coal [35]. In our simulation, the Berendsen thermostat [41] was applied to keep the temperature constant (310 K). The GCMC allowed us to calculate the amount of methane under different coal volume fraction. In other words, The number of methane molecules adsorbed in coal was computed under varying number of coal molecules at the reservoir pressure of 8.92 MPa and temperature of 310 K. This temperature and pressure conditions correspond to the geological depth of 900 m [42]. GCMC simulation was switched to MD simulation to calculate the diffusion coefficient of methane under different methane volume fraction. Details of the combined GCMC and MD method is available in [35].
The initial configuration consisted of disordered coal molecules. All molecules were randomly placed in an empty space of the periodic simulation box which has x, y, zdimensions of 2.88×3.10×2.88 nm 3 . Periodic boundary conditions were applied in three directions to reduce size effect. The equations of motion were integrated with a time step of 0.001 ps. The simulation included equilibrium run (10 ns) and production run (50 ns). To get good statistics, simulations should run sufficiently long enough. We run MPI parallel programming with four processors in the Raijin Cluster in the National Computational Infrastructure Australia. Figure 2 shows the relationship between fluid volume fraction, defined as the ratio of fluid volume to total volume, and fluid diffusion coefficient obtained by MD. Molecular simulation results indicates that at low fluid volume fraction (or low void fraction), fluid diffusion coefficient is small as fluid molecules are at adsorbed state. As fluid volume fraction increases (over than 45.2%), the fluid volume fraction (or void fraction) effect on the diffusion coefficient becomes significant. The diffusion coefficient increases with the void fraction.

Effective percolation fraction p f
The parameter p f for each voxel in equation (3) represents the effective volume fraction of the fluid which contributes to the flow for the voxel. It is generally lower than or equal to the voxel porosity as some pore may not fully contribute to the flow. It is an intrinsic parameter in (a) In the MD simulation, the fluid volume fraction is equivalent to pore volume fraction in coal as the pores in the coal sample are the space available for the fluid. (b) The relationship between diffusion coefficient and porosity derived by molecular simulation is applicable to the CT imaging voxel size as the diffusion coefficient obtained from MD simulation can represent the diffusion in meso-and macroscopic scales. (c) The pores in coal which contribute to diffusion also contribute to flow. The basis for this assumption is as follows: In the PP-LBM simulation, the fluid obeys Darcy's law, which can be expressed as [43]:  (4) and (5) have the same mathematical formulation, we assume that the connected pores in coal equally contribute to the diffusion to the flow. Although CH 4 diffusion and water flow are different physical phenomena, the effective percolation fraction p f would be equally applicable to both of them, as diffusion and flow use the same connected pores. After measured 47 different rock samples, Reimus et al [45] also found that the log of normalized diffusion coefficient (tracer matrix diffusion coefficients divided by the free-water diffusion coefficient) is positively correlated with the log of matrix permeability. Their results indicate that matrix diffusion coefficients, like matrix permeability, have a greater dependence on the interconnectedness of matrix porosity than on the matrix porosity itself.
Considering that the value of p f cannot be higher than 1, the fluid diffusion coefficient in figure 2 was normalized by its maximum value. Based on the above assumptions, the relationship between p f and the voxel porosity was obtained as shown in figure 3.
Using the DCM derived voxel porosity values v n where n is the positional coordinate of the voxel, the value of ( ) p f n for the voxel is estimated as:

PP-LBM simulation
After the values of ( ) p f n were obtained, the PP-LBM was used to simulate water flow through the coal sample. The speed distribution of the coal sample is shown in figure 4. In the simulation, a constant pressure difference 1.5×10 −5 Pa was applied on the two opposite external surfaces in z-direction (figure 4). Other parameters of the PP-LBM in the simulation process are the same as those in previous publication [32]. The simulation was completed after reaching the following convergence criterion: where k(t) represents the calculated bulk permeability at time t. The simulation converged after 3 607 000 iterations. The bulk permeability (k) of the sample was calculated based on the generated speed datasets using Darcy's law (equation (22) in our previous publication [32]). Figure 4 shows the computed 3D distribution of fluid speed. The average speed is 8.15×10 −12 m s −1 and the bulk permeability is calculated to be 5.30×10 −6 Darcy. They are at the same magnitude as the results reported by others [14,46,47]. It can be seen that most region in figure 4 are of high grayscale value, only a few areas are of low grayscale  value, indicating that the flow field is largely inhomogeneous. The bulk permeability is dominated by small portion of the sample. It should be noted that, as mentioned in [10], most of the pores in the coal sample used in this paper exist in the form of partial volume forms, and if the traditional fluid flow simulation method is used, the permeability would be a very low value. Figure 5 shows the computed two-dimensional (2D) distribution of porosity and fluid speed in the 17th slice. Comparing figures 5(a) and (b), it can be found that only a small portion of porosity contributed to the flow. Some areas with high porosity cannot be reached by fluid, while some areas with low porosity have a high fluid speed distribution. This indicate that the 2D distribution of fluid speed is not directly proportional to the 2D distribution of porosity. For this low permeability and low porosity coal samples, fluid permeability is more dependent on pore connectivity. Figure 6 illustrates the 3D spatial correlation between porosity and flow. The high porosity voxels are displayed as blue and the high fluid speed voxels are highlighted as red. It can be seen that not all the pores contribute to the flow equally. For example, although there is higher porosity in area marked as 'a' (surrounded by white line in figure 6), almost no flow through the region. This is an indication that the 3D spatial distribution of the pores would be an important factor in determining the bulk permeability. Figure 7 shows the 3D distribution of fluid speed (red area) and minerals (blue and green area) in the coal sample. It can be seen in figure 6 that most flow distributes at those pores  which locate at the edge of minerals. The reason for this is that higher porosity and better connectivity are found at the edge of minerals than other area of the sample [10].

Conclusions
Combined GCMC and MD simulations were carried out to determine the relationship between the fluid diffusion coefficient of a voxel and its porosity. Utilizing this relationship and the coal sample 3D structure characterized by DCM, a 3D distribution of effective percolation fraction has been estimated. PP-LBM was used to simulate the water flow in the  In the image, the fluid speed voxels are displayed as red, the light minerals (Illite, Quartz and Kaolinite) are displayed as blue, and the heavy minerals (Chlorite and Titania) are displayed as green. In a voxel the displayed intensity of red is proportional to the speed value. Volume fractions of the material phases in a voxel are proportional to their respective display intensities. Co-existence of multiple items is shown as a mixed color.
sample. The incorporated GCMC, MD, DCM, and PP-LBM method allows a multi-scale flow simulation for the sample. The 3D distribution of fluid speed has been obtained and our results have demonstrated that: (1) The sample has a low water permeability of 5.3×10 −6 Darcy; (2) Only a small portion of porosity contributed to the flow, and the flow field distribution is inhomogeneous; (3) For this low permeability and low porosity coal samples, fluid permeability is more dependent on the 3D connectivity than on the porosity; (4) The flow speed is high at the boundary of mineral-rich regions.
It should be noted that the effective percolation fraction ( ) p f n is a key parameter for multiscale flow simulation. As mentioned above, fluid flow in coal is affected by a variety of factors, and accurate simulation of fluid flow behavior faces multiple challenges. Among these challenges, multi-scale flow simulation is a fundamental problem. In this paper, considering that the flow of methane and water in coal has been widely concerned by many scholars, therefore methane is used as a fluid to obtain the relationship between void fraction and effective percolation fraction. For the multi-scale flow of water, the simulation results would be more accurate if the effective percolation fraction at each voxel was obtained using water as a fluid in MD simulation. More accurate simulations should also consider other factors such as nonlinearity in fluid flow, fluid-solid interactions, and composition fraction. If more factors affecting flow can be included in the effective percolation fraction, the accuracy of multi-scale simulation would be effectively improved.