Pore-Scale Lattice Boltzmann Simulation of Gas Di ﬀ usion–Adsorption Kinetics Considering Adsorption-Induced Di ﬀ usivity Change

: The di ﬀ usion–adsorption behavior of methane in coal is an important factor that both a ﬀ ecting the decay rate of gas production and the total gas production capacity. In this paper, we established a pore-scale Lattice Boltzmann (LB) model coupled with ﬂuid ﬂow, gas di ﬀ usion, and gas adsorption–desorption in the bi-dispersed porous media of coalbed methane. The Knudsen di ﬀ usion and dynamic adsorption–desorption of gas in clusters of coal particles were considered. use of a static isotherm adsorption equation may be incorrect.


Introduction
In recent years, the adsorption-desorption behavior of gases in porous media has gained wide attention, and many theories such as single-layer adsorption, multilayer adsorption, capillary condensation, and pore-filling have been successfully applied to many areas [1,2]. The adsorption effect of gas-solid in nanopores of the unconventional natural gas (UNG) reservoirs is an important factor that both affecting the decay rate of gas production and the total gas production capacity [3,4]. However, for the UNG, including shale gas, coalbed methane, due to the involved complex pore geometries and multi-scale spatial distribution, which makes it difficult to reveal the underlying characteristics of the gas adsorption and storage in the reservoir. Therefore, further research and deeper understanding in this area can better cope with the challenges in the development and production of UNG.
The coal matrix provides abundant adsorptive sites for gas molecules, which contributes greatly to methane storage. Methane in the coal reservoir is mainly present in three states: (1) dissolved in water, (2) as the free gas in the fractures or pores, and (3) adsorbed on the inner surface of micro-fractures or nanopores [2]. The amount of gas in the adsorption state accounts for more than 85% of the total gas [5]. In a typical UNG production process, free gases in fractures or large pores first diffuse out, which was related to the rapid production cycle of a gas productivity curve [6]. However, due to the limit of the gas diffusion rate and the desorption kinetics of gas molecules on the micropore surface [7], the desorbed gas will dominate the subsequent gas production process. Accordingly, understanding the gas adsorption/desorption behavior in the matrix cannot only obtain the main controlling factors on gas adsorption/desorption and provide a theoretical basis for developing new methods to enhance gas desorption but also can supply a reference for the optimization of production process design and evaluation of the total gas content in reservoirs.
The migration and storage of methane in the coal matrix are closely related to pore size [8]. The International Federation of Applied Chemistry (IUPAC) classified the pore into three types: micropore (less than 2 nm), mesopore (2-50 nm), and macropore (more than 50 nm) [9]. Typically, Coal is a typical tight porous media, with a wide distribution of pore size, which is distributed in the range of 1-100 nm and contains a great number of micropores, some mesopores, and a small number of macropores [10]. Meng [5] analyzed the pore distribution of coal samples (from the Xishan coal mining area, China) by using the static nitrogen adsorption capacity method (SY/ T6154-1995) and found that the average radius of pores in coal ranged from 7.729 to 15.338 nm, with the largest proportion of micropore reaching 65.4%. Since the wide distribution of pore size, it is difficult to characterize the pore structure of the coal matrix, which makes it complex to understand the behavior of gas flow in the reservoir. At present, the most common method for simplifying the coal matrix is to treat the matrix as a bi-dispersed porous media, which with the heterogeneous geometry and the homogeneous material, and regardless of the specific structure and distribution characteristics of micropores inside it, only considering its statistical characteristics, namely the porosity and the average pore diameter [11][12][13][14][15][16]. Zhao [17] tested six groups of coal samples with middle and high rank by applying synchrotron small-angle X-ray scattering (SAXS) and found the average pore sizes were 2.9, 5.1, 6.4, 6.9, 14.3, and 23.1 nm, respectively, and obtained the logarithmic normal distribution of pore sizes based on Beaucage [18].
Gas adsorption behavior, which involved some complicated adsorption mechanism, is a key factor in the coalbed methane exploitation process. However, most of the previous studies [19,20] for investigating the gas adsorption behavior in coal were mainly based on the isothermal adsorption equation, which can give an important reference for reservoir evaluation and numerical modeling. Nevertheless, the isothermal data only provide static and macroscopic information, ignoring the dynamic processes of adsorption-desorption and the detailed information of the nanoporous surface. Do and Wang [21] observed the adsorption delay of activated carbon in the process of measuring the adsorption and desorption equilibrium data. They believed that the adsorption of the pore surface was non-uniform, which resulted in gas adsorption occurs not instantaneously and took some time to Energies 2020, 13, 4927 3 of 18 reach equilibrium. These above indicated that the time scale between the process of gas adsorption and gas diffusion was comparable [21,22]. Besides, gas flow in coal reservoirs usually involves multiple physical processes, including fluid dynamics, thermodynamics, chemical kinetics, and electrodynamics (because most natural media surfaces are charged), all of which are ultimately governed by interface phenomena at the pore-scale [23]. Since the length scale of the coal matrix is much larger than the typical pore or mineral particle sizes, the pore size inside the matrix differs greatly from the outside. Therefore, the continuity formula based on the spatial average is often employed to characterize the matrix, and the spatial heterogeneity of the smaller scale can be ignored.
The present work is based on the Lattice Boltzmann method (LBM). LBM originated from classical statistical physics and is a mesoscopic method based on simplified dynamic equations. In LBM, fluid flow is simulated by a series of virtual particles that propagate and collide on discrete lattice domains. Through rigorous mathematical analysis, mesoscopic continuity and momentum equations can be obtained from this propagation and collision dynamics. Particle properties and local dynamics provide advantages for complex boundaries and parallel computing. Besides, the dynamic nature of LBM makes new physical considerations in the LBM framework easy, which is especially useful for building multi-physical phenomena. Using the multi-scale expansion technique (Chapman-Enskog), the LB equation can be recovered to the Navier-Stokes equation under the incompressible limit, and the LB equation of mass transport can be recovered to the advection-diffusion equation [24,25]. Nowadays, the LB method has been successfully applied to a variety of fluid flow and gas transport phenomena, such as fluid flow, turbulence, multiphase multicomponent flow, particle suspension, heat transfer, and diffusion-reaction flow in porous media [26][27][28][29][30][31].
The adsorption characteristics of coalbed methane in the reservoir are affected by many factors, such as fluid flow, gas transport, gas adsorption/desorption, etc. In this paper, a multi-component gas flow-diffusion-adsorption coupled LB model was established to investigate the effects of fluid flow, gas transport, adsorption/desorption on adsorption, and adsorption-induced matrix diffusion coefficients and porosity.

Lattice Boltzmann Method
Coal is a low-permeability porous media with pore systems ranging from fractures (cleats) to nanopores. The present work focuses on the modeling of multi-scale systems with considering the gas viscous flow in macropores and Fick diffusion in the matrix, besides, the Knudsen diffusion and adsorption kinetics in nanopores inside the matrix is also included in the simulation. Due to the gas velocity under typical reservoir condition is very low and the Mach number is far less than 0.3, it acceptable to describe the methane flow in the coal matrix under the incompressible limit [24].

The LB Equation for Fluid Flow
The LB method for gas viscous flow at pore-scale can be expressed as follows [32]: where f i is the discrete density distribution function, f eq i is the local equilibrium function, e i is the discrete velocity of the particle, c s is the lattice sound velocity, c is the sound velocity, τ is the relaxation time, ω i is the weight coefficient, and δ t is the time step.

The LB Equation for Gas Diffusion-Reaction
For the process of methane transport, it is acceptable to employ the passive scalar LB model due to the gas velocity in the reservoir is low enough and its effect on solution density and velocity can be negligible [25]. Therefore, gas transport in the reservoir with adsorption can be described as follows: where g i is a discrete concentration distribution function; g eq i is the corresponding local equilibrium distribution function; τ g is the diffusion-related relaxation time; R s is the source/sink term associated with gas adsorption/desorption.  (3) and (4) can recover the advection-diffusion equation with a source/sink term.
The diffusion coefficient is a key parameter for the gas transport process in porous media. According to the previous study [33], diffusion coefficients can be divided into Knudsen diffusion coefficients (KDC), corrected diffusion coefficients (CDC), and Fickian diffusion coefficients (FDC). The KDC can be employed to describe the phenomenon of random walk and the CDC corresponding to the chemical potential-driven diffusion, and the FDC is usually adopted to concentration-driven diffusion. The KDC is more closely related to the microscopic features of molecules, while FDC is more directly related to the macroscopic transport of mass, which is more experimentally available.
At pore-scale, gas diffusion coefficients related to methane migration in the coal matrix includes FDC and KDC. The FDC can be obtained by experiment, and based on the parallel capillary model of porous media, the KDC can be written as [34]: where R is the gas constant, T is the temperature, M is the molar mass of the gas, and r is the radius of the capillary. When adsorption is considered, since gas molecules are adsorbed on the inner surface of the nanopore, the cross-sectional area of the free gas transfer is reduced, and the porosity associated with the pore pressure in the coal matrix is decreased. Xiong et al. [35] proposed an effective capillary radius of the adsorption layer for the single capillary of porous media: where d m is the diameter of the methane molecule, e.g., d m = 0.38 nm [35]; p is the pressure, Pa; p L is the Langmuir pressure, Pa; and p/(p + p L ) is the adsorption saturation of the pore surface. The model is based on the Langmuir isotherm adsorption equation, which corrects the influence of gas molecular radius and adsorption saturation on the effective radius of the pore but ignores the time effect of concentration/adsorption in the adsorption process. In the present work, the adsorption Energies 2020, 13, 4927 5 of 18 saturation is introduced by the ratio of the instantaneous adsorption amount to the saturated adsorption amount, and the effective pore radius is then considered: where V m and V are the saturated adsorption amount and the adsorption amount at a current time, respectively. The modified effective porosity of the capillary can be written as: Due to the complex geometric composition of coal, the effects of porosity and tortuosity of the porous media also should be considered. Therefore, effective KDC can be given by: 3 where ε 0 is the initial porosity of the coal matrix and τ w is the tortuosity.

Langmuir Adsorption Kinetic Equation
The adsorption mechanism of gas in coal is very complicated, and various adsorption forms coexist [2], such as monolayer adsorption, multi-layer adsorption, capillary condensation, and capillary filling, etc. At present, many adsorption models have been established and applied to the methane-coal adsorption, among them, the Langmuir isothermal adsorption model is the most commonly used one [2,36,37]. In this paper, we employ the Langmuir adsorption kinetics equation to control the evolution of gas adsorption-desorption: where k a and k d are the adsorption and desorption rate constants, respectively. The methane-coal adsorption of coal is typical physical adsorption and is a reversible process [21], which means that there are simultaneous gas adsorption and desorption in adsorption sites. At a specific time, if ∂V/∂t < 0, the adsorption rate is lower than the desorption rate, the gas desorption from pores; if ∂V/∂t > 0, the adsorption rate is greater than the desorption rate, the net gas adsorption amount increases; when ∂V/∂t = 0, the right part of the Equation (10) can recover the classical Langmuir isotherm adsorption equation: where p is the gas pressure, p = CMc 2 s , Pa; C L is the concentration corresponding to the Langmuir pressure, It can be found that p L is only related to the ratio of k a and k d and independent of their magnitude. The value of k a , k d , p L can be determined based on the experiment, and anyone can be obtained from the other two of them.
To consider the gas-solid dynamic adsorption process, we developed a LB model to realize the adaptive conversion of gas in porous media between adsorption and desorption based on the model proposed by He [38]. In the model, each site in the porous media can be considered as a good adsorption position, the Langmuir rate equation can be incorporated into Equation (3) through a source/sink term: In the evolution of the LB equation, the amount of gas adsorption/desorption is updated with time, and the new amount can be obtained from a first-order difference scheme of Equation (10): Recently, the bi-dispersed porous media has been drawn wide attention due to the similar pore size distribution and fractal characteristics with the geometric characteristics of coal. The bi-dispersed porous structure of the coal matrix, as shown in Figure 1, is composed of clusters, which are agglomerated by small particles. In the bi-dispersed porous media, there are numerous intraparticle pores within the clusters, which contain micro-and mesopores, and some interparticle pores between the clusters, mainly macropores. The previous study has been reported that the micro-and mesopores will play an important role in methane adsorption, and the macropores act primarily as transport channels [2]. The fluid flow and gas diffusion in macropores are governed by the double distribution LB equation; due to the pore size of the clusters being extremely small, the effect of Knudsen diffusion and adsorption should be considered.
Energies 2020, 13, x FOR PEER REVIEW 6 of 18 Recently, the bi-dispersed porous media has been drawn wide attention due to the similar pore size distribution and fractal characteristics with the geometric characteristics of coal. The bi-dispersed porous structure of the coal matrix, as shown in Figure 1, is composed of clusters, which are agglomerated by small particles. In the bi-dispersed porous media, there are numerous intraparticle pores within the clusters, which contain micro-and mesopores, and some interparticle pores between the clusters, mainly macropores. The previous study has been reported that the micro-and mesopores will play an important role in methane adsorption, and the macropores act primarily as transport channels [2]. The fluid flow and gas diffusion in macropores are governed by the double distribution LB equation; due to the pore size of the clusters being extremely small, the effect of Knudsen diffusion and adsorption should be considered.

Physical Model and Verification
In this paper, the D2Q9-LB model is used to describe fluid flow and mass transfer. Firstly, the validity of the model is verified by two classic single-channel steady-state flow models, the diagram of the convection-diffusion system with surface adsorption reaction is shown in Figure 2. The size of the first model is set to 200 × 100, the parabolic flow rate is set at the inlet, umax = 0.06; the initial concentration in the field is 0, the gas with a diffusion coefficient of 1/6 spreads from the left, the inlet concentration is 1, the outlet boundary is set to / 0 n C ∂ ∂ = , and the Henry adsorption kinetic boundary condition is set at the bottom boundary [40,41]: The Lévesque analytical solution after the adsorption has stabilized is: The comparison between the LB results and the analytical solution are shown in Figure 3a. The

Physical Model and Verification
In this paper, the D2Q9-LB model is used to describe fluid flow and mass transfer. Firstly, the validity of the model is verified by two classic single-channel steady-state flow models, the diagram of the convection-diffusion system with surface adsorption reaction is shown in Figure 2. The size of the first model is set to 200 × 100, the parabolic flow rate is set at the inlet, u max = 0.06; the initial concentration in the field is 0, the gas with a diffusion coefficient of 1/6 spreads from the left, the inlet concentration is 1, the outlet boundary is set to ∂C/∂n = 0, and the Henry adsorption kinetic boundary condition is set at the bottom boundary [40,41]: Energies 2020, 13, 4927 The Lévesque analytical solution after the adsorption has stabilized is: The comparison between the LB results and the analytical solution are shown in Figure 3a. The LB result agrees well with the analytical solution except for slight deviations at the entrance due to the singularity.
Energies 2020, 13, x FOR PEER REVIEW 7 of 18 The size of the second model is set to 100 × 100. There is a 40 × 40 solid block in the center of the flow field, which represents a cluster of coal particles. Assuming that the solid block is a homogeneous material, the internal micro-and nanopores adsorb gas molecules follow the Langmuir adsorption rate equation [42]. The fluid is driven by a pressure gradient, 0.01 p ∇ = , ignoring the fluid velocity inside the solid block (u = 0, v = 0), the four boundaries are all periodic boundaries; the initial concentration in the field is C0 = 0, the gas diffuses from the left, the concentration at the inlet is unity, the diffusion coefficient inside and outside the solid is consistent, the diffusion coefficient is 1/6, and the other boundaries are set to / 0 n C ∂ ∂ = . Three sets of simulations with different adsorption constants were carried out, respectively. The simulation results are shown in Figure 3b. The LB results fit well with the Langmuir isotherm adsorption curves.    The size of the second model is set to 100 × 100. There is a 40 × 40 solid block in the center of the flow field, which represents a cluster of coal particles. Assuming that the solid block is a homogeneous material, the internal micro-and nanopores adsorb gas molecules follow the Langmuir adsorption rate equation [42]. The fluid is driven by a pressure gradient, 0.01 p ∇ = , ignoring the fluid velocity inside the solid block (u = 0, v = 0), the four boundaries are all periodic boundaries; the initial concentration in the field is C0 = 0, the gas diffuses from the left, the concentration at the inlet is unity, the diffusion coefficient inside and outside the solid is consistent, the diffusion coefficient is 1/6, and the other boundaries are set to / 0 n C ∂ ∂ = . Three sets of simulations with different adsorption constants were carried out, respectively. The simulation results are shown in Figure 3b. The LB results fit well with the Langmuir isotherm adsorption curves.   The size of the second model is set to 100 × 100. There is a 40 × 40 solid block in the center of the flow field, which represents a cluster of coal particles. Assuming that the solid block is a homogeneous material, the internal micro-and nanopores adsorb gas molecules follow the Langmuir adsorption rate Energies 2020, 13, 4927 8 of 18 equation [42]. The fluid is driven by a pressure gradient, ∇p = 0.01, ignoring the fluid velocity inside the solid block (u = 0, v = 0), the four boundaries are all periodic boundaries; the initial concentration in the field is C 0 = 0, the gas diffuses from the left, the concentration at the inlet is unity, the diffusion coefficient inside and outside the solid is consistent, the diffusion coefficient is 1/6, and the other boundaries are set to ∂C/∂n = 0. Three sets of simulations with different adsorption constants were carried out, respectively. The simulation results are shown in Figure 3b. The LB results fit well with the Langmuir isotherm adsorption curves.

Diffusion -Adsorption of Gases in Simple Porous Media
In this paper, a simplified bi-dispersed porous media model is taken as an example to investigate the effects of fluid flow, gas diffusion (Fickian diffusion), and gas adsorption/desorption on the coupling process. As shown in Figure 4, the simulated area is 200 × 200, and the resolution of each grid is 10 nm. There are 16 clusters of coal particles with a length of 30 × 30 distributed in the field.   The fluid is driven by the pressure gradient, and the boundaries are set as the periodic boundary conditions. Besides, the left entrance of the field is set to the fixed concentration boundary (C = C 0 ), the other boundary condition is set to ∂C/∂n = 0, and the initial concentration in the field is 0. Since the extremely low permeability and low porosity of the coal seam, fluid velocity inside the clusters of coal particles is low enough to assume the clusters as an impermeable material (u = 0). Gas transport inside the clusters in the form of the Knudsen diffusion and gas-solid adsorption/desorption only occurs in micropores inner the clusters. The local adsorption amount and concentration update with time. Due to the interaction between the adsorptive sites and the gas molecules that do not cause a sharp change during the period, the numerical stability in the adsorption process can be ensured. The specific parameters in the simulation are shown in Table 1.     Figure 5 shows the gas concentration and adsorbed amount distribution of the five cases, respectively. We can observe that: (i) The influence of concentration distribution for fluid flow and gas diffusion is non-uniform in the field; a large number of gases accumulate near the inlet and the effect on the rest area is limited by the gas diffusivity. Gas adsorption mainly occurs on the first column of the solid particles, and the maximum adsorbed amount appears near the solid wall facing the inlet boundary. The adsorbed amount on the solid particles along the flow direction decreases gradually. (ii) When the fluid velocity is reduced by 10 times and the other conditions are maintained unchanged, the performance of gas diffusion-adsorption is almost identical to that of the case (i), which indicates that the fluid flow rate is not a main controlling factor. (iii) As the gas diffusion rate increases by 10 times only, the extent of gas adsorption effects is increased, and a more uniform concentration gradient is presented along the flow direction. The appearance of gas adsorption is similar to gas diffusion, which indicates that this is an adsorption-limited diffusion process. (iv) The Pe and Da number are both small (the FDC is increased by 10 times, other conditions are maintained unchanged), the gas diffusion rate is relatively large, the gas concentration distribution of the whole system is more consistent, the gas adsorption of the organic particles expands evenly from the boundary, and the adsorption amount of all particles is almost identical, which means the diffusion process is limited by the adsorption rate and the adsorption rate is low enough to keep the concentration field uniform at all times so that adsorption on all solid walls occurs uniformly. The results of cases (ii)-(iv) indicate that the magnitude of the FDC determines the distribution of the gas concentration and the location of the adsorption, which has a significant effect on the adsorption. (v) When the Pe number is small, but the Da number is large (which means the adsorption constant is increased by 10 times compared with the former case, other conditions are maintained unchanged), a higher inlet concentration and adsorption rate can accelerate the concentration update in the field, the adsorption intensity is higher, the adsorption occurs uniformly, and the adsorption amount increases remarkably.
By comparing the results of the case (i) and (ii), it can be found that the fluid velocity is not the main controlling factor in the methane migration and adsorption process at the pore-scale. The difference among the results of the cases (ii)-(iv) indicates that the magnitude of the FDC controls the extent of the gas diffusion, affecting the position and form of adsorption (whether uniform  Figure 5 shows the gas concentration and adsorbed amount distribution of the five cases, respectively. We can observe that: (i) The influence of concentration distribution for fluid flow and gas diffusion is non-uniform in the field; a large number of gases accumulate near the inlet and the effect on the rest area is limited by the gas diffusivity. Gas adsorption mainly occurs on the first column of the solid particles, and the maximum adsorbed amount appears near the solid wall facing the inlet boundary. The adsorbed amount on the solid particles along the flow direction decreases gradually. (ii) When the fluid velocity is reduced by 10 times and the other conditions are maintained unchanged, the performance of gas diffusion-adsorption is almost identical to that of the case (i), which indicates that the fluid flow rate is not a main controlling factor. (iii) As the gas diffusion rate increases by 10 times only, the extent of gas adsorption effects is increased, and a more uniform concentration gradient is presented along the flow direction. The appearance of gas adsorption is similar to gas diffusion, which indicates that this is an adsorption-limited diffusion process. (iv) The Pe and Da number are both small (the FDC is increased by 10 times, other conditions are maintained unchanged), the gas diffusion rate is relatively large, the gas concentration distribution of the whole system is more consistent, the gas adsorption of the organic particles expands evenly from the boundary, and the adsorption amount of all particles is almost identical, which means the diffusion process is limited by the adsorption rate and the adsorption rate is low enough to keep the concentration field uniform at all times so that adsorption on all solid walls occurs uniformly. The results of cases (ii)-(iv) indicate that the magnitude of the FDC determines the distribution of the gas concentration and the location of the adsorption, which has a significant effect on the adsorption. (v) When the Pe number is small, but the Da number is large (which means the adsorption constant is increased by 10 times compared with the former case, other conditions are maintained unchanged), a higher inlet concentration and adsorption rate can accelerate the concentration update in the field, the adsorption intensity is higher, the adsorption occurs uniformly, and the adsorption amount increases remarkably.
By comparing the results of the case (i) and (ii), it can be found that the fluid velocity is not the main controlling factor in the methane migration and adsorption process at the pore-scale. The difference among the results of the cases (ii)-(iv) indicates that the magnitude of the FDC controls the extent of the gas diffusion, affecting the position and form of adsorption (whether uniform adsorption). The results of the case (iv) and case (v) show that the magnitude of the adsorption constant determines the adsorption strength/rate of the clusters of the coal particles.

Diffusion-Adsorption of Gas in 2D Reconstituted Porous Media
Coal has a complex pore structure and is difficult to characterize. In this section, the coal matrix is reconstructed based on the Random Circles Packing (RCP) algorithm. The reconstruction process is controlled by the random growth kernel and porosity, and the reconstructed image is shown in Figure 6. Due to the existence of the nanopores in the clusters of the coal particles, the effect of adsorbed gas molecules on internal pores of the clusters (adsorption layer effect) should be considered, and other conditions are set in the same way as in Section 4.1.
Energies 2020, 13, x FOR PEER REVIEW 11 of 18 adsorption). The results of the case (iv) and case (v) show that the magnitude of the adsorption constant determines the adsorption strength/rate of the clusters of the coal particles.

Diffusion-Adsorption of Gas in 2D Reconstituted Porous Media
Coal has a complex pore structure and is difficult to characterize. In this section, the coal matrix is reconstructed based on the Random Circles Packing (RCP) algorithm. The reconstruction process is controlled by the random growth kernel and porosity, and the reconstructed image is shown in Figure 6. Due to the existence of the nanopores in the clusters of the coal particles, the effect of adsorbed gas molecules on internal pores of the clusters (adsorption layer effect) should be considered, and other conditions are set in the same way as in Section 4.1. From the results of Section 4.1, it can be found that the adsorption is limited by the FDC and the adsorption constant, both of which are related to the time. This section employed the adsorption saturation (Equation (8) to update the effective pore radius and investigated the impact of the effective KDC and adsorption constant on the adsorption process and compared the corresponding results with Xiong's model [35]. The simulation results can be seen in Figures 7-9.  Figure 7 shows the simulation results at 500,000 steps as follows: (a) the fluid velocity distribution, (b) the gas concentration distribution, and (c) the gas adsorption amount, respectively. It illustrates that the gas diffusion behavior inside and outside of clusters of the coal particles is quite different. Outside the clusters, the diffusivity is large and gas diffusion is relatively uniform. While inside the clusters, due to the limitation of the pore size, the diffusion process is very slow, leading From the results of Section 4.1, it can be found that the adsorption is limited by the FDC and the adsorption constant, both of which are related to the time. This section employed the adsorption saturation (Equation (8) to update the effective pore radius and investigated the impact of the effective KDC and adsorption constant on the adsorption process and compared the corresponding results with Xiong's model [35]. The simulation results can be seen in Figures 7-9.
Energies 2020, 13, x FOR PEER REVIEW 11 of 18 adsorption). The results of the case (iv) and case (v) show that the magnitude of the adsorption constant determines the adsorption strength/rate of the clusters of the coal particles.

Diffusion-Adsorption of Gas in 2D Reconstituted Porous Media
Coal has a complex pore structure and is difficult to characterize. In this section, the coal matrix is reconstructed based on the Random Circles Packing (RCP) algorithm. The reconstruction process is controlled by the random growth kernel and porosity, and the reconstructed image is shown in Figure 6. Due to the existence of the nanopores in the clusters of the coal particles, the effect of adsorbed gas molecules on internal pores of the clusters (adsorption layer effect) should be considered, and other conditions are set in the same way as in Section 4.1. From the results of Section 4.1, it can be found that the adsorption is limited by the FDC and the adsorption constant, both of which are related to the time. This section employed the adsorption saturation (Equation (8) to update the effective pore radius and investigated the impact of the effective KDC and adsorption constant on the adsorption process and compared the corresponding results with Xiong's model [35]. The simulation results can be seen in Figures 7-9.  Figure 7 shows the simulation results at 500,000 steps as follows: (a) the fluid velocity distribution, (b) the gas concentration distribution, and (c) the gas adsorption amount, respectively. It illustrates that the gas diffusion behavior inside and outside of clusters of the coal particles is quite different. Outside the clusters, the diffusivity is large and gas diffusion is relatively uniform. While inside the clusters, due to the limitation of the pore size, the diffusion process is very slow, leading  that the gas diffusion behavior inside and outside of clusters of the coal particles is quite different. Outside the clusters, the diffusivity is large and gas diffusion is relatively uniform. While inside the clusters, due to the limitation of the pore size, the diffusion process is very slow, leading to the big concentration difference, and thus an amount of gas accumulated near the interface. Consequently, the gas adsorption first occurs at the periphery of the particle, and then slowly advances inward. In the simulation, the impact of the fluid flow on gas diffusion and adsorption is negligible due to its weak effect at the pore-scale, which has been proved in Section 4.1. Therefore, only the influence of the Pe and Da are discussed in the following. For the sake of comparison, we defined an average diffusion coefficient of the clusters of coal particles to describe the global gas diffusion characteristics: where D e f f (x, y) is the effective KDC at a position (x,y) at a certain time and W sum is the total area of clusters in the region. Figure 8a shows the variation of the global diffusion coefficient (GDC) for nine groups of Pe and Da (at 11.2 million steps), respectively. The GDC of the organic particle is normalized by the inherent KDC, D OM,0 (which ignores the effect of the gas adsorption layer). It can be found that with the Pe increased, the local gas diffusivity drops due to the gas adsorption was limited by diffusion, and the decay rate of the GDC becomes slow. When Pe remains unchanged, the local adsorption/desorption rate is accelerated as Da is increased, resulting in the gas molecules quickly filled the nanopores and the process of gas adsorption-desorption reaching dynamic equilibrium, and the GDC decay rate becomes faster. From the trend of groups 4-9 (red and black lines), we can find that the decay of the GDC will eventually reach a unified value. Groups 1-3 (blue lines) will also reach this unified value but will take a relatively long time.
To further explore the effect of the adsorption/desorption rate on the adsorption process, Figure 8b separately extracted three sets of data at Pe = 7.5 × 10 −4 and compared with the results obtained by the method based on the Langmuir isotherm adsorption equation (the red line). Due to the time for the gas diffuse to the adsorptive site was ignored, the GDC decay rate of the Langmuir equation reaches the maximum value. However, sufficient time is required for gas to diffuse to an adsorptive site, even the pressure of the site reaches a specific value, the adsorption amount will not immediately reach the corresponding amount calculated by the Langmuir isothermal adsorption equation, which is just a maximum amount of gas that can be adsorbed under the pressure conditions. For a coal matrix with a fixed space size, due to the limited gas adsorption amount, the results obtained by the kinetic diffusion-adsorption-desorption model will eventually return to the analytical solution of the static isotherm adsorption equation with sufficient time; however, coal is usually tight and has low permeability, so it is difficult to ensure that the gas diffusion and adsorption are sufficient, the direct use of static isotherm adsorption equation may be incorrect. Figure 8c,d show the effects of p L and V m on diffusivity, respectively. It shows that the greater the p L or V m , the slower the GDC decay. This is because the p L is related to the ability of gas-solid adsorption and desorption: the greater the p L means the faster the increment of the desorption rate than the adsorption rate, and the greater the V m , the greater the adsorption capacity. The pore radius is a key parameter that affects the diffusion performance of gas in clusters. This section assumes that the pore size distribution inside the clusters conforms to the logarithmic normal distribution law, and we employed 6 groups of the typical pore size distribution of coal to investigate the dynamic effect of pore size on gas diffusion, as shown in Figure 9-the corresponding results can be seen in Figure 10. The pore radius is a key parameter that affects the diffusion performance of gas in clusters. This section assumes that the pore size distribution inside the clusters conforms to the logarithmic normal distribution law, and we employed 6 groups of the typical pore size distribution of coal to investigate the dynamic effect of pore size on gas diffusion, as shown in Figure 9-the corresponding results can be seen in Figure 10.   Figure 10 shows the GDC variation of clusters with different pore size distributions. It can be found that as the median or logarithmic standard deviation σ increases, the GDC decreases more slowly, and the final value of the steady-state is smaller. The main reason for this phenomenon is that the greater of the median or σ means a wider range of pore radius distribution, which is with a higher proportion of pores with larger pore diameters in the particles. The blocking effect of the adsorbed gas molecules on the pore space will be weaker.
In Section 4.1, the characteristics of gas flow-diffusion-adsorption with different characteristic parameters Re, Pe, and Da were visualized through a simple bi-dispersed porous media. In Section   Figure 10 shows the GDC variation of clusters with different pore size distributions. It can be found that as the median or logarithmic standard deviation σ increases, the GDC decreases more slowly, and the final value of the steady-state is smaller. The main reason for this phenomenon is that the greater of the median or σ means a wider range of pore radius distribution, which is with a higher proportion of pores with larger pore diameters in the particles. The blocking effect of the adsorbed gas molecules on the pore space will be weaker.
In Section 4.1, the characteristics of gas flow-diffusion-adsorption with different characteristic parameters Re, Pe, and Da were visualized through a simple bi-dispersed porous media. In Section  Figure 10 shows the GDC variation of clusters with different pore size distributions. It can be found that as the median or logarithmic standard deviation σ increases, the GDC decreases more slowly, and the final value of the steady-state is smaller. The main reason for this phenomenon is that the greater of the median or σ means a wider range of pore radius distribution, which is with a higher proportion of pores with larger pore diameters in the particles. The blocking effect of the adsorbed gas molecules on the pore space will be weaker.
In Section 4.1, the characteristics of gas flow-diffusion-adsorption with different characteristic parameters Re, Pe, and Da were visualized through a simple bi-dispersed porous media. In Section 4.2, the influence of different parameters on the gas diffusion-adsorption process is analyzed. This paper is a continuation of our previous studies [45], which studied gas transport in a coal reservoir with dynamic adsorption. The results further explain the effects of Langmuir pressure and volume, matrix porosity on the methane diffusion-adsorption process at pore-scale. What's more, Chen [46] investigated the impact of various parameters on the production of coalbed methane at the macroscale and found that gas-production rate increases with Langmuir pressure, matrix porosity, and desorption rate, which is consistent with the results of this paper and verifies the correctness of the results.
In the research, many key parameters of coalbed methane, such as Langmuir pressure and volume, cover a relatively wide range. Moreover, due to many characteristics, including geological environment, complex components, fractures, etc., of coal are excluded at pore-scale, the results obtained can be widely applied in the coalbed methane exploitation.

Conclusions
In this paper, we established the pore-scale LB model coupled with fluid flow, gas diffusion, and gas adsorption-desorption in the bi-dispersed porous media of coalbed methane. The Knudsen diffusion and dynamic adsorption-desorption of gas in clusters of coal particles were considered. Firstly, the model was verified by two classical cases. Then, three dimensionless numbers, Re, Pe, and Da, were adopted to discuss the impact of fluid velocity, gas diffusivity, and adsorption/desorption rate on the gas flow-diffusion-adsorption process. The effect of the gas adsorption layer in micropores on the diffusion-adsorption-desorption process was considered, and a Langmuir isotherm adsorption theory-based method was developed to obtain the dynamic diffusion coefficient, which can capture the intermediate process when adsorption/desorption reaches equilibrium. The pore-scale bi-disperse porous media of coal matrix was generated based on the RCP algorithm, and the characteristics of gas diffusion and adsorption in the coal matrix with different Pe, Da, and pore size distribution were discussed. The conclusions were as follows: (1) The influence of fluid velocity on the diffusion-adsorption process of coalbed methane at the pore-scale is very small and can be ignored; the magnitude of the FDC affects the spread range of gas diffusion and the process of adsorption and determines the position where adsorption takes place preferentially. (2) The magnitude of the adsorption constant controls the strength/rate of gas adsorption. A larger FDC or greater adsorption constant can effectively enhance the adsorption rate, and the trend of gas concentration-adsorption is closer to the Langmuir isotherm adsorption curve. (3) The gas diffusion-adsorption-desorption process is affected by the adsorption properties of coal.
The specific performance is that the greater the p L or V m , the slower the GDC decay. This is because the p L is related to the ability of gas-solid adsorption and desorption, the greater the p L means the faster the increment of the desorption rate than the adsorption rate, and the greater the V m , the greater the adsorption capacity. (4) The effect of the gas molecular adsorption layer has a great impact on the kinetic process of gas diffusion-adsorption-desorption. For a coal matrix with a fixed space size, due to the limited gas adsorption amount, the results obtained by the kinetic diffusion-adsorption-desorption model will eventually return to the analytical solution of the static isotherm adsorption equation with sufficient time; however, coal is usually tight and has low permeability, so it is difficult to ensure that the gas diffusion and adsorption are sufficient, the direct use of static isotherm adsorption equation may be incorrect.
In the research, many key parameters of coalbed methane, such as Langmuir pressure and volume, cover a relatively wide range. Moreover, due to many characteristics, including geological environment, complex components, fractures, etc., of coal are excluded at pore-scale, the results obtained can be widely applied in the coalbed methane exploitation.