Molecular Simulation Study on Methane Adsorption in Amorphous Shale Structure

: Gas adsorption in the porous shale matrix is critical for gas-in-place (GIP) evaluation and exploration. Adsorption investigations beneﬁt signiﬁcantly from the use of molecular simulation. However, modelling adsorption in a realistic shale topology remains a constraint, and there is a need to study the adsorption behaviour using molecular models containing both organic and inorganic nanopores. Most simulations use a single component, either kerogen (organic composition) and quartz or clay (inorganic composition), to represent the shale surface. In this work, the molecular dynamic (MD) and grand conical Monte Carlo (GCMC) simulations were utilised to provide insight into methane adsorption behaviour. Amorphous shale structures composed of kerogen and quartz were constructed. The kerogen content was varied to replicate the shale with 2 wt.% and 5 wt.% Total Organic Carbon (TOC) content with 5 nm pore size. The simulated densities of the shale structures showed consistent values with actual shale from the Montney, Antrim, and Eagle Ford formations, with 2.52 g/cm 3 and 2.44 g/cm 3, respectively. The Average Error Analysis (ARE) was used to assess the applicability of the proposed amorphous shale model to replicate the laboratory adsorption isotherm measurements of actual shale. The ARE function showed that the amorphous shale shows good agreement with experimental measurements of all Barnett shale samples with an average of 5.0% error and slightly higher for the Haynesville samples with 8.0% error. The differences between the experimental adsorption measurement and simulation resulted from the amorphous packing, and actual shales have more minerals than the simulated model.


Introduction
In recent years, the increasing demand for cleaner energy to sustain the global economy has attracted significant interest in the exploration and development of natural gas. As a result, the Energy Information Administration (EIA) estimated approximately 7300 tcf of recoverable shale gas, and 32% of the total natural gas resources are in shale formations [1]. Shale gas comprises a mixture of four occurring gases: Methane, which makes up 70%-90% of natural gas, followed by ethane, butane, and propane [2]. Unlike conventional reservoirs, shale formations have low porosity and nano-Darcy (nD) permeability characteristics [3]. Shale rocks are inherently heterogenous, and the shale matrix pores can be divided into organic and inorganic pores [4]. The pore structure of shale is highly complex, with different sizes classified as micropores (<2 nm), mesopores (2 nm-50 nm), and macropores (<50 nm) [5].
There are three ways gas was stored in the shale-free gas in pores and natural fractures, namely, adsorbed on the surface of organic and clay pore walls and absorbed in organic matter and connate water [6]. The adsorbed gas in shale accounts for up to 85% of the Original Gas in Place (OGIP) [7,8]. Adsorption in shale occurs through the mechanism of physical adsorption or physisorption [9]. Adsorption is caused by van der Waals and electrostatic contact forces between gas molecules and the shale surface. In addition, adsorption is facilitated by the abundance of nano-microscale pores that provide a large internal surface area and restrictions that force the fluid-solid phase close enough for a weak dipole-dipole interaction. In these nanopores, the phase behaviour of methane deviates significantly from conventional reservoirs [10].
Furthermore, in such microscopic scales, the effect of the walls on the molecules becomes more pronounced, inducing different behaviour of the fluid than from the bulk state in the absence of the walls [11,12]. The nanopore confinement effects include capillary pressure, a fluid-wall interaction, and molecular adsorption. Therefore, the result of adsorption is more significant inside shale nanopores. Therefore, a production forecast without adsorbed gas will be less accurate when the effect of gas adsorption is ignored.
Gas adsorption on shale matrix particles belongs to physical adsorption. Consequently, the organic matter in the shale matrix is a critical parameter that influences gas adsorption characteristics in shale gas reservoirs. The pore structure of organic matter in shale is predominately micropores and mesopores, with a total pore volume of 77% to 92%, providing a large specific surface area for the gas to be adsorbed [13][14][15][16].
Shale is a dense, fine-grained, laminated sedimentary rock that contains organic and inorganic matter. It should be highlighted that extensive experiments have observed the presence of intercalated organic and inorganic nanocomposites of the mineral matrix in shale [17,18].
A realistic shale organic matrix is composed of kerogen molecules, possessing complex amorphous structures and different surface attributes, inevitably impacting surfacemethane interactions and methane adsorption behaviour. Gonciaruk et al. (2021) revealed that kerogen is in an omnipresent amorphous phase [19]. A study by Xin et al. (2012) showed the kerogen structure in the Huadian shale model has very little polycyclic aromatic structure [20]. Most aromatic units are separated by various bridge bonds and long aliphatic chains, almost without a co-planar in space. In a separate work, Tong et al. (2011) also found that most methylene straight chains in Huadian kerogen cannot form a crystalline but amorphous structure [21].
Recent advances in computational modelling have aided in developing a reliable molecular computational model for kerogen. The realistic kerogen model was first proposed by Ungerer et al. (2015) [22]. Using experimental data, they used molecular dynamics and quantum mechanics to analyze different kerogen types (based on their biogenic origin). The second approach is to use molecular dynamics to generate the amorphous model. Using the reactive force field, disordered carbons were transformed into an amorphous solid structure in this method. The models have been referred to and reused in various studies, including investigating gas adsorption behaviour and kerogen swelling [19,23,24].
Studies on the adsorption behaviour of gas in slit-like amorphous silica nanopores with the GCMC method have been of huge interest among academicians from various fields. The adsorption behaviours within amorphous silica are more complicated than those in a well-defined crystalline structure, reflecting the heterogenous surfaces observed with shale [25]. Recent work has also seen a number of simulations on methane adsorption in amorphous silica to simulate adsorption in inorganic shale compared to a crystalline quartz structure [26][27][28].
Despite the rising scientific interest in exploring the adsorption characteristic of methane in unconventional formations with tight porous networks, there is a limited calibrated understanding of the influence of the chemistry of solid surfaces and nanoscale confinement in shale [29][30][31][32]. For example, some researchers on shale gas adsorption have explored the role of minerals and kerogen-interfaced partitioning on gas adsorption and nanoscale transport behaviour [31,33]. However, the minerals and kerogens were positioned at different slit ends, indicating that neither component is mixed and interacting. Notably, the copresence of amorphous organic and inorganic interfaces was rarely discussed [33]. The ultimate goal is to model amorphous shale with kerogen in the proximity of all the present minerals through molecular dynamics simulations.
Molecular dynamic simulation has been used to determine the interaction of methane and shale at a molecular level. However, the molecular simulation was usually simplified with adsorption measurements conducted with kerogen cells, and long computational work is required. Hence, a molecular dynamics simulation with an intrinsic shale composition arrangement is needed.
This work uses a molecular dynamics simulation with a new proposed intrinsic shale composition arrangement by varying kerogen cells following the TOC composition and incorporating quartz as inorganic components. The adsorption isotherm measurements were conducted with GCMC simulation to determine the applicability of this new proposed model to represent the actual adsorption capacity in shale. The results were compared with the actual shale adsorption isotherm measurements conducted experimentally in the literature.  [22,34,35]. The molecular structure is shown in Figure 1. discussed [33]. The ultimate goal is to model amorphous shale with kerogen in the proximity of all the present minerals through molecular dynamics simulations. Molecular dynamic simulation has been used to determine the interaction of methane and shale at a molecular level. However, the molecular simulation was usually simplified with adsorption measurements conducted with kerogen cells, and long computational work is required. Hence, a molecular dynamics simulation with an intrinsic shale composition arrangement is needed.

Materials and Methods
This work uses a molecular dynamics simulation with a new proposed intrinsic shale composition arrangement by varying kerogen cells following the TOC composition and incorporating quartz as inorganic components. The adsorption isotherm measurements were conducted with GCMC simulation to determine the applicability of this new proposed model to represent the actual adsorption capacity in shale. The results were compared with the actual shale adsorption isotherm measurements conducted experimentally in the literature.

Kerogen
The organic models used in this work rely on the kerogen monomers with the chemical formula C242H219013N5S2 initially developed by Ungerer et al. (2015) [22]. This kerogen is a Type II-C kerogen belonging to the oil/gas condensate window. The construction of the kerogen models derives from the work of Yiannourakou [22,34,35]. The molecular structure is shown in Figure 1.

Quartz
In shale, the organic matter is primarily kerogen, while the inorganic components are minerals. The minerals are typically silicates, carbonates, and pyrites [17,36]. In this work, quartz minerals represent the inorganic comp. Quartz or silica oxide (SiO2) naturally occurs in crystalline and amorphous forms. Quartz can occur in four different polymorphs (α, β, tridymite, and cristobalite) dependent upon its surrounding environment [17]. The most abundant form of silica is α-quartz, and the term quartz is often used in place of the general term crystalline silica. The α-quartz crystal structure from the Material Studio software packages structure database was used to construct the quartz unit model ( Figure  2) [37].

Quartz
In shale, the organic matter is primarily kerogen, while the inorganic components are minerals. The minerals are typically silicates, carbonates, and pyrites [17,36]. In this work, quartz minerals represent the inorganic comp. Quartz or silica oxide (SiO 2 ) naturally occurs in crystalline and amorphous forms. Quartz can occur in four different polymorphs (α, β, tridymite, and cristobalite) dependent upon its surrounding environment [17]. The most abundant form of silica is α-quartz, and the term quartz is often used in place of the general term crystalline silica. The α-quartz crystal structure from the Material Studio software packages structure database was used to construct the quartz unit model ( Figure 2) [37].

Construction of Shale Structure
The shale structures were generated using Amorphous Cell and Forcite modules in DS BIOVIA Material Studio 2020 software packages [37]. The shale matrix is generated by performing MD simulations in the conical ensemble (NVT) and the isobaric isothermal

Construction of Shale Structure
The shale structures were generated using Amorphous Cell and Forcite modules in DS BIOVIA Material Studio 2020 software packages [37]. The shale matrix is generated by performing MD simulations in the conical ensemble (NVT) and the isobaric isothermal ensemble (NPT). Initially, the structures of kerogen units are relaxed by geometry optimisation and annealing simulation [20]. A SMART minimisation algorithm with a fine convergence criterion optimised the structure geometry. The non-bonded interactions, which include Coulomb and Van der Waals interactions, are calculated using an atom and a fine cut-off distance of 12.5 Å. Ten annealing cycles with increasing temperatures from 298.15 K to 1200 K were performed for the annealing method. The simulations are run using the canonical ensemble (NVT) with a fixed molecular number, box volume, and system temperature. To bring the structures to the lowest energy state, a total simulation time of 400 ps is used at each stage. The simulation time for each relaxation stage is determined by referring to previously reported work by   [38,39]. The kerogen molecule's geometry optimization and annealing process flow chart is shown in Figure 3.

Construction of Shale Structure
The shale structures were generated using Amorphous Cell and Forcite modules i DS BIOVIA Material Studio 2020 software packages [37]. The shale matrix is generated b performing MD simulations in the conical ensemble (NVT) and the isobaric isotherma ensemble (NPT). Initially, the structures of kerogen units are relaxed by geometry optimi sation and annealing simulation [20]. A SMART minimisation algorithm with a fine con vergence criterion optimised the structure geometry. The non-bonded interactions, whic include Coulomb and Van der Waals interactions, are calculated using an atom and a fin cut-off distance of 12.5 Å. Ten annealing cycles with increasing temperatures from 298.1 K to 1200 K were performed for the annealing method. The simulations are run using th canonical ensemble (NVT) with a fixed molecular number, box volume, and system tem perature. To bring the structures to the lowest energy state, a total simulation time of 40 ps is used at each stage. The simulation time for each relaxation stage is determined b referring to previously reported work by   [38,39 The kerogen moleculeʹs geometry optimization and annealing process flow chart is show in Figure 3. The molecular modelling framework has the advantage of controlling the petrophys ical properties of kerogen during the initiation phase, where different kerogen types bonding strengths, and other chemical properties can be changed to isolate their impac independently. A realistic shale structure is required to serve as nanoporous media. I this research, the molecular modelling approach was followed to construct a 3D nanopo rous model. To develop nanoporous shale structures, first, 5-10 units of optimised kero gen molecules and 1399 units of quartz were randomly placed in a relatively large ce with periodic conditions to generate the initial configurations of shale models, with th target density being 0.1 g/cm 3 .
The initialisation of the system occurred at a low density to avoid any stability issues A sequence of MD simulations relaxes these configurations. A single stage of the iso choric-isothermal NVT ensemble is adopted to relax the configuration at 900 K for 250 ps Then subsequent NPT (a fixed molecular number, system temperature, and pressure) sim ulations are performed at 20 MPa with a stepwise decreasing temperature from 900 K t 350 K. The simulation time increases from 250 ps, sufficient for the shale models' densit convergence. The targeted temperature and pressure of the final configuration were se lected to be 350 K and 20 MPa, respectively, to represent typical reservoir conditions. Th time step for all MD simulations is 1 fs. The gradual decrease in temperature to the fina targeted level assured structural stability during convergence. The final structures ha densities of 2.3-3.14 g/cm 3 , with approximately ± 5 % differences from the values of quart The molecular modelling framework has the advantage of controlling the petrophysical properties of kerogen during the initiation phase, where different kerogen types, bonding strengths, and other chemical properties can be changed to isolate their impact independently. A realistic shale structure is required to serve as nanoporous media. In this research, the molecular modelling approach was followed to construct a 3D nanoporous model. To develop nanoporous shale structures, first, 5-10 units of optimised kerogen molecules and 1399 units of quartz were randomly placed in a relatively large cell with periodic conditions to generate the initial configurations of shale models, with the target density being 0.1 g/cm 3 .
The initialisation of the system occurred at a low density to avoid any stability issues. A sequence of MD simulations relaxes these configurations. A single stage of the isochoricisothermal NVT ensemble is adopted to relax the configuration at 900 K for 250 ps. Then subsequent NPT (a fixed molecular number, system temperature, and pressure) simulations are performed at 20 MPa with a stepwise decreasing temperature from 900 K to 350 K. The simulation time increases from 250 ps, sufficient for the shale models' density convergence. The targeted temperature and pressure of the final configuration were selected to be 350 K and 20 MPa, respectively, to represent typical reservoir conditions. The time step for all MD simulations is 1 fs. The gradual decrease in temperature to the final targeted level assured structural stability during convergence. The final structures had densities of 2.3-3.14 g/cm 3 , with approximately ±5% differences from the values of quartz shale samples measured experimentally [40]. A summary of the molecular instruction protocol is given in Figures 4 and 5. shale samples measured experimentally [40]. A summary of the molecular instruction protocol is given in Figures 4 and 5.

Model Validation
The simulation results depend primarily on the molecular model. In this work, the amorphous shale matrix was constructed to approximate actual shale, so the rationality of the model must be verified. shale samples measured experimentally [40]. A summary of the molecular instruction protocol is given in Figures 4 and 5.

Model Validation
The simulation results depend primarily on the molecular model. In this work, the amorphous shale matrix was constructed to approximate actual shale, so the rationality of the model must be verified.

Model Validation
The simulation results depend primarily on the molecular model. In this work, the amorphous shale matrix was constructed to approximate actual shale, so the rationality of the model must be verified.
In this work, the density of the amorphous model was compared with the densities of shale from three different formations, Montney, Antrim, and Eagle Ford, with TOC contents of 2 wt.% and 5 wt.% [41,42]. Physical density is a significant indicator of the rationality of the molecular model. Therefore, the simulated densities of shale were compared with available test data from the literature. The discrepancies between the simulated and experimental values are shown in Section 3.

Simulation Methods
The shale pore system is found to exist at multiple scales, with its size ranging from a few nanometres up to hundreds of micrometres. Capillary pressure analysis conducted on Barnett shale by Loucks et al. (2009) [43] concluded that the size of most pores is in the 5-15 nm range. Typical pore throat sizes of siliciclastic shale taken from Nelson et al.
(2009) vary from 5 nm to 100 nm [44]. A multiscale characterisation analysis on Longmaxi shale samples by Chen et al. (2021) showed the shale pore sizes range from almost 5 nm to 100 µm, with a pore size smaller than 300 nm representing the nanopore system [45].
The previous study has shown that the GCMC simulations of methane adsorption in kerogen slit suggested that gas adsorption capacity is closely related to the slit aperture and the average density of confined methane [46]. In this work, we used the slit shape pore size of 5 nm to conduct the adsorption isotherm measurements.
The Visualiser module of Material Studio software was used to build the slit pore models. First, the shale matrix was stacked parallel to create a slit-shaped pore. Next, the pore height (H) was determined by simulating a vacuum to obtain pore sizes of 5 nm [47]. Schematic representations of the shale slit pore structures are presented in Figure 6. In this work, the density of the amorphous model was compared with the densities of shale from three different formations, Montney, Antrim, and Eagle Ford, with TOC contents of 2 wt.% and 5 wt.% [41,42]. Physical density is a significant indicator of the rationality of the molecular model. Therefore, the simulated densities of shale were compared with available test data from the literature. The discrepancies between the simulated and experimental values are shown in Section 3.

Simulation Methods
The shale pore system is found to exist at multiple scales, with its size ranging from a few nanometres up to hundreds of micrometres. Capillary pressure analysis conducted on Barnett shale by Loucks et al. (2009) [43] concluded that the size of most pores is in the 5-15 nm range. Typical pore throat sizes of siliciclastic shale taken from Nelson et al. (2009) vary from 5 nm to 100 nm [44]. A multiscale characterisation analysis on Longmaxi shale samples by Chen et al. (2021) showed the shale pore sizes range from almost 5 nm to 100 μm, with a pore size smaller than 300 nm representing the nanopore system [45].
The previous study has shown that the GCMC simulations of methane adsorption in kerogen slit suggested that gas adsorption capacity is closely related to the slit aperture and the average density of confined methane [46]. In this work, we used the slit shape pore size of 5 nm to conduct the adsorption isotherm measurements.
The Visualiser module of Material Studio software was used to build the slit pore models. First, the shale matrix was stacked parallel to create a slit-shaped pore. Next, the pore height (H) was determined by simulating a vacuum to obtain pore sizes of 5 nm [47]. Schematic representations of the shale slit pore structures are presented in Figure 6. Monte Carlo (MC) simulations following the Metropolis scheme were employed to mimic the migration of gas molecules into a solid structure. This work employed the BIO-VIA Materials Studio 2020 simulation package for all MC simulations. Several different MC moves were permitted during simulations, mimicking the behaviour of actual gas molecules undergoing sorption into a solid lattice. The equilibration process is achieved by performing insertion, deletion, and translation moves for the methane molecules. For methane molecules, rotational moves are also applied.
Comparisons of the properties of methane from experimental measurements and simulation were gathered from the Energy calculation through the VAMP module in Material Studio software. The differences between simulated and experimental properties of methane are tabulated in Table 1. Materials Studio VAMP can rapidly predict many physical and chemical properties for molecular organic and inorganic systems using a semiempirical molecular orbital method and is an ideal intermediate approach between forcefield and first-principles methods. Monte Carlo (MC) simulations following the Metropolis scheme were employed to mimic the migration of gas molecules into a solid structure. This work employed the BIOVIA Materials Studio 2020 simulation package for all MC simulations. Several different MC moves were permitted during simulations, mimicking the behaviour of actual gas molecules undergoing sorption into a solid lattice. The equilibration process is achieved by performing insertion, deletion, and translation moves for the methane molecules. For methane molecules, rotational moves are also applied.
Comparisons of the properties of methane from experimental measurements and simulation were gathered from the Energy calculation through the VAMP module in Material Studio software. The differences between simulated and experimental properties of methane are tabulated in Table 1. Materials Studio VAMP can rapidly predict many physical and chemical properties for molecular organic and inorganic systems using a semi-empirical molecular orbital method and is an ideal intermediate approach between forcefield and first-principles methods. The grand conical ensemble (µVT), in which there is the chemical potential of gas (µ), volume (V), and temperature (T), describes a system in open contact with a reservoir. The grand conical ensemble is close to the adsorption case by performing a series of insertions and deletions of adsorbate gas molecules. In the ensemble, the µ and the T are fixed, and the equilibrium status is that the adsorbate molecules have the same chemical potential and temperature in the reservoir. In the Metropolis scheme, a trial configuration is compared to the original configuration.
The "sorption isotherm" module in Sorption of Material Studio with Grand Canonical Monte Carlo (GCMC) and a COMPASS forcefield was used to simulate the adsorption capacity of gases on the shale molecular cells [51]. Methane adsorption on shale models is simulated at 90°C (363.15 K) and 120°C (393.15 K) with pressure up to 30 MPa. The Anderson thermal bath was used to keep the temperature constant, while the cut-off radius was set as 12.5 Å. The configuration of adsorption isotherm measurements is described in Table 2. A schematic representation of methane adsorption in shale structure is described in Figure 7.  The grand conical ensemble ( ), in which there is the chemical potential of gas (μ), volume (V), and temperature (T), describes a system in open contact with a reservoir. The grand conical ensemble is close to the adsorption case by performing a series of insertions and deletions of adsorbate gas molecules. In the ensemble, the μ and the T are fixed, and the equilibrium status is that the adsorbate molecules have the same chemical potential and temperature in the reservoir. In the Metropolis scheme, a trial configuration is compared to the original configuration. The "sorption isotherm" module in Sorption of Material Studio with Grand Canonical Monte Carlo (GCMC) and a COMPASS forcefield was used to simulate the adsorption capacity of gases on the shale molecular cells [51]. Methane adsorption on shale models is simulated at 90 ℃ (363.15 K) and 120 ℃ (393.15 K) with pressure up to 30 MPa. The Anderson thermal bath was used to keep the temperature constant, while the cut-off radius was set as 12.5 Å. The configuration of adsorption isotherm measurements is described in Table 2. A schematic representation of methane adsorption in shale structure is described in Figure 7.
In the GCMC simulations, gases inside the slit pores are assumed to be in equilibrium with an external bulk reservoir under the same temperature and chemical potentials.   In the GCMC simulations, gases inside the slit pores are assumed to be in equilibrium with an external bulk reservoir under the same temperature and chemical potentials.

GCMC Adsorption Calculations
The fixed parameters in the GCMC ensemble were the volume, temperature, and chemical potential of the adsorbed molecules. The chemical potential is substituted by pressure, and the fugacity coefficient of an ideal gas reservoir and the fugacity of gas molecules can be calculated by the equation of state, such as the Peng-Robinson equation [52,53].
The correlation between the fugacity and pressure can be expressed as: where ( f ) is the fugacity in MPa, P is the pressure, R is the gas constant, T is the temperature in K, and a and b are EOS constants. The a and b can be expressed as follows: where T c is the critical temperature = 190.55 K and P c is the critical pressure = 4:59 MPa. α(T) can be expressed as: where T r is the reduced temperature, T r = T/T c , and k can be expressed as: ω is the eccentric factor for methane, which is 0.0113 [53]. The total loading capacity of adsorbate gas molecules can be obtained using Equation (6) [54].
where q is the absolute adsorption in mmol/g, m s is the adsorbent mass in the simulation box (g), n is the ensemble of adsorbing methane molecules in the simulation box, and N A is Avogadro's constant (6.0221 × 10 23 mol −1 ).

Model Validation
Physical density is an essential criterion for evaluating the reasonability of the molecular model. Therefore, the simulated densities of shale models are compared with documented experimental data from three shale formations. The comparison of simulation densities and experimental data is tabulated in Table 3. The simulated densities of shale structures in this work are 2.52 g/cm 3 and 2.44 g/cm 3 for TOC 2 wt.% and 5 wt.%, respectively, consistent with the density of actual shale from the Montney, Antrim, and Eagle Ford formations. The simulated densities of the shale models created in this work are consistent with the data in the literature, with the density decreasing as TOC content increases. A lower density value was observed in the 5 wt.% TOC shale model compared to 2 wt.% as a result of the theoretical basis that the kerogen density in shale is significantly lower (1.03-1.10 g/cm 3 ) compared to the density of inorganic minerals such as quartz (2.65 g/cm 3 ) and clay (2.0-2.9 g/cm 3 ) [56].

Comparison with Experimental Data
The ability of the new proposed shale molecular model's methane adsorption isotherms was tested by comparing the simulated adsorption results to the literature. The experimental values were gathered from Gasparik et al. (2014) [57]. For a clear comparison, the explicit boundary for adsorption isotherms measurement was conducted at 90 • C ± 10 • C and 120 • C ± 10 • C, and TOC contents between 2 wt.% and 5 wt. % were selected. The proximities of both techniques were graphically described in Figure 8.

Comparison with Experimental Data
The ability of the new proposed shale molecular model's methane adsorption isotherms was tested by comparing the simulated adsorption results to the literature. The experimental values were gathered from Gasparik et al. (2014) [57]. For a clear comparison, the explicit boundary for adsorption isotherms measurement was conducted at 90 C  10 C and 120 C  10 C, and TOC contents between 2 wt.% and 5 wt. % were selected. The proximities of both techniques were graphically described in Figure 8.  [57]. "Adapted with permission from Ref. [57]. 2023, Elsevier [License No. 5479231488060]".
The simulated adsorption isotherms measurements were conducted at 90 ℃ and 120 ℃ with pressure up to 30 MPa. In this study, the adsorption isotherms were calculated using GCMC simulations with a 5 nm pore size. The adsorption isotherm shows that the adsorption capacity increases as the TOC increases. This trend was observed in most adsorption isotherm measurements conducted in the literature [58][59][60].
Constructed shale models in this work also successfully tabulated the effect of temperature. The adsorption capacity in nanopores decreases as the temperature increases 23 . As adsorption is an exothermic process, the gas required higher energy to adsorb on the shale surface; thus, a lower adsorption capacity was observed. This relation was first reported for adsorption studies in Devonian shales by Lu et al. (1995) [61]. Their studies indicated that the sorption capacity decreases with temperature. This observation also concurs with other experimental values [62][63][64].
The applicability of simulated shale models to determine the adsorption capacity of shales was analysed via a comparison with the experimental data. The Average Relative Error (ARE) function was adopted to determine simulation values' deviations from experimental data (Equation (7)). Table 4 shows the ARE analysis values of the simulated adsorption with experimental data.  [57]. "Adapted with permission from Ref. [57]. 2023, Elsevier [License No. 5479231488060]".
The simulated adsorption isotherms measurements were conducted at 90°C and 120°C with pressure up to 30 MPa. In this study, the adsorption isotherms were calculated using GCMC simulations with a 5 nm pore size. The adsorption isotherm shows that the adsorption capacity increases as the TOC increases. This trend was observed in most adsorption isotherm measurements conducted in the literature [58][59][60].
Constructed shale models in this work also successfully tabulated the effect of temperature. The adsorption capacity in nanopores decreases as the temperature increases [23]. As adsorption is an exothermic process, the gas required higher energy to adsorb on the shale surface; thus, a lower adsorption capacity was observed. This relation was first reported for adsorption studies in Devonian shales by Lu et al. (1995) [61]. Their studies indicated that the sorption capacity decreases with temperature. This observation also concurs with other experimental values [62][63][64].
The applicability of simulated shale models to determine the adsorption capacity of shales was analysed via a comparison with the experimental data. The Average Relative Error (ARE) function was adopted to determine simulation values' deviations from ex-perimental data (Equation (7)). Table 4 shows the ARE analysis values of the simulated adsorption with experimental data.
q e,simulated − q e,experimental q e,experimental i (7) The Average Error Analysis (ARE) was used to assess the applicability of the proposed amorphous shale model to replicate the laboratory adsorption isotherm measurements of actual shale. The ARE function showed that the amorphous shale shows good agreement with experimental measurements of all Barnett shale samples with an average of 5.0% error and slightly higher for the Haynesville samples with 8.0% error. The lowest ARE value of the simulated with experimental values is shown in bold in Table 3.
Several studies have shown that kerogen type is another important parameter controlling the adsorption capacity in shale [65,66]. The methane adsorption capacities of kerogen increased according to the type of kerogen. The adsorption capacity is the highest in Type III kerogen, followed by Type II and Type I. The type of kerogen results from thermal maturity. Thus, the effect of different types of kerogen in the amorphous shale must be considered in constructing a molecular model for shale.
The differences in experimental and simulation values may be due to the amorphous packing of kerogen and quartz that it had to undergo to form the shale model. Moreover, actual shale formations have more minerals than presented in this work. Several studies on the formation lithologies of Haynesville and Barnett consist of black siliceous shale, limestone, and minor dolomite. Mineralogical characterisation using X-ray diffraction (XRD) analysis conducted by Sone and Zoback (2010) on Haynesville and Barnett samples observed that different minerals are found in the samples [67].A summary of mineral compositions from each sample is summarised in Table 5. The intrinsic mineralogical composition possessed by actual shale shows the importance of developing an improved realistic shale model with a structure attainable in reservoir conditions. In this work, we developed an approach to creating a reproducible shale model for adsorption analysis using molecular dynamic analysis.

Conclusions
This work uses molecular dynamics simulation with a new proposed intrinsic shale composition arrangement by varying kerogen cells following the TOC composition and incorporating quartz as inorganic components. Two amorphous shale structures composed of kerogen and quartz were constructed. The kerogen content was varied to replicate the shale with 2 wt.% and 5 wt.% Total Organic Carbon (TOC) content with a 5 nm pore size. The simulated densities of shale structures in this work are 2.52 g/cm 3 and 2.44 g/cm 3 for TOC 2 wt.% and 5 wt.%, respectively, consistent with the density of actual shale from the Montney, Antrim, and Eagle Ford formations. The Average Error Analysis (ARE) was used to assess the applicability of the proposed amorphous shale model to replicate the laboratory adsorption isotherm measurements of actual shale. The ARE function showed that the amorphous shale shows good agreement with experimental measurements of all Barnett shale samples with an average of 5.0% error and slightly higher for the Haynesville samples with 8.0% error. The differences between the experimental adsorption measurement and simulation resulted from the amorphous packing, and actual shales have more minerals than the simulated model.  Data Availability Statement: Data available on request due to privacy restrictions. The data presented in this study are available on request from the corresponding author. The data are not publicly available due to the part of the thesis dissertation own by corresponding author.

Conflicts of Interest:
The authors declare no conflict of interest.