Understanding Confined Fluids in Shale Gas Systems

Given the complexity of shale gas at high pressures, researchers aim to characterize the thermodynamic properties of confined fluids using a mixture of experimental, modeling, and simulation techniques. In this work we frequently use the predictive capabilities of simulation to couple the property results to models. The overall results are then compared to experimental data for verification purposes. We employ a Monte Carlo simulation technique to ensure that a simple linear mixing rule for internal energies of departure holds thereby allowing pure component data to extend to mixtures. The results are coupled to the Gibbs-Helmholtz Constrained equation of state allowing for bulk-scale bubble point reduction predictions. In addition, the sensitivity of the results is determined. Adsorption of n-alkanes at high pressure conditions are studied as a function of carbon chain length, temperature, and pore throat size (14.2 Å to 19.88 Å) to give an overall picture of shale gas behavior at reservoir conditions. A simple model is shown to provide a reasonable estimate of the isotherms at high pressures up to 500 bar and a temperature range of 300 K to 550 K. Under the assumption of ideal site-site interactions, mixtures are predicted for methane/ethane and methane/ethane/propane systems and compared to work in the literature. An important aspect of this work is the verification to experimental data; we expand on recent work by characterizing the experimental to simulation data in a robust manner. Quantitative agreement is achieved when estimating the surface area and void volume of the porous material.

capabilities of simulation to couple the property results to models. The overall results are then compared to experimental data for verification purposes.     Excess and Net adsorption isotherms compared to experimental data [15] . statements about what will occur in shale gas markets. Since markets are unpredictable, in the event of a continued oil price downturn throughout 2040, shale gas production will still remain relevant with an annual production of 40 billion cubic feet per day [4]. Depending on the amount extracted, this could lead to 20 -40 years of global energy, which could be used to offset the depleting conventional oil resources. Technological increases are leading to a reduction in drilling costs and an increased drilling efficiency in major reservoirs such as the Bakken, Marcellus, and Eagle Ford formations. For example, in China great progress has been achieved with the first commercial horizontal well to recover 16.7x10 4 m 3 /day after 15 stages of fracturing [5]. As a result, shale gas has the potential to remain a top energy resource for the near future as conventional supplies diminish. However, this success has been stymied by the complex nature of adsorbed shale gas not limited to size, shape, surface, chemical characteristics of individual pores. Finally, shale gas production could be limited by regulatory policies that limit of CO2 emissions.

Brief Adsorption Literature Background
In order to develop a theoretical model for the amount of gas adsorbed into the shale (adsorption) we must first look at available experimental investigations and theoretical models currently available in the literature (see Fig. 1 [7]. At higher pressures than experimentally available, the results have been difficult to obtain due to the complexity of the interactions (pore throat size, chemical composition, fluid-fluid interactions, etc.) of adsorbed gases inside shale rock [8]. To circumvent this, models have been developed to account for these interactions at high pressures. These models are not limited by varying pore throat sizes and shape distribution in the nanometer regime, chemical composition, surface area profile, etc.
but they typically rely on empirical methods for adsorption modeling.
Since the adsorption literature is vast, we only highlight methods relevant to this thesis.

1) Traditional Adsorption Models
Traditional adsorption models are used to describe shale gas in the low-pressure regime by fitting the model to available data (less than 1 bar  [9]. While these models provide excellent agreement with experimental data in the low-pressure regime, the authors note that they typically do not predict high-pressure adsorption due to large interactive forces between shale gas molecules, pore-filling geometry, and complex gas-rock interactions. In addition, they are not suited for extrapolation at high pressures if they are not properly fitted or no data is present. These models can be extended to describe mixtures with ideal adsorption solution theory (IAST) [10]. As with all models, there are assumptions and limitations, the adsorbed molecule in IAST is assumed to be well mixed and each molecule has access to the same surface area.

2) Equation of states (EOS)
There are numerous equations of state developed for the description of adsorption phenomena applied to shale gas reservoirs. These equations describe the pressure, temperature, volume, and compositions relationship of fluid mixtures and can be applied at high pressures. For example, a modified Peng Robinson equation was used to describe meso-porous materials MCM-41 and 13X for methane and carbon dioxide up to 1.2 bar by using a 16-constant expression fitted to simulated bulk fluid densities in thermodynamic equilibrium with confined media [11].
Again, most of these equations are regressed to experimental data in order to describe low-pressure adsorption (less than 100 bar). Myers [12] used a rigorous approach involving desorption functions and provides an example for ethylene adsorbed on NaX zeolite up to 1.4 bar. Earlier, more empirical work by Myers was based on the use of solution theory by treating the adsorbent as the solvent and small molecules at low concentrations as the solute [13] for adsorption pressures up to 1000 bar. The Elliott-Suresh-Donohue (ESD) equation has been applied to gas adsorption on activated carbon for a number of components including acetylene, propylene, and ethylene for pressures up to 2 bar by incorporating a simplified local density model gas adsorption [14]. Finally, the Bender EOS also uses empirical relationships to describe adsorption up to 500 bar for nitrogen and methane [15].

3) Molecular simulation
In this thesis, adsorption is studied using molecular simulation [16] in order to estimate thermodynamic properties where experimental data are not available and/or not easily obtained. Simulations have been used to model diffusion processes in porous materials (e.g. water diffusing through graphene-based nanopores) and aid in coarse grained reservoir simulations, which require transport phenomena parameters [17]. For example, recent work has been used to describe the selectivity (i.e., the extent of preferred adsorption onto a framework) of CO2/CH4 onto an organic-rich shale framework [18]. Similar recently published work shows the interplay between CO2 and CH4 adsorption observed in experimental samples for Barnett 26-Ha, Haynesville GU1-2, and Haynesville TWG3-1 rock [19].
To be fair, traditional numerical adsorption techniques in the literature were developed for zeolite and metal-organic framework applications in the low-pressure regime used for describing gas dehydration, small molecule separation, oxygen generation, etc [20]. However, with the recent interest in shale gas, high-pressures applications of the traditional model geometries of adsorption (e.g. slit pore) have led to a range of results from over-predictions to negative adsorption [21]. In addition, the use of inaccurate density calculations by equations of state have led to overcorrections when attempting to link simulation to experimental results [22]. Another objective of this thesis was to correct these limitations in the literature by using the GHC equation of state framework.

Advancement of Knowledge
Yet another goal of this thesis is to provide the research community at large with shale gas adsorption data at high pressures over a wide range of temperatures and pore throat diameters. These simulation efforts will focus on areas where experimental data is not available. The main contribution of this thesis is the study of internal energies of departure at the nano-scale and the subsequent up-scaling of this molecular information to the bulk scale using the GHC EOS. This work will, in the author's opinion, be quite useful and focuses specifically on molecular simulation to describe This mixing rule has been shown to save computational time and easily allow for property and phase equilibrium computations for mixtures. For example, for a mixture of water and methane, molecular simulation would be used to generate information for pure water and for methane, and then the mixture rule would be applied to estimate water/methane mixture at any composition.
Therefore, the components of this research are:

1)
Validate the linear mixing rule for confined fluids. Pure component internal energy of departure data for methane, carbon dioxide and n-alkanes and mixtures of methane/n-alkane & CO2/n-alkane in confined spaces are gathered using molecular simulation. The linear mixing rule will be applied and compared to results for mixtures. Corresponding percent errors as well as internal energies of departure sensitivities are reported.

2)
In the context of adsorption, molecular information is difficult to compute at high pressures. In this work, isothermal adsorption behavior at high-pressure is simulated for mixtures over a wide range of pressures and temperatures. This work investigates the impact of confinement with varying pore throat diameters and hydrocarbon chain length leading to a discussion on their impact on fluid internal energies of departure.

3)
Where available, recent adsorption methods [22] are applied and compared to experimental data. This method uses a probing molecule to estimate the accessible pore volume and surface area in order to normalize simulation data for comparison to experimental data. This aspect of the work provides much needed insight for the thermodynamic community.

EQUATION OF STATE COMPUTATIONS FOR CONFINED FLUIDS
The following manuscript is published as a special issue in Computers and Chemical Engineering.

Literature Survey
Interest in physical properties and phase behavior of shale gas and LTO is relatively recent and the literature on the subject is somewhat sparse. Early studies in reservoir and petroleum engineering from the 1940's to 2000 [5] suggested that capillary effects on phase behavior were negligible. However, all recent studies, which are largely numerical in nature, include interfacial tension between immiscible phases as part of the model. See [6]- [9].
The current approach to modeling fluid properties and phase equilibrium in tight pores in reservoir simulation consists of 1) An equation of state, e.g., [10]- [12] 2) A correlation (e.g., the parachor equation [8] or MacLeod-Sugden correlation [13] to calculate interfacial tension, s. 3) An estimate of capillary pressure, = 2 (e.g., using the Young-Laplace equation or Leverett J functions) [14].  [3] for the details used in computing internal energies of departure in the NPT ensemble using the MCCCS Towhee software system, version 7.10 [15].
However, in this work, the more recent RASPA software was used for all Monte Carlo simulations [1]. Therefore, the first issue to be resolved is to show that the same statistical results for unconfined NTP Monte Carlo simulations can be obtained for mixtures studied by Kelly and Lucia [3] using RASPA.  For these unconfined NPT simulations, volume, translation, and rotation move frequencies were set to 2.48%, 48.78%, and 48.78% respectively while radial cutoff distances were adjusted to include all interactions in the system. For electrostatic forces, Coulomb interactions were calculated using Ewald summations as defined in Dubbeldam (p. 81, [16]).
For larger alkane molecules, Configurations Bias Monte Carlo (CBMC) summations are needed and the corresponding volume, translation, rotation, and CBMC frequencies were set to 1.23%, 24.69%, 24.69%, and 24.69% respectively. Larger molecules also have torsion and have an ideal gas contribution to the internal energy of departure. For this, additional simulations were performed using a single molecule 20 in the canonical (NVT) ensemble with CMBC moves. The pressure in the NPT ensemble was measured by computing the negative of the stress tensor as defined in Dubbeldam et al. ([1], [17]) For mixtures, an identity switch move was used with a frequency of 5% to ensure higher mixing probability. Finally, all unconfined simulations were also performed without analytical tail cutoff corrections, and at the same temperature and pressure for the NPT ensemble, for the purpose of comparing confined and unconfined simulation results.

The Confined Canonical Ensemble.
Confined canonical (NVT) ensemble Monte Carlo simulations in this study were performed using the screening study of Dubbeldam [18] and the more recent energy Framework information can be found in supporting information (S178) of Dubbeldam et al. [18]. Also, while the energies of the graphite walls and inter-and intramolecular energies of the molecules in the nano-channel were computed, the and results reported for the confined cases only include the energies of the molecules inside the graphite square channels.  There may be several reasons for this. First, mixing within nano-channels can be problematic due to restricted particle movement and the presence of the nano-channel walls, especially for larger molecules.
The sensitivity of molar volume to changes in must be computed by directly solving the equation of state for its volume or density roots. In Kelly and Lucia [3], the corresponding relative sensitivities of and molar volume (or density) to 5% uncertainty in were less than 4% and 1.5% respectively for unconfined fluids with tail corrections.   Note that the maximum sensitivity of to 5% uncertainty in does not exceed 3.2% for both the unconfined and confined cases. With the exception of the waterhexane mixture, the sensitivity in in the unconfined and confined cases is less than 0.5%. For water/n-hexane in unconfined and confined space the sensitivity of is approximately 2.5% and 3.2% respectively. In addition, the mixture water-hexane exhibits the greatest difference in relative sensitivity (~ 0.4%) while all relative sensitivities of to a 5% uncertainty in for both the confined and unconfined mixtures is less than 1%.

Bubble Point
It is straightforward to estimate the sensitivity of bubble point pressure to changes in and illustrate that confinement results in a decrease in bubble point pressure.
Let be any bubble point pressure at fixed and . The pressure expression for the GHC equation is We assume that is insensitive to for liquids since previous numerical experiments clearly show it is less than 1%. See also Lucia et al. (Fig. 7, p. 85, [2]).
Therefore the partial derivative of pressure with respect to is given by Using the chain rule we have that Note that for the existence of a bubble point requires that < and therefore the Therefore the change in bubble point pressure due to confinement is given by where ∆ represents the change in internal energy of departure of a fluid mixture in confined space minus that in unconfined space at the same x and T. That is, ∆ in Eq. 2.8 is defined as Where , and , are the confined and unconfined pure component internal energies of departure respectively. Substituting Eq. 2.9 into Eq. 2.8 gives.
Equation 2.10 clearly shows that bubble point reduction is a complex function of composition, temperature, and pore radius because all pure component , are functions of temperature and pore radius.
All of the confined data presented in this article thus far corresponds to a pore radius of 0.25 nm which, while valid, is quite small. However, in order to provide bubble point reduction estimates that are more representative of pore sizes encountered in practice, confined NVT simulations were run for a number of pore radii. As pore radius increases, confined mixture internal energies of departure should increase and approach the unconfined in the limit. Table 2.8 shows this effect of pore radius on the reduction in bubble point pressure for a 50-50 mol% mixture of methane and octane at 300 K.

HIGH PRESSURES
The following manuscript is submitted to the Journal of Petroleum Science and Engineering. Finally, Ideal Adsorbed Solution Theory (IAST) is applied to the resulting Langmuir isotherms in order to provide some insight into the phase equilibria for applicable shale gas mixtures.

Introduction
Fluids in confinement generally have different properties than traditional bulk fluid properties [1]. These unique properties include (1) highly structured geometry [2], (2) decreased mobility in confined directions, which strongly effects sampling pressures [3], and (3) order-disorder transformations in slit-like pores [4]. Applying computational modeling to characterize n-alkanes in confinement has the potential to build fundamental understanding of Light Tight Oil (LTO) and Shale system [5] properties and phase behavior.

Literature Survey
The recent boom in shale gas technology has led to an increase in global production efforts [6]. In 2013 the United States Energy Information Administration (EIA) estimated the total unproven technically recoverable shale gas reserves across the globe [7] and clearly show that shale gas will play a major role in the world energy approximately 117.2 bar [9]. Additional data can be found for methane, nitrogen, and carbon dioxide on Woodford shale from the Payne, Hancok, and Caney county formations up to 125 bar [10]. At higher pressures, the results have been difficult to obtain due to the complexity of the interactions of adsorbed gases inside confined pores such as nanoporous (< 50 nm) shale rock [3]. To circumvent this, models have been developed to account for these interactions at high pressures. These models are not limited by varying pore throat sizes and shape distribution in the nanometer regime, chemical composition, surface area profile, etc. and typically rely on empirical methods for adsorption modeling. In this work, we propose a methodology for predicting adsorption isotherms of shale gas molecules at high pressures using a slit pore model. Finally, since the adsorption literature is vast, we only highlight methods relevant to our study.

Adsorption Models:
Recent work has reviewed the Henry, Freundlich, Langmuir, Dubinin-Radushkevich, Radke-Prausnitz, Toth, Langmuir-Freundlich models applied to shale gas adsorption for methane and carbon dioxide up to 140 bar [11]. While these models provide excellent agreement with experimental data in the low-pressure regime, the authors note that they typically do not predict high-pressure adsorption due to large interactive forces between adsorbates, pore-filling geometry, and complex adsorbate-wall interactions. Pure component adsorption results have been extended to multicomponent systems [e.g. Ideal Adsorption Solution Theory (IAST)] using fitting models [12], [13]. Well-fitted pure component models are essential for predicting multi-component adsorption.

Equation of states (EOS):
There are numerous equations of state derived for the description of adsorption onto a carbon framework that may be used for shale gases. However, most of these equations are regressed to experimental data in order to describe low-pressure adsorption (less than 100 bar). Myers [14] used a rigorous approach involving desorption functions to derived explicit expressions for G(T,P), in which the author gives an example for ethylene adsorbed on NaX zeolite up to 1.4 bar. Earlier, more empirical work by Myers was based on the use of solution theory by treating the adsorbent as the solvent and small molecules at low concentrations as the solute [15] for adsorption pressures up to 1000 bar. The Elliott-Suresh-Donohue (ESD) equation has been applied to gas adsorption on activated carbon for a number of components including acetylene, propylene, and ethylene for pressures up to 2 bar by incorporating a simplified local density model gas adsorption [16]. A modified Peng-Robinson EOS was used to describe meso-porous materials MCM-41 and 13X for methane and carbon dioxide up to 1.2 bar by using a 16-constant expression [17] fitted to simulated bulk fluid densities in thermodynamic equilibrium with confined media. Finally, the Bender EOS also uses empirical relationships to describe adsorption. In this work, high-pressure adsorption up to 500 bar experimental data is described for nitrogen and methane [18].

Molecular simulation:
Adsorption can be modeled by molecular simulation (Monte Carlo or molecular dynamics) and used to estimate thermodynamic and transport properties where experimental data is not available. Molecular dynamic simulations have been used to obtain diffusion processes in porous materials (e.g. water diffusing through graphenebased nanopores) and aid in coarse grain reservoir simulators which require transport phenomena parameters [19]. Adsorption processes can also be described using Monte Carlo simulations in the grand canonical ensemble. For example, recent work has been used to describe the selectivity (i.e., the extent of preferred adsorption onto the framework) of CO2/CH4 onto an organic-rich shale framework [20]. Similar recently published work shows the interplay between CO2 and CH4 adsorption observed in experimental samples for Barnett 26-Ha, Haynesville GU1-2, and Haynesville TWG3-1 rock [21].
There are clearly a limited number of adsorption models applicable to shale gas at high pressures over a wide range of temperatures and pore throat diameters. One novel aspect of this work is the development of a framework based on the Canonical (NVT) Ensemble to compute gas adsorption in regions of high pressures for n-alkanes ranging from methane to hexadecane. The primary goal of this study is to provide a library of adsorption and internal energy of departure data for n-alkanes typically found in the shale gases and light tight oils. In our opinion, this data would be very useful for all research communities that employ multi-adsorption models such as revised multi-component Langmuir models, IAST, and EOS. Furthermore, other recent work from our group has used results for internal energies of departure in confined spaces and a simple EOS linear mixing rule for mixture internal energies of departure to predict bubble point reduction using the Gibbs-Helmholtz constrained (GHC) equation [22]. The proposed work is quite useful for reservoir simulators that desire a predictive and multi-scale thermodynamic approach for the description of shale gas phenomena.
In this work the choice of ensemble is the Canonical approach. It widely known that high-pressure adsorption in the Grand Canonical Monte Carlo (GCMC) ensemble is challenging [23]- [25] and often leads to inaccurate estimates of pore pressure.
Difficulties arise when the system becomes denser and, as a result, the acceptance of insertion moves diminishes, which in turn effects convergence. Low acceptance rates can be computational challenging to resolve, requiring grid interpolation techniques in the GCMC ensemble [12]. In contrast, in the NVT ensemble the pressure can be estimated using the negative of the stress tensor and general fluctuation expressions [26], both of which are implemented in the RASPA software system [27] for both Monte Carlo and Molecular Dynamics applications. However, it is important to note that computing pressure in this manner often requires that the adsorbed fluid be in full contact with the framework [28], which essentially means that there cannot be any (GHC) EOS [22]. As a result, a subsequent goal of this work is to determine the impact of adsorption on computed U D in confinement within the GHC EOS framework [22], [29].

Computational Methodology
In this work the nano-channel framework shown in Fig. 3  The frequencies translation and rotation moves without CBMC moves were set at 50% each. For longer n-alkane chains, starting with butane, the corresponding frequencies for translation, rotation, and CBMC moves were 33.33% each. Reported averages are the result of using a 5-block average during the production phase of the simulations.
Simulations were run on three custom built computers with AMD 1090T 3.2 GHz and AMD FX8300 processors in double precision arithmetic using the GNU compiler. An example snapshot of the simulation is shown in Fig. 3.1 for n-hexadecane. Additional nano-channel framework specifications can be found in the supporting information of Dubbeldam et al. [30], in the nano-channel framework section, and have been reproduced here for convenience. Also, while the energies between the graphite walls and adsorbates were included, U D results only include energies between adsorbates.  [30] in nano-channel specifications An overview of the steps needed to compute adsorption isotherms is given in Fig. 3 predict adsorption for ternary mixtures [13]. The computational procedure is repeated in confined space for a range of pore throat diameters for each n-alkane. Benchmark simulations were used to establish a comparison between the methodology proposed in Fig. 3.2 and the procedure in Dubbeldam et al. [30]. Results of this comparison for U D are given in Table 3.2 for the NVT and GCMC ensembles. The corresponding pressures for each ensemble can be found in columns 2 and 3 with the standard deviation for the NVT ensemble computed pressures shown in parentheses. It is important to note that the same graphite nano-channel framework (i.e., in supporting information of Dubbeldam et al. [30]) was used for both ensembles. Also, the computation of U D is a required component for determining ∆H i ,ads , which is described in Dubbeldam et al. [30].
From Table 3.2, it is clear that the results are in very good agreement with a percent error ranging from 1.15 to 3.51%. With the exception of 25 bar, the standard deviations are higher for Dubbeldam et al. [30] than in this work.

Langmuir Equation & Ideal Adsorbed Solution Theory
In this portion of the work, the python package, pyIAST, was used to perform Ideal Adsorbed Solution Theory (IAST) calculations for applicable shale gas mixtures [13].
As with all models, there are assumptions and limitations. In this case, a ridged uniform framework was used in the Canonical ensemble, which is valid for IAST.
However, a rigid framework implies that thermodynamic properties are independent of volume during sampling because volume is constant [13]. Finally, shale gas in the pore is assumed to be well mixed and each adsorbate has access to the same surface area. Given these assumptions, the following well-known Langmuir equation was

Simulation Results
Adsorption results for methane, ethane, propane, butane, octane, and hexadecane over wide ranges of temperature, pressure, and pore throat diameter are given in Figs. 3.3-3.14. Results are reported as adsorption isotherms and values of U D .
Moreover, the extent of ordering generally depends on the specific adsorbate and adsorption material under consideration.

Case 1:
Comparison of simulated methane adsorption as a function of pore throat to experimental data from Heller and Zoback [9].
Experimental data used for comparison is taken from Heller and Zoback [9] who report methane and carbon dioxide adsorption on shale gas samples taken from the Barnett [14]. Note that the experimental results compared in this work are for methane at relatively low pressure and extend up to approximately 120 bars while simulation results go up to 450 bar. Also note the proposed simulation methodology provides results that are in qualitative agreement with the experimental data. Again, it is important to stress that our simulation methodology is predictive and not fitted to any adsorption data whatsoever.
The results in Fig. 3.18 clearly demonstrate that a simple slit pore model may be used to predict shale gas adsorption.  In that work, the authors obtained the data by fixing the temperature at 298.15K and varying pressure up to 60 bar. The data in Table 3 demonstrated that a slit pore model has the potential to predict shale gas adsorption, the next step is to fit Langmuir isotherms to all adsorption data in this work (found in Appendix 6.2). IAST calculations can then extend our simulation work to mixtures by using the Langmuir parameters as an input.
IAST calculations were generated to describe the Marcellus and Barnett Shale data in Table 3.

2015)
Finally, we used a rigid framework for the nanoporous material (see section 3.4). This makes it possible to use IAST calculations to study the behavior of any multicomponent mixtures using the pure component n-alkane Langmuir parameters given in Appendix 6.2, provided we invoke the assumptions given in Simon et al. [13].
Sharma et al. [33] assumed that the 4-5 mol% propane given in Table 3  3.19 A clearly shows that there is a relatively large amount of propane uptake of approximately 1.6 mmol/g for the given methane/ethane/propane mixture at 313K and 200 bar. What is important to emphasize here is that our adsorption isotherms can be used to provide quick mixture estimations in cases where the composition of propane cannot be neglected. On a broader note, the use of IAST coupled with the Langmuir parameters given in Appendix 6.2 provides a means of quantifying mixing effects of n-alkanes up to hexadecane, which in turn, allows for estimations and references of adsorption conditions.

Conclusion
In this work adsorption isotherms and UD of n-alkanes in a graphite nano-channel were studied over a range of pressures in the high adsorption regime. Results clearly showed that differences in n-alkane adsorption isotherms decrease overall as the carbon chain length increases at any given pore throat diameter. This behavior was attributed to occupancy limitations for larger n-alkanes. In addition, adsorption data was fit to Langmuir isotherms and corresponding parameters were determined. Pure component simulation results exhibit similar trends to experimental data (Heller and Zoback [9]) for methane adsorption on shale gas samples taken from the Barnett

Abstract
In order to compare absolute adsorption simulation data to experimental excess and net adsorption data for various adsorption processes, a conversion technique that uses the Gibbs dividing surface derivation to define the upper and lower limits of adsorption phenomena is used [1]. In the context of shale gas adsorption over prediction, a more recent conversion technique has been suggested to normalize both excess and net adsorption data by using the framework surface area to volume ratio [2], which when combined with the Gibbs dividing surface has been shown to be effective when linking simulation to experimental data [2]. The surface area to volume ratio framework has only been employed in the Grand Canonical (GCMC) ensemble. Finally, the conversion adsorption results are compared to existing methods and experimental data with exceptional agreement when compared to traditional methods.

Introduction
Production of non-conventional hydrocarbon resources is expected to increase throughout 2050 in order to compensate for diminishing conventional reservoir supplies [3]. This has led to an increase in interest in non-conventional reservoirs such as shale gas and light tight oil (LTO) worldwide. Coupled with this effort is the reduction in costs due to technological advances like fracking for major reservoir locations such as the Bakken, Marcellus, and Eagle Ford reservoirs.
The United States remains the world's top producer of natural gas and is in a position to continue to grow the shale gas market [3]. Outside of the United States, production of shale gas in China and Canada are expected to grow from 0.5 to 22 billion and 5 to 8 billion cubic feet per day respectively [3]. It is important to recognize that the growth of shale gas will depend on the market conditions for natural gas that, in turn, are dependent on many other economic factors such as the price of oil. Short term forecasting by the U.S. Energy Information Administration (EIA) predicts oil prices will recover around 2019. However, shale and LTO projections are predicted to increase 1.3 million barrels per day. In the event of a continued oil price downturn, shale gas production will still increase by 35 billion cubic feet per day from 2015 to 2017 [4]. In the long term, shale gas will play a major role in the world energy portfolio and is predicted to increase 70 billion cubic feet per day from 2015 to 2040 [4]. This it is clear that regardless of market conditions, Shale gas will remain a top energy resource in the near future as conventional supplies diminish.
To maximize the efficiency of shale gas production, a comprehensive understanding of the multi-scale problem ranging from the nano-scale to bulk scale is needed. In this work, the focus is on understanding the impact of the nano-scale since it has been demonstrated that the potential maximum efficiency has been hindered by the complex nature of fluid-pore interactions involving pore size, shape, chemical, and surface area distribution of individual pore frameworks [5].
In the nano-scale shale media (< 50 nm), the pore is occupied by a heterogeneous mixture of hydrocarbons and usually the main component is methane ranging in composition from 60 to 80% [6]. The choice of the GCMC is convenient because the chemical potential is constant and adsorption is studied over a range of chemical potentials requiring multiple simulations. The bulk pressure from a GCMC ensemble can be obtained by relating chemical potential to an equation of state such as the Peng-Robinson equation [9].
However, adsorption at high pressure conditions in the GCMC ensemble is challenging when the fluid densities is high ( [10], [11]. If an EOS is used to relate chemical potential to bulk pressure, there can also be inaccurate estimates of the bulk pressures and densities due to the accuracy of the EOS. Finally, low acceptance rates are reported at high pressures in the GCMC ensemble due to limited particle transfer moves in dense regions ( [7]). Therefore, following recent work of Thomas and Lucia [12] the canonical ensemble will be used in this work since high-pressure adsorption does not rely on particle transfer moves. The overall approach is described in Thomas and Lucia [12]; here only a brief overview is given in section 6.1.
Once adsorption is obtained through a given Monte Carlo ensemble, the results must be compared to experimental data. This work focuses on relating NVT simulation results to experimental data. In the next sub-section three different methods used for converting absolute adsorption results to excess, net, and surface area adsorption are discussed and compared. Subsequently, these conversion methods are combined and applied to a canonical (NVT) system and then compared to experimental data for a high-pressure shale gas system.
The determination of for high pressure shale gas adsorption must be carefully considered when using an equation of state [13] since recent work has clearly shown that the Peng-Robinson (PR) equation of state (EOS) over predicts at high pressures. This, in turn, can lead to negative nex (see Fig. 5 of Chen et al. [2].
Because nex is sensitive to values of NPT simulations are used to determine gas density at high pressure. In the case of methane there have been many studies that use the PR EOS to convert absolute to excess adsorptions but none of these studies discusses the impact of the accuracy of PR at high pressures [9] on adsorption.
Additional simulations are required to calculate V (cm 3 /mol), which is typically referred to as the helium void (or dead space) volume as described in Myers and Monson [14]. In this context, helium simulation acts as the reference state for V when determining the pore space that is not occupied by the framework. Gumma and Talu [1] demonstrate that this reference volume can be calculated in simulation by use of a configurational integral for helium where ɸ is the collision diameter describing the helium-solid interactions, is the mass of the framework, and Boltzmann constant k.
There remains debate in the adsorption community over the correct values of ɸ for helium as any perturbation in ɸ will result in a different value of V thereby affecting the . In fact, some recent work has suggested that the adsorbing molecule be used as the reference molecule [2]. This approach showed a small change in accessible pore volume found in Table 1 of [2] for a Na-Montmorillonite (Na-MMT) simulation cell probed by methane and helium yielding a volume of 40.80 nm 3 and 40.94 nm 3 respectively. Nonetheless, even this small difference in accessible pore volume can lead to noticeable differences in excess adsorption as shown in Fig. 4 of Chen et al. [2]. However, this difference can be attributed to geometric interactions between the framework and the probing molecules because helium is a smaller molecule, than methane and tends to occupy spaces between the surface molecules [2]. An illustrative diagram of the geometric considerations of probing molecules is given in Fig. 6 of Chen et al. [2].
In section 4.5.4, the method of using the adsorbing molecule as the reference molecule for a graphite-nanochannel system is explored further and the effects of helium and methane probing molecules for the determination of V are studied. More specifically, for shale gas applications of interest in this work, methane is used as the probing molecule because methane only occupies space on top of the adsorbent surface and not interstitially, thereby occupying less space.

Net Adsorption:
Although the adsorption community typically uses Eq the reader is referred to Gumma and Talu [1].
is better suited for framework materials that have an explicit volume such as slit pores and is not suitable for flat surfaces without an explicitly defined volume.
Previous applications of the net adsorption approach involve more complex structures such as the metal-organic framework HKUST-1 [1]. As mentioned earlier, great care must be taken for the determination of with an EOS.

Surface Area Comparison:
Key findings by Chen et al. [2] attempt to link results of GCMC simulations to experimental data by converting absolute to excess adsorption isotherms. In this work, the authors showed that while bulk fluid densities (ρb) and accessible pore volume (Vfree) are important variables that influence adsorption conversion, the specific surface area (SSA) of the framework is key and demonstrated this for a Na-Montmorillonite system and qualitatively matched experimental results. However, one potential drawback of this method is that it is strongly dependent on the availability of experimental framework SSA data. Furthermore, the SSA approach is also heavily dependent on the type and location of the reservoir. For example, Ji et al. [5] performed an analysis with a Beckman Coulter SA3100 SSA analyzer from a number of shale gas reservoirs. They reported SSAs for monmorillonite (76.4 m 2 /g), I-S mixed layer clay (30.8 m 2 /g), kaolinite (15.3 m 2 /g), chlorite (11.7 m 2 /g), and illite (7.1 m 2 /g).
Fan et al. [15] studied the adsorption of a highly matured sample in Longmaxi (China) of shale gas system and cited a SSA of 15.10 m 2 /g. In this work the surface area of Silurian Longmaxi Formation shales is used, which have areas that range from 17.83 to 29.49 m 2 /g [16].

Overview of proposed work
In this work, NVT simulations are performed for multiple isotherms corresponding to existing experimental data in the open literature [15]. bath temperature in the sample holder given in Fan et al. [15].
The objective of the proposed numerical procedure is to compare simulation results to the experimental results reported by Chen et al. [2]. The proposed computational methodology addresses the current limitations in the literature that lead to large differences between computationally determined and experimental adsorption profiles.
It is important to emphasize that this limitation has been resolved for the grand canonical ensemble but not for the canonical ensemble [2].  The framework used for the NVT simulations is shown in Fig. 4.2 and is a graphite nanochannel used to represent the adsorbent material consisting of 3776 atoms.
Specifications for the absorbent can be found in Table 4.1 [18]. The United Atom TraPPE force field was used for methane. The Lennard-Jones United Atom TraPPE force field was used to model methane [19]. Helium reference state simulations Lennard-Jones parameters from Bolboli Nojini et al. [20]. Parameters for the adsorbates are given in Table 4.2. Framework parameters can be found in Dubbeldam et al. [18]. Tail cut off corrections are generally used to estimate molecular interactions at very large distances with no walls; however they are not applicable for confinement.
Simulations were started with two unit cells: one containing an empty framework and the other containing a number of methane molecules on the outside of the slit pore on the outside of the pore throat. Since the United Atom TraPPE force field was used for methane, there was no need for rotation moves. Instead, a short Monte Carlo preequilibration step was run for 10,000 steps with 100 % translation frequency and a Molecular Dynamics step every 10 steps was employed to promote equilibration (see [21]). Only Monte Carlo translation moves were performed for 400,000 equilibration and production Monte Carlo cycles. The reported number of particles (N) was determined using a five block average thereby giving a 95% confidence interval. Once equilibrium (and adsorption) was achieved between the outer methane particles and those in the slit pore space, the bulk pressure was determined by using the standard virial pressure applied only to the bulk methane particles outside the framework [22], [23]. Thermodynamic equilibrium in the NVT ensemble was satisfied by using periodic boundary conditions in the directions of the pore throat and by placing additional atoms on the outside of the framework to provide equilibrium between the particle in the bulk and in the confined space. Since flexibility is not a strict requirement for thermodynamic equilibrium [24], a rigid framework was used to simulate small nanopores representing conditions in a shale gas reservoir.

Isothermal-Isobaric Ensemble (NPT)
In order to compute the bulk density of methane, NPT simulations were performed with 500 methane particles and runs were set to 400,000 equilibration and production Monte Carlo cycles. A radial cut off distance of 12 Å was used with tail cut off corrections applied. The frequencies for translation and volume moves were each set to 50%. Further simulation details can be found in Allen and Tildesley [23] and previously applied for different systems (e.g. water, hexane, CO2) [25]. Forcefield information can be found in Table 4.2.

Density Computations
Since the conversion of absolute to excess adsorption in Eq

Free Volume and Surface Area Computations
The free volume of the graphite nano-channel was determined using Eq. 4.2 which requires only particle insertion moves with the desired probing molecule to estimate the second virial coefficient [1]. Since the free volume is dependent on the geometry of the pore structure, it is not a function of pressure. However, the procedure must be performed at experimental temperature conditions to maintain a consistent reference point. Free volumes for methane and helium at a reference temperature of 298 K were determined. The void fraction of the slit pore (i.e., the empty space divided by the total volume) was estimated by particle insertions [21].
Surface area was computed using an auxiliary method provided by Dubbeldam [21].
The surface area computation consisted of rolling a probing molecules (e.g. nitrogen, helium, argon, etc) over the desired framework. Each framework atom location was assigned atom points that generate a sphere around individual framework atoms where the amount of overlap is computed. The probing atom was then rolled onto the surface of the framework atoms and the corresponding overlap was computed for the frameworkprobing molecule interactions. Finally, the fraction of overlap was multiplied by the area of the sphere resulting in the geometric surface area. Since the focus of this paper is on current methods for comparing simulation adsorption results to experimental data, the reader is referred to Connolly [26] for a rigorous description of this method for determining surface area.

Results and Discussion
In this section, the presented results demonstrate the estimations of bulk methane densities in 4.6. Use of a bulk fluid reference is described in section 4.  data [17] In Fig. 4 We have demonstrated the ability to indirectly compute the bulk fluid densities using a Monte Carlo and EOS approach as seen in Fig. 4 [27]. Now as with all techniques, there are advantages and disadvantages. The advantage is that the density can be directly estimated which saves computational time, this approach can be quite useful as geometric complexity of the framework is increased. However, the drawback is that the virial pressure that corresponds to the bulk pressure has error associated with it.
Fortunately, our previous work shows that if a multi-scale approach is desired, we can use the Gibbs-Helmholtz constrained EOS in this case. Specifically, if bubble point reduction estimations are desired, a sensitivity to a 5% uncertainty of confined fluid molar volume is less than 1% [12]. Since this drawback has been previous investigated and quantified, we focus on the adsorption aspect of confined fluids. In this work, we focus on the direct and indirect methods to compute bulk densities in the NVT ensemble and compare their impact on excess and net adsorption to experimental data (e.g. PR EOS, NPT ensemble, etc).

Comparison of excess and net adsorption to experimental data using conventional approach with NPT computed densities
Since simulation work employs smaller pore wall thickness than experimental conditions there will most likely be a larger specific surface area up to a few orders of magnitude at high pressures ( [2], [28]). The same key concept holds true for this work and provides an explanation for the large differences between the simulation and experimental conditions found in Fig. 4.5. The excess and net NPT adsorption curves in Fig. 4.5 are the benchmark for comparison purposes due to the accuracy of the NPT Monte Carlo methane densities as observed in Fig 4.3. Also, it should be noted that the net adsorption should be markedly lower than excess adsorption because the entire volume of the system is subtracted from the absolution adsorption contribution (see Eq. 4.3). On the contrast, the excess adsorption is typically higher than net adsorption at very high pressures in the supercritical region for light gases [14]. This is due to the lower amount of slit pore volume that is considered for excess adsorption (see Eq. 4.1). The reference pore volumes for the net and excess conversions can be found in Table 4.3. Figure 4.5: Excess and Net adsorption isotherms compared to experimental data [15] There is also debate in the literature concerning the negative adsorption for net adsorption conversion. As made clear in the literature the net adsorption terminology in the adsorption conversion literature lead to considerable confusion on the subject [1]. They further state that several thermodynamic properties such as departure functions are negative and still considered rigorous. Therefore, the assertion that the net adsorption should never be negative is incorrect [1]. An example of a negative net adsorption can be seen in Fig. 5 of Gumma and Talu [1] which compares absolute, excess, and net adsorption for methane on Norit R1 Extra at 298 K.

Surface area approach for linking excess and net adsorption curves
The intent here is to consider the surface area of both the simulation unit cell and experimental sample. Coupling the surface area information provides the thermodynamic community a link between the two methods in the NVT ensemble.  Fig. 4.5 that only consider the pore volume. Since there is a non-uniform SSA distribution for field conditions, the agreement with simulation results can be found in the upper and lower SSA range for the 363.15 K adsorption curves [15].
From inspection, there is better agreement with NPT excess adsorption curves with experimental data [15]. On the contrary, the PR EOS over predicts bulk methane densities at high pressures a further decrease in adsorption as expected. Furthermore, the excess adsorption curves capture a reasonable trend beyond the experimental pressure range.  The net adsorption NPT results in Fig. 4.7 are in better agreement with the experimental data than excess results in Fig. 4   1) The direct estimation of the bulk pressure on the outside portion of the slit pore leads to comparable excess and net adsorption profile curves.
2) Net adsorption appears to provide a more reasonable comparison to experimental data in the minimum SSA and lower temperature region.
3) This method should only be employed if a direct measurement is desired. A more precise estimation is the report NPT net adsorption results found in Fig.   4.7.

Impact of pore free volume between methane and helium molecules
In this section, the impact of assessable pore volume methane and helium probe molecules is explored. The procedure for determining assessable pore volume is presented in section 4.5.3. Since there is minimal impact of pore volume, we only present the benchmark NPT pressure method example thereby avoiding redundancy.
Pore volume occupancy information at room temperature is given in Table 4.3 where the total volume of a unit cell of a 1.42 nm graphite nanochannel pore throat is 71.428 nm 3 , methane at 66.31 nm 3 , and helium at 66.42 nm 3 . It is expected that helium pore volume is greater than methane due to the fact that helium is smaller than methane.
Another example of the difference in accessible pore volume can be seen in Table 1 of Chen et al. [2] which in their case leads to a noticeable difference in adsorption due to their choice of framework namely Na-Montmorillonite. The difference between helium and methane probing molecules is observed in Fig.   4.10. Since there is minimal difference between the two molecules on excess adsorption, we recommend using helium as the probing molecule since it is standard practice for experimentalist [14].  There is a noticeable difference between excess and net adsorption results when only considering the volume of the simulation unit cell and experimental sample as seen in Fig. 4 Key links were investigated between simulation and experimental data through a rigorous means. The experimental and simulation data must be converted using a combination of simulation and experimental techniques. Once converted, surface area data of the nanoporous material is used to normalize the data with exceptional results.
The estimation of void pockets of the porous material was shown to have a minimal dependence on the specific probing molecule used in this work.
A framework for bridging molecular information to the bulk scale is provided through simulation and modeling. This work provides a template for further molecular study of slit pores as well as meaningful information for research that seeks to capture molecular information on the bulk scale.

APPENDICES
6. All internal energies of departure in Appendices 6.1.1-6 are in units of cm 3 bar/mol and the numbers in parentheses correspond to standard deviations.