**1. 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 a ffecting 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 di fficult 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 di ffuse out, which was related to the rapid production cycle of a gas productivity curve [6]. However, due to the limit of the gas di ffusion 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 grea<sup>t</sup> 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 di fficult 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–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

reach equilibrium. These above indicated that the time scale between the process of gas adsorption and gas di ffusion 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 di ffers 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–di ffusion 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 di ffusion–reaction flow in porous media [26–31].

The adsorption characteristics of coalbed methane in the reservoir are a ffected by many factors, such as fluid flow, gas transport, gas adsorption/desorption, etc. In this paper, a multi-component gas flow–di ffusion–adsorption coupled LB model was established to investigate the e ffects of fluid flow, gas transport, adsorption/desorption on adsorption, and adsorption-induced matrix di ffusion coe fficients and porosity.

#### **2. 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 di ffusion in the matrix, besides, the Knudsen di ffusion 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].

#### *2.1. The LB Equation for Fluid Flow*

The LB method for gas viscous flow at pore-scale can be expressed as follows [32]:

$$f\_l(\mathbf{x} + \mathbf{e}\_l \delta\_{l'} t + \delta\_l) - f\_l(\mathbf{x}, t) = -\frac{1}{\pi} [f\_l(\mathbf{x}, t) - f\_l^{eq}(\mathbf{x}, t)] \tag{1}$$

$$f\_i^{\text{eq}} = \rho \omega\_i [1 + \frac{\mathbf{e}\_i \cdot \mathbf{u}}{c\_s^2} + \frac{\mathbf{u} \mathbf{u} : (\mathbf{e}\_i \mathbf{e}\_i - c\_s^2 \mathbf{I})}{2c\_s^4}] \tag{2}$$

where *fi* is the discrete density distribution function, *feq i* is the local equilibrium function, *ei* is the discrete velocity of the particle, *cs* is the lattice sound velocity, *c* is the sound velocity, τ is the relaxation time, ω*i* is the weight coe fficient, and δ*t* is the time step.

#### *2.2. The LB Equation for Gas Di*ff*usion–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:

$$\mathcal{g}\_i(\mathbf{x} + \mathbf{e}\_i \delta\_t, t + \delta\_t) - \mathcal{g}\_i(\mathbf{x}, t) = -\frac{1}{\tau\_\mathcal{g}} [\mathcal{g}\_i(\mathbf{x}, t) - \mathcal{g}\_i^{\text{eq}}(\mathbf{x}, t)] + \omega\_i \mathcal{R}\_s \delta\_t \tag{3}$$

$$g\_i^{\alpha \eta} = \mathbb{C}\_s a \nu\_i [1 + \frac{\mathbf{e}\_i \cdot \mathbf{u}}{c\_s^2}] \tag{4}$$

where *gi* is a discrete concentration distribution function; *geqi* is the corresponding local equilibrium distribution function; <sup>τ</sup>*g* is the diffusion-related relaxation time; *Rs* is the source/sink term associated with gas adsorption/desorption.

For the standard LB equation, the macroscopic density and velocity can be defined as ρ = &*i fi*, ρ**u** = & *i fi***<sup>e</sup>***i*, the concentration can be defined as *Cs* = & *i gi*. Using the Chapman–Enskog technique to expand the LB equation of viscous flow, and with *p* = *<sup>c</sup>*2*s*ρ, μ = *<sup>c</sup>*2*s*(<sup>τ</sup> − 0.5)<sup>δ</sup>*t* the Equations (1) and (2) can recover the Navier–Stokes equation under the incompressible limit. Similarly, the LB equation of diffusion-reaction with *Ds* = *<sup>c</sup>*2*s*<sup>τ</sup>*<sup>g</sup>* − 0.5<sup>δ</sup>*<sup>t</sup>*, the Equations (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]:

$$D\_K = \frac{2r}{3} \sqrt{\frac{8RT}{\pi M}}\tag{5}$$

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:

$$r\_{eff} = r - d\_m \frac{p}{p + p\_L} \tag{6}$$

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 saturation is introduced by the ratio of the instantaneous adsorption amount to the saturated adsorption amount, and the e ffective pore radius is then considered:

$$r\_{eff} = r - d\_m \frac{V(\mathbf{x}, t)}{V\_m} \tag{7}$$

where *V* m and *V* are the saturated adsorption amount and the adsorption amount at a current time, respectively.

The modified e ffective porosity of the capillary can be written as:

$$
\varepsilon\_{eff} = \varepsilon\_0 \frac{r\_{eff}^2}{r^2} \tag{8}
$$

Due to the complex geometric composition of coal, the e ffects of porosity and tortuosity of the porous media also should be considered. Therefore, e ffective KDC can be given by:

$$D\_{eff} = \frac{\varepsilon}{\tau\_w} D\chi = \frac{\varepsilon\_0}{\tau\_w} \frac{r\_{eff}^2}{r^2} \frac{2r}{3} \frac{r\_{eff}}{r} \sqrt{\frac{8RT}{\pi M}} = \frac{\varepsilon\_0}{\tau\_w} \frac{r\_{eff}^3}{r^3} D\chi \tag{9}$$

where ε0 is the initial porosity of the coal matrix and τw is the tortuosity.

#### *2.3. 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:

$$\frac{\partial V}{\partial t} = k\_d \mathbb{C}(V\_m - V) - k\_d V \tag{10}$$

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:

$$V = \frac{V\_{mp}}{p + P\_L} = \frac{V\_m \mathcal{C}}{\mathcal{C} + \mathcal{C}\_L} \tag{11}$$

where *p* is the gas pressure, *p* = *CMc*<sup>2</sup> *s*, Pa; *C*L is the concentration corresponding to the Langmuir pressure, *CL* = *kd*/*ka*, mol/m3.

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:

$$R\_8 = \frac{\partial V}{\partial t} = k\_4 \mathbb{C}(V\_m - V) - k\_d V \tag{12}$$

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 di fference scheme of Equation (10):

$$V\_{t+1} - V = \delta\_l [k\_d \mathcal{C}(V\_m - V) - k\_d V] \tag{13}$$

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 di ffusion in macropores are governed by the double distribution LB equation; due to the pore size of the clusters being extremely small, the e ffect of Knudsen di ffusion and adsorption should be considered.

**Figure 1.** Schematic of the bi-dispersed porous structure of the coal matrix. Note that the scanning electron microscope (SEM) image of coal is copied from reference [39].

#### **3. 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–di ffusion 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 di ffusion coe fficient 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]:

$$D\_s \frac{\partial \mathcal{C}}{\partial \mathbf{n}} = k \mathcal{C} \tag{14}$$

The <sup>L</sup>évesque analytical solution after the adsorption has stabilized is:

$$\frac{L}{C\_0} \frac{\partial \mathcal{C}}{\partial \mathfrak{u}} = 0.854 (\frac{\mu\_{\text{max}} L^2}{\mathfrak{x} D\_s})^{1/3} \tag{15}$$

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.

**Figure 2.** Diagram of the convection–diffusion system with surface adsorption reaction.

**Figure 3.** Comparison of the present LB results and the analytical solution. (**a**) Comparison of the LB results and the analytical solution of Henry adsorption. (**b**) Comparison of the LB results and the analytical solution of the Langmuir isotherm adsorption.

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, ∇*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.

#### **4. Results and Discussion**

#### *4.1. Di*ff*usion –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.

**Figure 4.** The simplified 2D reconstructed geometry of the coal matrix.

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.


**Table 1.** The input parameters in the simulation.

Note: the relevant parameters in the table are taken from typical geological exploration data of coal reservoir in Qinshui Basin [43]; the adsorption and desorption constant is from reference [41,44].

For the sake of analysis, three dimensionless numbers, namely the Reynolds number (*Re*), Péclet number (*Pe*), and Damkohler number (*Da*), choose to quantitatively analyze the e ffects of fluid flow, gas di ffusion, and adsorption/desorption on the coupling process. This cannot only quantify the simulation results but also extends the results to other similar situations based on the characteristics of the dimensionless number. Where *Re* = *UL*/<sup>ν</sup>, *Pe* = *UL*/*Ds*, *Da* = *kdL*2/*Ds*, and *U*, *L* is the characteristic velocity and length of the system, respectively. In the simulation, the *Re* characterizes fluid flow, the *Pe* represents the relative proportion of advection and di ffusion, and the *Da* describes the relative time scale of adsorption and gas di ffusion in the system. In the present work, five kinds of combinations of *Re*, *Pe*, and *Da* are used to investigate the e ffects of the combination of fluid flow, gas di ffusion, and gas–solid adsorption in the adsorption process: (i) *Re* = 5.25 × <sup>10</sup>−2, *Pe* = 7.5 × <sup>10</sup>−2, *Da* = 1.5, (ii) *Re* = 5.25 × <sup>10</sup>−3, *Pe* = 7.5 × <sup>10</sup>−3, *Da* = 1.5, (iii) *Re* = 5.25 × <sup>10</sup>−3, *Pe* = 7.5 × <sup>10</sup>−4, *Da* = 0.15, (iv) *Re* = 5.25 × <sup>10</sup>−3, *Pe* = 7.5 × <sup>10</sup>−5, *Da* = 0.015, (v) *Re* = 5.25 × <sup>10</sup>−3, *Pe* = 7.5 × <sup>10</sup>−5, *Da* = 0.15. The simulation results are shown in Figure 5.

**Figure 5.** *Cont*.

**Figure 5.** Concentration and adsorbed amount magnitude distribution. (**a**) Gas concentration distribution and (**b**) adsorbed amount distribution.

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.

#### *4.2. Di*ff*usion–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.

**Figure 6.** The reconstructed bi-dispersed porous media based on the Random Circles Packing (RCP) algorithm.

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.** Simulation results at 500,000 steps. (**a**) The fluid velocity distribution, (**b**) the gas concentration distribution, and (**c**) the gas adsorption amount.

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 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:

$$D\_{\rm OM}(t) = \sum\_{\mathbf{x}} \sum\_{\mathbf{y}} D\_{eff}(\mathbf{x}, \mathbf{y}) / \mathcal{W}\_{\rm sum} \tag{16}$$

where *Deff*(*<sup>x</sup>*, *y*) is the effective KDC at a position (*<sup>x</sup>*,*y*) at a certain time and *Wsum* 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, *<sup>D</sup>*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−<sup>4</sup> 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.

**Figure 8.** Changes of the global diffusion coefficient (GDC) with the adsorbed layer. (**a**) The variation of the GDC for 9 groups of *Pe* and *Da*, (**b**) the variation of the GDC at *Pe* = 7.5 × <sup>10</sup>−4, (**c**) the effect of *p*L on diffusivity, and (**d**) the effect of *V*m on diffusivity.

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 9.** Pore size distribution with different median and logarithmic standard deviation.

**Figure 10.** The variation of the GDC with the different pore size distribution.

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 e ffect of the adsorbed gas molecules on the pore space will be weaker.

In Section 4.1, the characteristics of gas flow–di ffusion–adsorption with di fferent characteristic parameters *Re*, *Pe*, and *Da* were visualized through a simple bi-dispersed porous media. In Section 4.2, the influence of di fferent parameters on the gas di ffusion–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 e ffects of Langmuir pressure and volume, matrix porosity on the methane di ffusion–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.
