## 2 Materials and Methods  
### DFT Calculations  
We used plane-wave density functional theory (DFT) to predict the adsorption behaviors of methane and water on the coal surfaces at the molecular level. All calculations were performed using the DMol\({}^{3}\) program in the Materials Studio software package. (Version 2017 R2, Dassault Systemes BIOVIA, San Diego, 5005 Wateridge Vista Drive, CA, USA) [21; 22]. The Perdew-Burke-Ernzerhof (PBE) [23] form of the generalized gradient approximation (GGA) was used as the exchange-correlation functional. [PERSON]'s DFT-D method was adopted to account for the dispersion effect [24]. We also tested other forms of functionals, including revised PBE (RPBE) [25] and Becke-Lee-Yang-Parr (BLYP) [26] forms of GGA functionals, as well as two typical higher-level meta-GGA functionals, namely m11-L [27] and revtpss [28], but found that PBE coupled with DFT-D can closely reproduce the well-established adsorption energy value of methane on the surface of graphene, while keeping the computational cost reasonably low. A detailed comparison between functionals is given in Section 3.1.1. The triple numerical plus polarization (TNP) basis set (version 4.4) was used in all DFT calculations [21]. All other parameters related to the electronic and geometric steps were set to the default values at the fine level provided by the DMol\({}^{3}\) program. Coal is a chemical complex consists of C, O, H, N, and S. It does not have a fixed molecular structure and can be subdivided into different types, including anthracite coals, bituminous coals, subbituminous coals, and lignite coals, depending on the concentration of the C element [29]. In this study, we focused on bituminous coals and selected three classical molecular models, which are shown in Figure 1a-c, to simulate the coalbed. In this paper, these three models in Figure 1a-c are referred to as Given type 1 model (G1) [30], Given type 2 model (G2) [31], and Fuchs-Sandoff (FS) model [32], respectively, according to their original references. We chose these three models because they contain a wide range of structural elements and functional groups, including but not limited to aromatic planes, alkane chains, hydroxyl groups, carbonyl groups, and pyridine. These models can represent different structures of coal molecules. As shown later in Section 3.1, although our calculations were based on these three bituminous coals molecules, the results are general and can be extended to other types of coals. lower part of the amorphous cell was separated by 20 A, 60 A and 100 A, respectively, to create a gap which represents micropores in the coalbed where CH\({}_{4}\) and H\({}_{2}\)O can be adsorbed. Methane and water were represented by full atomic models. The COMPASS forcefield was used in this work [33]. Electrostatic interactions were calculated by using forcefield assigned charges. The summation methods for electrostatic and van der Waals interactions were Ewald and atom-based, respectively. Grand canonical Monte Carlo (GCMC) simulations were used to calculate the adsorption isotherms of CH\({}_{4}\), H\({}_{2}\)O and the mixture of CH\({}_{4}\) and H\({}_{2}\)O at various pressures in coalbed with different pore sizes. The temperature was fixed at 318.15 K (45 \({}^{\circ}\)C) since this is a typical coalbed temperature in coalbed gas extraction processes. For each case, 10\({}^{4}\) configurations at equilibrium were generated and the configuration with the lowest energy was used as the initial input in a MD simulation. The MD simulations were performed in the constant-temperature, constant-volume (NVT) ensemble with a Nose-Hoover thermostat to calculate the diffusivities [34]. The cutoff radius of 15.5 A and time step of 1 fs were used. During simulations, the coalbed model atoms were kept fixed. The diffusion coefficients were calculated by the mean squared displacement (MSD) based on the [PERSON] relation given by Equation (1).