Here, the applicability of molecular dynamics (MD) to the study of coalbed methane geology is examined, specifically the role in visualizing the multimember behaviors accompanying adsorption/diffusion of gaseous molecules, gas-induced coal swelling, pyrolysis, and combustion. A review of MD applications to the study of geoscience is available (Cygan and Kubicki, 2018). Its applications to estimate the microscopic behaviors in coalbed methane geology have yet to be reviewed however. The earliest work for the application of MD aimed to reveal the dynamics damage of the materials exposed to radiation (Gibson et al., 1960), highlighting the non-equilibrium boundary conditions for particles with continuous forces. This preliminary work led to qualitative estimation of thermodynamic properties and to quantitative study of microscopic processes. In 1964, Rahman simulated liquid hydrogen with an MD method in periodic boundary conditions (PBCs) and surprisingly discovered that even a few particles represented the thermodynamic properties of the whole system (Rahman, 1964). Since then, MD simulation of condensed matter has become available, attracting worldwide physicochemical scholars to invest in this research work, taking advantage of the zero approximation nature, tracking particle trajectory, and accurate simulation results. The application of MD to visualize coal physico-chemical processes can date from 1992 as “The minimum-energy conformations obtained after MD calculations have provided insights into the types of bonding that provide the rigid structure observed for coals.” using Given, Wiser, Solomon, and Shinn’s molecular representations of coal in Carlson’s work (Carlson, 1992). Since then, the coal literature, over the last 29 years, has generated a surprisingly large number of MD applications in coal studies. The MD work using those predominantly classical molecular representations of coal slices remains state-of-the-art until advances in the analytical techniques as well as modeling software can be obtained (Mathews and Chaffee, 2012).

The jumps in computation power have also helped us to overcome some of the simulation challenges in model construction and calculation scale, and also aided scholar confidence in the microscopic representations, bridging the macroscopic characteristics and microscopic behaviors of coal. Figure 1A depicts the gradual emergence of journal publications utilizing MD for coal or coal-related research over the last ∼20 years, covering disciplines such as engineering, chemistry, energy fuels, environmental science ecology, and physics (Figure 1B). Although, this is likely an underestimation of the work since it was assessed from only the Web of Science Core Collection. It is expected that this simulation method will become more common for exploring coal behaviors. This paper provides a dedicated review of the history advances, and state-of-the-art availability for MD utility, along with applications in adsorption/diffusion of gaseous molecules, gas-induced coal swelling, pyrolysis, and combustion.

The molecular mechanism of microscopic interactions between CH₄/CO₂/H₂O and coal is the base to understanding the occurrence state of these gaseous molecules and the induced swelling (Xiang et al., 2014). No experimental technology has been found to operate perfectly, especially for the microscopic interactions of the molecular scale. The traditional experimentations have patent limitations to detect the molecular mechanism of microscopic interactions between the gaseous molecules and coal. MD simulation has been a powerful tool to reveal the relationship between the structure and properties of coal and the interaction mechanism of physical chemistry system compared with the various instrumental analysis technologies. The flow chart for the adsorption process was designed based on the minimum energy principle (Figure 4). Combining GCMC and MD simulation could provide the quantitative understanding of gaseous molecules adsorbing onto coal. The different contributions from van der Waals force, electrostatic force, and hydrogen bond force to the reduction in non-bonding energy suggest that the interaction of coal-CH₄ typically exhibits as physical adsorption (Song et al., 2017a). However, the interaction between coal and CO₂ is dominated by physical adsorption and incorporates weak chemical adsorption; for coal-H₂O, the case encounters the simultaneous existence of physical and chemical adsorption (Xiang et al., 2014).

The initial efforts of simulating the adsorption of gaseous molecules onto coal were based on perfect graphite or defective graphite (Thierfelder et al., 2011). Also, the methane adsorption of MD calculation has been reported to amplify the gas storage in open-ended single-walled carbon nanotubes, where the hybrid isotherm model of Langmuir and Sips was nicely in line with the calculated results (Shokri et al., 2010). The present implementation of MD over-binds about as much as bare DFT (density function theory) calculations under bind, but yields a meaningful adsorption height. Thierfelder et al. (2011) conducted methane adsorption on graphene from the first principles including dispersion interaction, and presented a microscopic picture of the reliability of approximate schemes in which van der Waals forces were limited. For quite a long time, the MD simulation of methane adsorption on graphite (0 0 0 1) was used as the platform to evaluate the gas behaviors onto coal in the molecular scale via state-of-the-art surface science analysis and computation approaches (Dash., 1975; Vidali et al., 1991; Bruch et al., 2007). Mosher et al. (2013) investigated methane adsorption in carbon-based 3D pore networks and found that the simulation estimates overpredict excess adsorption. This discrepancy is potentially induced from the limitation of the experimental model. Liu et al. (2013) simulated the CH₄ adsorption on N- and B-containing carbon models of coal and proposed that the physical adsorption of CH₄ on coal depends slightly upon the adsorption sites, orientation of CH₄, and the electronegativity of the dopant in coal. Hu et al. (2010) calculated the adsorption properties of CH₄ and CO₂ onto coal in an attempt to shed light on the adsorption at the atomistic level, demonstrating the successful utilization of MD in the study of the adsorption field. The MD simulation indicates the stronger adsorption affinity of CO₂ over CH₄ almost in all the published literature to date. Compared with the experimental-based characterization, MD can be understood as a “virtual experiment” at the molecular level, serving as an interface between laboratory experimentation and physicochemical theories (Liu et al., 2013). From the last decade on, MD work has been extended to more complex cases, for example, the contributions from water and oxidation on the adsorption process.

The coal macromolecular representation used in Hu et al. (2010) was the Wiser model and the MD results proved that the coal could absorb more CO₂ than CH₄ at the given temperature, providing an alternative technique to clarify the interactions between coal and adsorbates even for practically unreachable conditions. The general outcomes of MD are frequently consistent for different coal macromolecular representations, proving MD to be a promising tool incorporating the coincidences between the experimental and simulation results. Since the last decade, the coal geology community has witnessed a blowout increase in the generation of coal macromolecular representations with the gradually popularizing analysis and modeling technologies. The combined Monte Carlo and molecular dynamics is the most favored method to reveal the single and competitive adsorption of gaseous adsorbates onto coal. Zhang et al. (2014) conducted MC and MD simulations at a temperature ranging from 308 to 370 K to investigate the methane adsorption on dry and moist coal, providing a quantitative insight into the impacts of moisture and temperature on methane adsorption. The experimentations of Goodman et al. (2007) suggested good agreement of adsorption amount for pressure up to 8 MPa, however, for higher pressures, the results of different laboratories diverged distinctly in ranks, maceral, ash content, porosity, surface functional groups, temperature, and moisture.

Aiming at this issue, Liu and Wilcox (2011), Liu and Wilcox (2012a), Liu and Wilcox (2012b) conducted detailed simulation to articulate the CO₂ adsorption onto microporous carbons of heterogeneous surface functional groups, and the results suggested that the oxygen functional groups were favored for CO₂ adsorption, where the opposite effects can be indicated for the hydrated graphite. Based on the CS1000 coal model, the competitive adsorption behavior of CO₂ and CH₄ in coal micropores was also clarified based on the poromechanical model and MD and proposed a safe and long-term carbon sequestration resulting from the preferential adsorption of CO₂ to CH₄ (Brochard et al., 2012). Figure 5a∼fdepict the coal macromolecular representations which were frequently used in the publications. The consistent results among different macromolecular representations suggest the rationality of the MD method to depict the adsorption behaviors of gaseous adsorbates onto coal. Recently, MD has been prevalently used to optimize macromolecular representations prior to the adsorption process. Compared with the GCMC method, the adsorption process calculated from the MD simulation could provide more force and energy state parameters both for the initial and final status of the coal-gas compounds. However, MD requires much more computing memory and time than the MC calculation to acquire the normal termination status (Zhao et al., 2016a).

As increasingly dependable MD models are being developed for coal macromolecular representations, the adsorption processes of CH₄/CO₂/N₂/H₂O have been prevailing using the MD theory and method, providing a microscopic insight of high fidelity to present the structural heterogeneity of coal. In addition, the single maceral has been separated from the whole coal to represent the compositional complexity, as well as the macromolecular representations (Song et al., 2017a; Song et al., 2017b). For common gases such as CH₄, CO₂, H₂O, and N₂, the adsorption priority follows the order of CO₂> H₂O > CH₄> N_2, independent of the coal macromolecular representations used. For the outcomes of MD, Song et al. (2017a) exhibited the typical results and a comparison with the experimental results. As seen in Figures 6A–C, the maximum adsorption amounts of MD were close to the experimental result using the IS-100 isothermal adsorption desorption instrument (at 298 K), indicating the rationality of MD. However, their variation laws differ significantly, especially for the low pressure region. The adsorption amount of MD first increases sharply at low pressures and then remains stable with the increasing pressure. However, the experimental results gradually increase with the increasing pressure until 5 MPa. This was induced from the existence of ash in coal and may be the limitation of the coal macromolecular representations for the ash-related analysis. The D-A model fit better than the Langmuir model, which perhaps indicates the microporous filling properties for adsorption processes. For CO₂, the experimental adsorption amount was slightly lower than the MD results and their difference tended to decrease with the increasing pressure. When approaching the saturation state, the relative error of MD was reduced to 1.2%, maintaining an acceptable error (Song et al., 2018).

For the U.S., CO₂ emissions from fossil fuel combustion accounts for about 94% (Liu and Wilcox, 2011), and worldwide economic stability and development require energy, which today is also largely dominated by the direct combustion of fossil fuels (Markewitz et al., 2012; Middleton et al., 2012; Moore., 2012). As indicated in almost all the literature, coal has the advantage of adsorption of CO₂ over CH₄, thus, the injection of CO₂ into coal seams could not only improve the recovery of the coalbed methane (Xiang et al., 2014) but also minimize global emissions (Chen et al., 2021). This provides a promising continuable technique for carbon capture and storage. Although the recovery rate of CBM can be enhanced to 95% after injecting CO₂ compared with the original 77% (Reeves et al., 2004; Wang K. et al., 2020, Wang et al., 2021a), the adsorption of CO₂ onto coal could induce the differential swelling of the coal matrix, which will narrow (or even lead to the closure of the fracture network) the seepage channel and result in the reduction of permeability by ∼1 order of magnitude next to the injection well (Brochard et al., 2011; Wang R. et al., 2020).

The differential swelling of coal is caused by the reduction of mechanical strength after absorbing the gaseous adsorbates (Fan et al., 2020). For the last four decades, the swelling deformation of coal has been attracting increasing attention both via the experimental and mimetic angle. Green et al. (1984) firstly applied the volumetric swelling procedure to examine the coal swelling of the coal matrix and compared the results with the gravimetric swelling experiments, verifying the consistent results between the volumetric and gravimetric procedure. Reucroft and Sethuraman (1987) conducted dilatometric studies on swelling deformation induced from CO₂ adsorption and quantified that the swelling deformation induced from CO₂ adsorption can account for 20–50% of the specific surface area estimated from the low pressure CO₂ adsorption experiments. Aida et al. (1991) examined the steric effects on the solvent swelling of Illinois No. 6 coal and discovered good correlations between cross-link densities and swelling deformation behaviors. Then, Pan and Connell, 2007 established a theoretical model to depict adsorption-induced coal swelling and strain equilibrium, successfully describing the differences in swelling behaviors among different gas species. There also exist anisotropic properties for coal swelling, where more swelling occurs in the bedding direction than that parallel to the bedding (Pan and Connell, 2011). Larsen, 2014 proposed that the dissolved CO₂ serves as a plasticizer which can enable the physical structure rearrangements and make the coal into a more associated form in which fluids will be less soluble. Romanov et al. (2006) calculated the errors induced from coal swelling for CO₂ adsorption measurements and found that the true adsorption amount is at least 10% more than the observed values for the gravimetric experiment. For the swelling of pressure dependence, the swelling rate firstly increases with the increasing pressure until 8–10 MPa and remains stable (Day et al., 2008), indicating the existence of the swelling limit. The rank and sorbate used also impact the swelling rate in coals and these have been experimentally clarified by Astashov et al. (2008) and Majewska et al., 2010. Niu et al. (2017) estimated the adsorption swelling of natural and reconstituted anthracite coals using an experimental platform and discovered the swelling hysteresis phenomenon.

As mentioned above, the experimental methodology has achieved the quantification of the swelling rate and the influencing factors with the development of the physical analog. However, the experimentation could not reveal the microscopic mechanism of the swelling deformation in essence. MD was another efficient method to reproduce the swelling behaviors of gaseous adsorbates into coal. Narkiewicz and Mathews (2009) calculated the CO₂-induced swelling of coal via MD and proposed that the anisotropic swelling was greater perpendicular to the bedding plane. The deformation of coal induced by CH₄ adsorption was also quantified via the quenched solid density functional theory (QSDFT) and clarified the adsorption-induced deformation effect (Yang et al., 2010). Compared with MD, two qualitatively different types of swelling mechanisms can be identified depending upon the pore size. However, the QSDFT requires the minimization of the grand thermodynamic potential, which needs a lot of computation power. Based on the poromechanical model, Brochard et al. (2011) conducted effective MD work for swelling deformation behaviors of the CS1000 model induced from the adsorption of CO₂ and CH₄ in micropores, clarifying that the differential swelling was independent of the geological temperatures and pressures. This work highlights the foundation of the research for swelling deformation, and works as a guide for future coal scholars. Then, the solvent-coal interactions during solvent swelling were examined using the self-built macromolecular representations of Permian-aged South African coals and the results indicated that the nonbonding interactions play an important role during coal-solvent swelling. Xiang et al. (2014) conducted MD simulation for the competitive adsorption of CH₄/CO₂/H₂O and swelling of Yanzhou coal macromolecular representation. This GCMC and MD work clarifies that the reduction in valence electronic energy suggests the molecular mechanism of the swelling deformation for the “stressed coal” formed from the stress (Xiang et al., 2014). Using the extended poromechanical model developed by Brochard et al. (2011), the MC work for the Wiser model suggested that the swelling deformation enhances with the increasing mole fraction of CO₂ (Zhang et al., 2015). The calculation flow of the swelling deformation using GCMC and MD is depicted in Figures 8A–D using the coal vitrinite macromolecular representations built in our previous work (Song et al., 2017a). The adsorption process for the CO₂ molecules can be simulated via GCMC. Also, the pore volume, specific surface area, and pore size distribution can be calculated using the Atom Volumes & Surface Tool in the Materials Studio software. Then, the swelling ratio of coal vitrinite at this given condition was 5.4%, indicating the loss of porosity correspondingly. This provides a promising method to estimate the swelling deformation of coal and has achieved consistent outcomes (Song et al., 2019).

Zhao et al., 2016b proposed that the swelling deformation of coal was negligible at relatively low pressures and becomes significant with increasing pressure. We also calculated the volumetric change using the GCMC and MD based on the coal vitrinite macromolecular representation and proposed that higher temperature is conducive for the swelling equilibrium (Song et al., 2018), which is considered to be related to chemical thermodynamics. Frequently, the swelling rate, always following the order of [CO₂] > [CH₄] > [N₂], first increases sharply and then slowly when reaching a temperature >348 K independent of the adsorbate type and adsorption state. This reduction in the increase rate of swelling deformation was orientated from the structural variation under high temperatures. For the Materials Studio software environment, the swelling rate can be calculated via the difference in the configuration volume before and after MD. The research synchronously using the GCMC and MD suggests that the swelling rate is almost proportional to the adsorption quantity of CH₄/CO₂/H₂O before the adsorption limit (Song et al., 2018). The occurrence of the swelling deformation should overcome the non-bond potential energy, and thus the proneness of the swelling deformation can be estimated from the composition of molecular potential energy, which can be referred to in detail in Eq. 3. For the pressure dependence, the swelling rate increases slowly at low pressures and sharply at high pressures, in line with the results that higher pressures lead to the closure of the fractures and thus delay the gas diffusion speed (Zhao et al., 2016a; Song et al., 2018).

As the temperature and pressure dependence have been clarified, the differentiated swelling deformation of coal could be revealed in response to the geological conditions using the GCMC and MD (Chen et al., 2021; Wang H. et al., 2021, Wang et al., 2021c, Wang et al., 2021 J.). The work of Chen et al. (2021) has reached the quantification of differentiated swelling deformation for the binary mixture. For the binary mixture of n CO₂+ m CH₄ (binary mixture of N₂+ CH₄ with the molar ratio of n: m) and n N₂+m CH₄, the differentiated swelling deformation was dominated by the mole fraction of CO₂ and CH₄, respectively (Figure 9). It is noteworthy that, as the geological conditions were actually determined by the pressure and temperature, the variations of the swelling deformation for the coal matrix in response to the geological conditions is of great significance for the extraction of CBM and CO₂ capture and storage, especially for the deep un-minable coalbeds (Van Niekerk and Mathews, 2011). Based on the consistent results for the temperature and pressure dependence of swelling deformation, the variation of swelling deformation can be clarified using MD and a poromechanical model incorporating the geothermal gradient law. Chen et al. (2021) investigated the differential swelling of coal vitrinite at geological conditions and proposed that the differential swelling ratio increases regularly with the increase of the burial depth. Although both the swelling ratios for the binary mixtures of n CO₂+m CH₄ and the n N₂+m CH₄ increase with the increasing burial depth, there exist distinctly different variations between them (Figure 10). For the former, the differential swelling ratio first increases and then remains stable with the increasing bulk mole fraction of CO₂. However, for the latter, it gradually increases with the increasing bulk mole fraction of N₂ despite of the low increase rate at a higher bulk mole fraction of N₂. With the development of computational power and physical examination level, the simulation sophistication and improvement in MD, GCMC, and other numerical models will provide more opportunities to go beyond our current informed approach, gaining our confidence in the engagement in the estimation of coal swelling deformation behaviors.

There also exist properties for matrix shrinkage for coal. Durucan et al. (2009) conducted matrix deformation experiments on a range of different rank coals subjected to carbon dioxide and methane sorption at representative in situ reservoir pressures under unconfined conditions. These experiments were complemented by a set of simultaneous gas permeability and strain measurements, relating to coal rank and the mechanical and elastic properties. The matrix volume of coal shrinks when occluded gases desorb from its structure. In coalbed gas reservoirs, matrix shrinkage could cause the fracture aperture width to increase, causing an increase in permeability. Levine (1996) developed a computer model to evaluate the potential effect of matrix shrinkage on the absolute permeability and the results indicated that permeability changes as large as +250 mD are predicted for the upper case value of matrix shrinkage. Harpalani and Schraufnagel (1990) established gas pressure-volumetric strain relationships and proposed that the volume of the coal matrix shrinks by ≈0.4% when the gas pressure falls from 6.9 MPa to atmospheric pressure. The permeability increase as a result of matrix shrinkage was also observed in Bowen Basin coal where the matrix shrinkage becomes a dominant factor on cleat permeability over the effective horizontal stress. Liu and Harpalani (2013) reviewed the experimental and theoretical studies on coal matrix shrinkage and also developed a new theoretical approach to model sorption-induced coal shrinkage or swelling. Connell (2016) established a new interpretation of the response of coal permeability to changes in matrix shrinkage and found that considerable fines production at lower pore pressures with this reservoir could contribute to permeability damage, a mechanism not considered in the current model or the Shi-Durucan model.

As the first and necessary reaction step, coal pyrolysis, also termed as coal dry distillation, occurs at almost all the thermal upgrading processes of coal, producing CH₄, CO₂, H₂O, etc. (Hou et al., 2017; Wang et al., 2021a). Although the methodology for pyrolysis includes the physical method, chemical method, physical chemistry method, and computer-aided molecular design (CAMD), the traditional physical methodology faces great challenges when monitoring complex reaction processes and variations for free radicals (Khare et al., 2011; Ashraf et al., 2019). Mathews et al. (2011) proposed the combined experimental examinations and MD of ReaxFF to reveal the pyrolysis mechanisms of coal. ReaxFF MD can estimate the variations of free radicals in the pyrolysis process and adequately consider the weak interaction of chemical bond formation, Coulomb force, and van der Waals force. Therefore, the method of ReaxFF MD is more and more widely used in the coal pyrolysis process (Chenoweth et al., 2008; Zheng et al., 2013; Gao et al., 2018). The ReaxFF force field was proposed by Van Duin et al. (2001) to simulate the break and generation of the chemical bonds between atoms, thus systematically analyzing the continuous chemical reaction process involved in the pyrolysis process.

The simulation process for the ReaxFF MD includes 1) selection or construction of the coal macromolecular representations, 2) geometrical and energetic optimization for the coal macromolecular representations, and 3) simulation of the ReaxFF force field (Yang et al., 2020, Yang et al., 2021). The appropriate selection of the ReaxFF force field should enable the estimation and capture of the products and free radicals resulting from the pyrolysis reaction to obtain the product distribution and elementary reaction of the reaction (Kitsuka et al., 2007; Ashraf et al., 2019). The coal pyrolysis process can be simulated by ReaxFF MD, as well as the effects of pyrolysis temperature, which is of great significance for further study of the coal pyrolysis mechanism. ReaxFF MD can also be used to reproduce the large scale (above 1,000 of atoms) reactive system. For the pyrolysis/combustion process, ReaxFF MD can produce the results 100 times faster than the methods of quantum mechanics (QM) (Bhoi et al., 2014). Although kinetic chemical models have been established to depict the pyrolysis/combustion process of coal, they face utility challenges due to very expensive validation with experimental results and their time-consuming process. In light of this, MD simulation could provide a fundamental reaction mechanism based on the coal macromolecular representation. As mentioned above, the force field characterizing coal pyrolysis/combustion could depict the generation, transformation, and dissociation of the chemical bonds (Weismiller et al., 2010), whose results were in line with those from QM. Bhoi et al. (2014) conducted the MD simulation of spontaneous combustion using brown coal representation and proposed that the combustion of brown coal was initialized by thermal degradation at simulated temperatures high enough to enable chemical reaction within a computationally affordable time. Zheng et al. (2013) investigated the pyrolysis reaction of bituminous coal at 1,000–2000 K (250 ps in duration) based on the Wister and Shinn coal macromolecular representations and highlighted the consistency of the three types of pyrolysis products between MD and experimental results. The comparison of ReaxFF MD and thermogravimetric experimental results could verify the coal macromolecular representation constructed via the elemental analysis, nuclear magnetic resonance (NMR), and infrared spectroscopy (Wang et al., 2017).

The appropriate utilization of the force field enables us to clarify the complex reactions of hydrocarbons. As an extremely complex mechanism, the temperature of ReaxFF MD for the coal pyrolysis reaction mechanism was significantly higher than the actual reaction temperature, where the quantitative relation between the simulated and actual temperature was via the simple harmonic transition state theory. The background processing time of ReaxFF MD can be shortened by programming. Combined with the advantages of the ReaxFF reaction MD in exploring the mechanisms and behaviors of coal pyrolysis, the numbers of calculation cores and force field parameters should be furtherly optimized to improve the calculation efficiency, chemical analysis ability, authenticity, and conclusion accuracy of ReaxFF MD in simulating the pyrolysis mechanisms. Since the parameters of the reactive force fields (RFFs) were frequently derived from QM results, the reactive MD was a powerful tool for exploring complex reactive systems (Castro-Marcano et al., 2012). With the computational resources of the state of the art, the time scale required for ReaxFF MD simulation is of several orders of magnitude shorter than that required for the experiments (Castro-Marcano et al., 2012). Despite the traditional ReaxFF MD, the newly trained ReaxFF which incorporates sulfur-containing hydrocarbons has been found to correctly simulate pyrolysis and combustion dynamics of coal (Castro-Marcano et al., 2012), which provides an improved char construction methodology and demonstrates that the chemical events occurred during the combustion as well (Figure 11). ReaxFF MD simulation has been successfully utilized to reveal the atomic-scaled mechanism of the combustion of an Illinois No. 6 coal char combining the appropriate char representation and ReaxFF reactive force field. These reactive MD works have clarified the feasibility and capability of ReaxFF to describe the initial reactive events for coal pyrolysis and combustion.

Factors affecting the pyrolysis process such as rank, pulverized coal particle size, temperature, heating rate, pyrolysis atmosphere, residence duration, and catalyst could interact and mutually restrain, dividing the pyrolysis process into three stages, i.e., primary pyrolysis reaction, secondary pyrolysis reaction, and later condensation reaction (Kitsuka et al., 2007; Casal et al., 2018). The outcomes of the hydrocarbon yield with the temperature is depicted in Figure 12. There are three peaks during the pyrolysis process, which exhibit at 350–380°C, 380–470°C, and 480–520°C. There, the yields are low (C₁: 1.98–4.39 mg/g*TOC, C₂–C₅: 2.35–2.93 mg/g*TOC) (TOC, the content of the total organic carbon) for the first stage. The majority of C₂–C₅ in this stage comes from the initial cracking of aliphatic chains and oxygen functional groups. The instantaneous rate was also low at this stage, which was perhaps induced from limited heat accumulation to trigger the thermal cracking of coal macromolecules (Bhoi et al., 2014). For the second stage, production of C₁ and C₂–C₅ stably increased with temperature, where C₁ likely comes from the initial cracking of coal macromolecules but also from the secondary cracking of C₂–C₅ molecules. The instantaneous rates of C₂–C₅ decrease with the increasing temperature. The increasing yield of C₂–C₅ indicates that the generation rate of C₂–C₅ was higher than their cracking rates, which is consistent with the results of ReaxFF MD (Castro-Marcano et al., 2012). The third peak is from 480 to 520°C with a high hydrocarbon yield (C₁: 30.23–40.38 mg/g*TOC, C₂–C₅: 7.65–7.71 mg/g*TOC), as well as generation rates. In this stage, the faster rate of reactive dynamics time of MD enables us to initiate chemical reactions within a computational affordable time. The ReaxFF could detect the types of free radical intermediates to obtain the distribution and elementary reaction of general products during coal pyrolysis, which is of wider significance for explaining the mechanism of pyrolysis behaviors. The corresponding temperatures are representative of the loss of short chain aliphatic structures such as -CH₃ and -CH₂, and the condensation of aromatic structures.

The temperature programmed pyrolysis work is of great interest in coal utilization processes where there is a slow-heating rate. The pyrolysis products of vitrinite are mainly gaseous hydrocarbon C_1–5. Gaseous hydrocarbon yield and change law of functional groups were obtained. Light hydrocarbon (C₆–C₁₄) content of M1 (R_{o, max} = 0.58%), M2 (R_{o, max} = 0.80%), M3 (R_{o, max} = 1.59%), M4 (R_{o, max} = 1.61%), and M5 (R_{o, max} = 2.46%) calculated by experiment results of the closed system and open system are respectively 4 mg/g·TOC, 110 mg/g·TOC, 111 mg/g·TOC, 104 mg/g·TOC, and 38 mg/g·TOC. Thus, the yields of C₆–C₁₄ first increase and then decrease with the increasing rank. During pyrolysis, the higher carboxyl content and H-bonding were an essential contributor to the stability of coal macromolecules with the different lost rate of oxygen-containing functional groups and alkyl side chains with the increasing coal rank (Li et al., 2013).

As mentioned above, the main components of CBM were generated from the thermal evolution of the organic matters. Thus, the ReaxFF MD is an appropriate method to clarify the gas generation mechanism. The enrichment in coal is closely related to the occurrence of CH₄, which primarily can be divided into three states, i.e., adsorption, free, and dissolved states. These were related to the microscopic interactions between coal macromolecules and CH₄. Using MD we clarified the dominance of physical adsorption for CH₄ from different contributions of van der Waals force, electrostatic force, and hydrogen bond force. For CBM industry, CH₄ first desorbed from the micropores surface, then diffused into fractures of the coalbed, and finally was extracted through permeation flow. The hypothesis of the diffusion investigation of MD was the identification of the adsorption process, which assumes that the atom motions follow Newton’s first law and can be described via the empirical potential functions. Thus, MD has long been an appropriate approach to explore the underlying interaction mechanism of gas-coal and gas-gas during the diffusion process. MD could offer one approach to simulate the nanoscale diffusion process even for binary and ternary mixtures. An excellent review on the diffusion coefficients calculated from MD for binary and ternary mixtures is currently available (Liu et al., 2013). There also exist properties for matrix swelling and shrinkage for coal during the adsorption and desorption of CH₄. MD was another efficient method to reproduce the swelling behaviors of gaseous adsorbates into coal. Narkiewicz and Mathews (2009) calculated the CO₂-induced swelling of coal via MD and proposed that anisotropic swelling was greater perpendicular to the bedding plane. Based on the poromechanical model, Brochard et al. (2011) conducted effective MD work for swelling deformation behaviors of the CS1000 model induced from the adsorption of CO₂ and CH₄ in micropores, clarifying that the differential swelling was independent of the geological temperatures and pressures. These basic enrichment and transportation processes for coalbed methane geology are closely related to the interactions among coal and gaseous adsorbates, where MD can be appropriately used to reproduce these processes.

MD has applicability over a wide range of applications in coalbed methane geology. While much of the early work was categorized into “biological systems” that were of high performance and large-scaled, later work has been extended to quantify and rationalize coal’s complex behaviors. Compared with the experimental work, MD simulation can provide insights into the molecular mechanism of microscopic interactions between CH₄/CO₂/H₂O/N₂ and coal. With the development of the computation power, we are entering a period where MD will be used for higher level applications such as 1) adsorption of gaseous molecules onto coal, 2) diffusion of gaseous molecules into coal, 3) gas induced coal swelling, and 4) coal pyrolysis and combustion.

Combining GCMC and MD simulation can provide a quantitative understanding of gaseous molecules onto coal. MD simulation indicates the stronger adsorption affinity of CO₂ over CH₄ almost for all the published literature to date. Compared with the experimental-based characterization, MD can be understood as a “virtual experiment” at the molecular level, serving as an interface between laboratory experimentation and physicochemical theories. Since the last decade, the MD work has been extended to more complex cases, for example, the contributions from water and oxidation on the adsorption process. The hypothesis of the diffusion investigation of MD was identified to the adsorption process, assuming that the atom motions follow Newton’s first law and can be described via the empirical potential functions. Thus, MD has long been an appropriate approach to explore the underlying interaction mechanism of gas-coal and gas-gas during the diffusion process. The experimental observations face significant challenges when encountering the nanoscale diffusion process due to the heterogeneity of coal structure. Today, all types of the diffusion coefficients, such as the self, corrected, and transport diffusion coefficients can be calculated based on MD and the Peng-Robinson equation. To date, MD simulation of a single component has reached a situation of unprecedented success after Song et al. (2018). To the best of our knowledge, the current research origin is the most microscopically appropriate method to determine the diffusion properties in coal using the realistic coal model accounting for the heterogeneities of the physical and chemical structures in natural coal. Future efforts will clarify the impacts from the guest-host and guest-guest interactions on the transport properties in coal macromolecular representations.

Although the recovery rate of CBM can be enhanced after injecting CO₂ compared with the original 77%, the adsorption of CO₂ onto coal could induce the differential swelling of the coal matrix, which is caused by the reduction of mechanical strength after absorbing the gaseous adsorbates. For the last four decades, the swelling deformation of coal has been attracting increasing attention both via the experimental and mimetic angle. Based on the consistent results for the temperature and pressure dependence of swelling deformation, the variation of swelling deformation can be clarified using MD and a poromechanical model incorporating the geothermal gradient law. With the development of computational power and physical examination level, simulation sophistication and improvement in MD, GCMC, and other numerical models will provide more opportunities to go beyond our current informed approach, gaining our confidence in the engagement in the estimation of coal swelling deformation behaviors. Coal pyrolysis, also termed as coal dry distillation, occurs at almost all the thermal upgrading processes of coal. The faster rate of reactive dynamics time of MD enables us to initiate chemical reactions within an affordable computational time. The background processing time of ReaxFF MD can be shortened by programming. These reactive MD works have clarified the feasibility and capability of ReaxFF to descript the initial reactive events for coal pyrolysis and combustion. Combined with the advantages of the ReaxFF reaction MD in exploring the mechanisms and behaviors of coal pyrolysis, the representation construction method, auxiliary construction program, pyrolysis simulation temperature range, and pyrolysis time should be furtherly optimized to improve the calculation efficiency, chemical analysis ability, authenticity, and conclusion accuracy of ReaxFF MD in simulating the pyrolysis mechanisms. Advancing MD simulation (primarily characterized by a force field) will also allow the exploration of the more complex reaction process. The reaction mechanism of pyrolysis and spontaneous combustion should be a positive trend. The potential of MD for both visualization and microscopic mechanisms is also an interesting possible direction. Thus, it is expected that the availability of MD will continue to increase and will be added to the extensive list of advanced analytical approaches to explore the behaviors of multi-scale mechanisms in coalbed methane geology.