- Email: [email protected]
Microstructure and adsorption properties of organic matter in Chinese Cambrian gas shale: Experimental characterization, molecular modeling and molecular simulation
Accepted Manuscript Microstructure and adsorption properties of organic matter in Chinese Cambrian gas shale: Experimental characterization, molecular modeling and molecular simulation
Liang Huang, Zhengfu Ning, Qing Wang, Hongtao Ye, Zizheng Wang, Zheng Sun, Huibo Qin PII: DOI: Reference:
S0166-5162(18)30551-2 doi:10.1016/j.coal.2018.09.001 COGEL 3075
To appear in:
International Journal of Coal Geology
Received date: Revised date: Accepted date:
7 June 2018 31 August 2018 1 September 2018
Please cite this article as: Liang Huang, Zhengfu Ning, Qing Wang, Hongtao Ye, Zizheng Wang, Zheng Sun, Huibo Qin , Microstructure and adsorption properties of organic matter in Chinese Cambrian gas shale: Experimental characterization, molecular modeling and molecular simulation. Cogel (2018), doi:10.1016/j.coal.2018.09.001
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
ACCEPTED MANUSCRIPT Microstructure and adsorption properties of organic matter in Chinese Cambrian gas shale: Experimental characterization, molecular modeling and molecular simulation Liang Huang a,b , Zhengfu Ning a,c* , Qing Wang a
a,c
, Hongtao Ye
a,c
, Zizheng Wang c, Zheng Sun a,c, Huibo Qin d
State Key Laboratory of Petroleum Resources and Prospecting, China University of Petroleum (Beijing), Beijing
b
CR IP T
102249, China. Department of Chemical and Biomolecular Engineering, University of California, Berkeley, Berk eley, CA
94720-1462, United States.
Department of Petroleum Engineering, China University of Petroleum (Beijing), Beijing 102249, China.
d
State Key Laboratory of Heavy Oil Processing, China University of Petroleum (Beijing), Beijing, 102249, China.
AN
*Corresponding author: e-mail address: [email protected] .
US
c
ABSTRACT
M
A representative molecular model of shale organic matter (OM) is prerequisite to further theoretic investigation on the
ED
fundamental mechanisms governing storage, transport and recovery of shale gas. In this work, a systematic strategy to prepare structural and compositional properties of OM is reported first, and then a realistic molecular model of Chinese Cambrian OM is
PT
generated based on analytical experimental data. Microstructure characterization and adsorption simulation are further performed using molecular dynamics simulation and grand canonical Monte Carlo simulations, respectively. The OM model, composed of
CE
kerogen macromolecules, bitumen components and residual lighter components, shows a reasonable representation of realistic Cambrian OM with respect to structural parameters, generic compositions, physical density and porosity. The OM porous network
AC
consists of dominant ultra-micropores and limited micropores. Compared with heavier components, lighter components are more inclined to occupy accessible pores. Interestingly, lighter components are observed to appear in pairs due to competitive adsorption around heteroatom groups. Water molecules are scattered in the system because of abundant functional groups and poor pore connectivity. The OM skeleton represents the adsorption behaviors of methane, carbon dioxide and nitrogen well. The adsorption capacity is carbon dioxide > methane > nitrogen. A higher adsorption capacity corresponds to a lower pressure when the excess isotherm reaches the maximum. The adsorption behaviors of heavier hydrocarbon species (ethane, propane and n-butane) cannot be represented in the OM skeleton with ultra-micropores and limited micropores, and the effect of molecular sizes of these species cannot be neglected. This work reports a systematic construction process for realistic molecular model of shale OM, and the representative OM model can serve as a starting point to explore gas adsorption and transport mechanism in shale organic pores at
1
ACCEPTED MANUSCRIPT microscopic scale. Keywords: Organic matter; molecular model; kerogen; experimental characterization; shale; adsorption
1 Introduction Shale gas has been considered as one of the most important substitutes for conventional hydrocarbon resources due to its considerable resource abundance and high utilization efficiency (Weijermars, 2014; Yuan et al., 2015a; Paylor, 2017). However, the recovery of most shale gas reservoirs remains below 30 % because of complex gas transport mechanisms and tight reservoir
CR IP T
characteristics, as well as depleted development regime (Weijermars, 2013; Kim et al., 2017). With the rising demand for clean fossil fuels, potential technologies for enhanced shale gas recovery like carbon dioxide sequestration have attracted enormous attentions worldwide (Middleton et al., 2015; Huang et al., 2017a, 2018a, 2018b).
Shale gas in subsurface rock can be divided into adsorbed gas on various pore surfaces, free gas in fractures and pores and dissolved gas in organic matter (OM). (Wu et al., 2016a). A large proportion of shale gas is located in nanoporous nodules as a
US
form of adsorbed gas, which, for example, can take up 60-85 % in Lewis Shale (Curtis, 2002). Adsorbed gas exploitation is the key
AN
issue of enhanced shale gas recovery (Wu et al., 2015, 2016b). Adsorbed gas is traditionally believed to be distributed in organic matter and clay mineral (Curtis, 2002). However, recent investigations indicate that the contribution of clay minerals to adsorbed
M
gas content is much limited because of their hydrophilic nature, which greatly hinders the gas adsorption (Li et al., 2016a, 2016b, 2017). Accordingly, most of adsorbed gas in shale is stored in the OM nanopores, which have large specific surface areas and
ED
strong adsorption potentials (Ross and Bustin, 2009). A better understanding of gas adsorption and transport behaviors in shale OM nanopores can help improve the performance of enhanced shale gas recovery technologies. Compared with experimental
PT
measurements, molecular simulations including molecular dynamics (MD) and grand canonical Monte Carlo (GCMC) simulations
CE
are more powerful tools to gain insights into the microbehaviors of fluid mixtures in complex systems at high pressure and temperature conditions (Zhang et al., 2015). To further study the microscopic mechanisms of enhanced shale gas recovery with
AC
molecular simulation method, it is of fundamental significance to c haracterize the structure of shale OM and construct representative molecular model.
To date, various organic models related to shale OM fractions have been reported. However, most of these organic models are simplified carbon nanopore models such as carbon nanotubes (Zhu and Zhao, 2014; Yuan et al., 2015b) and graphene slits (Kurniawan et al., 2006; Kazemi and Takbiri-Borujen, 2016). Some attempts have been performed to improve the complexity of these simplified organic models. Wang et al. (Wang et al., 2016) constructed a composite pillar-layer model to represent the OM pore structure in Chinese Lower Silurian Longmaxi shale based on experimental characterization. Feng and Akkutlu (Feng and Akkutlu, 2015) built an improved carbon nanotube model with roughness by setting bumps and trenches in the nanotube surfaces. These simplified models fail to capture the heterogeneity of OM structure. Some investigations have also been conducted to consider the surface heterogeneity of organic models. Lu et al. (Lu et al., 2015) and Sun et al. (Sun et al., 2017) added various 2
ACCEPTED MANUSCRIPT heteroatom functional groups in the edge of carbon clusters to represent the chemical surface heterogeneity. Liu et al. (Liu et al., 2016b) generated physically heterogeneous surface model by adsorbing carbon adatoms on graphene. In addition to surface heterogeneity, realistic shale OM also has heterogeneous porous structure. Ungerer et al. (Ungerer et al., 2014) developed realistic kerogen molecular models with different organic types and thermal maturities based on documented experimental data by Kelemen and coworkers (Kelemen et al., 2007). The heterogeneous chemical structures of these kerogen models are constrained by analytical results, while the porous structures are generated by the phase equilibria using the MD simulations. To optimize the
CR IP T
porous network, Collell et al. (Collell et al. 2014a) inserted dummy particles into a mature type II kerogen matrix model, generating heterogeneous porous network with both artificial pores and kerogen matrix pores. Zhou et al. (Zhou et al., 2016) inserted cutter clusters instead of dummy particles into kerogen matrix models to further improve the pore shape of artificial pores. However, it should be noted that authentic shale OM is one kind of mixture of multi-components, containing bitumen components
US
and residual lighter components, as well as kerogen macromolecules. These kerogen models fail to consider the effect of OM fragments on the bulk OM skeleton structure. To our best knowledge, except in our recent work which reported the construction of
AN
a realistic OM molecular model of Chinese Silurian shale (Huang et al., 2017b), there is only one published report on bulk molecular model of shale OM. Collell et al. (Collell et al., 2014b) constructed a bulk type II shale OM model in the middle of the
M
oil formation window constrained by reported experimental data. They adopted the assumption that hydrocarbons produced in the maturation process were unexpelled, which cannot be applied to OM in gas shale with high maturity, because produced
ED
hydrocarbons are inclined to expel during the maturation process. Also, inclusion of bitumen components in their OM model was
PT
not constrained by analytical data. In our recent work (Huang et al., 2017b), we generated a bulk OM molecular model representative of Chinese Silurian shale. Nevertheless, most of structural parameters of the OM skeleton were derived from
structural parameters.
CE
experiment-based correlations rather than direct characterization results, which should be further improved with detailed analytical
AC
Laboratory experiments serve as fundamental approaches to characterize the compositional and structural properties of OM components, which are essential for the generation of molecular models of OM model (Vandenbroucke and Largeau, 2007; Bous ige et al., 2016). So far, various experimental methods including chemical approaches and physical approaches have been performed to characterize the structural parameters of OM components. Chemical approaches such as chemical degradation (Bajc et al., 2008; Budinova et al., 2009), thermal analys is (Li and Yue, 2003; Wang et al., 2010) and pyrolysis–gas chromatography– mass spectroscopy (Ru et al., 2012) require the decomposition of experimental sample. By contrast, physical approaches such as Fourier transform infrared (FTIR) spectroscopy (Zeng and Wu, 2007), Raman spectroscopy (Zeng and Wu, 2007; Zhang et al., 2017), solid-state
13
13
C nuclear magnetic resonance ( C NMR) spectroscopy (Wei et al., 2005; Tong et al., 2016) and X-ray
photoelectron spectroscopy (XPS) (Kelemen et al., 2007; Zhang et al., 2017) are nondestructive direct methods for structural characterization. Compared with destructive chemical approaches, direct physical characterization experiments can provide more 3
ACCEPTED MANUSCRIPT preserving structural information and less associated uncertainties. However, only averaged structural information can be obtained through direct characterization, and specific components and elementary compositions still need to be quantified by destructive methods (Kelemen et al., 2007). Therefore, to develop representative molecular models of complex organic samples, it is better to integrate experimental results from both indirect and direct characterization methods. In this work, the OM of Cambrian shale in Sichuan Basin, currently the most significant target zone for shale gas exploration and exploitation in China, was targeted for analysis. Based on our systematic construction strategy, we aimed to generate a realistic
CR IP T
bulk OM model containing kerogen macromolecules, bitumen components and residual lighter components by a combination of experimental characterization and molecular simulation. The chemical structure of Cambrian kerogen was characterized by elementary analysis, FTIR spectroscopy,
13
C NMR spectroscopy and XPS technique. The bitumen components were constrained
by the experimental group separation results. The residual lighter components resulted from the atomic balance assumption,
US
bitumen group separation and gas chromatography of produced gas. These OM components were integrated to generate the bulk OM molecular model using the molecular simulation method. Microstructure characterizations including steric structure and
AN
porous structure of the OM model were then conducted with the MD simulations. The gases adsorption characteristics including absolute and excess adsorption amount and isosteric heat of gas adsorption in the OM skeleton were further explored with the
M
GCMC simulations. This work provides a systematic construction process of bulk OM molecular model, which can facilitate the representative reproduce of realistic OM structure at molecular level. The proposed OM molecular model, along with its
ED
adsorption characteristics, is expected to serve as a starting point for further theoretical investigations on gas adsorption and
2 Materials and methodology
CE
2.1 Materials
PT
transport mechanism in organic pores at microscopic scale.
Two black shale outcrops were sampled from the Lower Cambrian Niuititang Formation in Sichuan Basin, Southwest China.
AC
The two samples were collected from different sampling sites located in the same district. The Lower Cambrian Niuititang Formation is deposited in a deep-water shelf environment. Pelagic algae, fungi and zooplankton are the main sources of Cambrian OM (Yang et al., 2016). A typical reservoir condition (381 K and 33.6 MPa) can be determined for the Cambrian OM according to the reported temperature and pressure gradients (Wu et al., 2014; Liu et al., 2016a; Wang et al., 2016). Detailed information on the geology, stratigraphy and resource potential for Cambrian shale has been previously reported (Yang et al., 2013).
2.2 Characterization experiments Large chunks of shale samples were crushed first, and then grounded by a ball mill. The milled powder samples were dry sieved by a 200 mesh sieve to collect particles smaller than 75 µm for subsequent measurement. A LECO CS230 carbon/sulfur analyzer was utilized to measure the total organic carbon (TOC) of investigated samples. The milled shale powders were treated with HCl to dissolve carbonates first, and then pyrolyzed up to 813.15 K. Vitrinite-like macerals 4
ACCEPTED MANUSCRIPT were targeted for thermal maturity evaluation because of the absence of vitrinite in shale of Niutitang Formation. We used a MPV-SP microscope equipped with oil-immersion objective lens and photometer to measure the reflectance of vitrinite-like materials. A reported mathematical conversion equation was then utilized to convert the measured reflectance into equivalent vitrinite reflectance (Chen et al., 2011). Kerogen samples analyzed below were isolated from the shale powders by a combination of acid treatment and Soxhlet extraction process. The shale powders were demineralized by a succession of HCL, and HF treatments. Carbonates were removed
CR IP T
by 20 wt% HCl acid at 343.15 K. Silicates were removed by a mixture of HCl and HF acids (20 wt% HCl + 40 wt% HF) at 343.15 K. Note that pyrites may not be removed from the samples, since there was no HNO 3 treatment in the acid treatment process. Therefore, it is important to check the possible effect of pyrites on kerogen characterization. The remains after each acid treatment were filtered and washed with deionized water and oven-dried at 353.15 K. Thereafter, Soxhlet extraction with toluene was utilized
US
to extract soluble bitumen components. The extracted bitumen components were then treated with chromatography column for group separation. The obtained kerogen samples were oven-dried under vacuum first, and then treated by a Thermo Scientific
In this work, FTIR spectroscopy,
13
AN
FLASH 2000 CHNS-O Analyzer for elemental analysis.
C NMR spectroscopy and XPS techniques were utilized to characterize the chemical
M
structure of the studied kerogen samples. The FTIR spectra were measured by a MAGNA-IR56 FTIR spectrometer. Sample powders of 0.5 mg were mixed with 200 mg dry KBr powders, then pressed under vacuum. FTIR spectra measurement was
13
C NMR spectra were measured at 100.64 MHz by a Bruker AVANCEIII TM-600 MHz spectrometer with the MAS magic
PT
The
ED
performed at ambient temperature, and spectra were recorded between 4000 and 400 cm -1 with 50 scans at a resolution of 4 cm -1.
angle spinning and using a ZrO2 rotor (4 mm diameter). A recycle delay time of 3 s and a contact time of 4 ms were adopted in the
CE
cross-polarization process. The XPS measurement on kerogen powders was conducted on a Thermo Scientific ESCALAB 250Xi 2
spectrometer using monochromatic Al Kα radiation (1486.6 ev). Ar+ scanning gun was used on an area of 8×8 mm . The sputtering
AC
rate was 4 nm/min, and the transmit current was 10 mA. To specify and quantify the functional groups and elemental forms, Lorentz peak-fitting was performed on the measured spectra curves.
2.3 Molecular modeling method In this section, the construction process of average molecular model for Cambrian OM is discussed. The OM model is integrated from individual fractions, including kerogen, bitumen fractions and lighter components. Kerogen. Kerogen has complex macromolecule structure. It is composed of polyaromatic units with cross-linkages such as alkyl, ether and sulfur bridge. Kerogen constitutes the basic skeleton of OM, and can represent most of the physical and chemical 13
characteristics of OM. In this work, we adopt a combination of elementary analysis, FTIR, C NMR and XPS techniques to obtain detailed analytical data of elemental compositions, structural parameters and functional groups for Cambrian kerogen. These analytic data are utilized to derive the kerogen chemical structure. The building procedure of the kerogen model follows the 5
ACCEPTED MANUSCRIPT process in our previous work (Huang et al., 2017b). In this study, carbon structural parameters of Cambrian kerogen are derived from the
13
C NMR spectrum. Only organic
carbon forms including aromatic carbons and aliphatic carbons are considered. The functional group forms of O, N and S are analyzed based on the XPS results. For oxygenated groups, we exclude these components in residual minerals and water. In particular, carbonyl groups, alcoholic groups and ether groups are considered. For nitrogen-bearing groups, pyridine, pyrrole, amine, as well as nitrogen oxide are considered. In the case of sulfur-bearing groups, we not only focus on the organic sulfur
CR IP T
groups including thiophene, sulfoxide and sulfoether, but also check the possible influence of pyrite minerals. The FTIR spectra are utilized to qualitatively characterize the functional groups of Cambrian kerogen, which can serve as a self-consistent check 13
with the results of C NMR and XPS techniques.
The kerogen structure topology is largely inspired from previous kerogen fractions in the literature (Shinn, 1984). A unit of
US
100-260 carbon atoms for each kerogen molecule has been reported as the minimized size requirement to capture the major structural features of kerogen (Orendt et al., 2013; Ungerer et al., 2014). A total number of 210 carbon atoms for single Cambrian
AN
kerogen molecule is adopted to meet the size requirement. In addition to linear alkyl fragments, branched alkyl, oxygen-, sulfurand nitrogen-bridge segments are also introduced as cross-links between aromatic units. Moreover, apart from aromatic ring units,
M
monocyclic and polycyclic saturated rings are also added to match the H/C ratio of Cambrian kerogen. The initial version of kerogen molecule is designed to be a flat-shape structure with minimized extension so as to approach the equilibrium conformation.
ED
Thereafter, adjustments on the functional group locations and polyaromatic unit shapes are performed to improve the consistency
PT
between the simulated and experimental 13 C NMR spectra. The pcff force field is adopted throughout this work. We adopt the pcff force field to assign atomic types and charges for kerogen atoms so as to prepare the kerogen unit for further structure relaxation.
CE
The reason that the pcff force field is used will be discussed later. Bitumen. Bitumen is defined as the component of OM that can be soluble in dichloromethane (Vandenbroucke, 2003). It can
AC
be divided into asphaltene, resin, aromatic hydrocarbons, and saturated hydrocarbons. The asphaltene and resin molecules are composed of polyaromatic units with polar heteroatoms. Their structures are similar to that of kerogen, but with a lower molecular weight. Asphaltene is soluble in polar solvent but insoluble in n-alkanes, while resin is soluble in liquid n-alkanes. Compared with asphaltene, resin has a lower aromaticity with more polar fractions and alkane tails. The structure sizes of asphaltene and resin are determined to ensure an economic and reasonable simulation system, as reported in our latest work (Huang et al., 2017b). To be consistent with the conclusion of Yen (Yen, 1979), the nitrogen group is represented by quinoline unit, while the oxygen group is in the form of ether unit. In the case of aromatic hydrocarbons and saturated hydrocarbons, a limited number of representative components are targeted in this study. Toluene is chosen as the representative aromatic component, as adopted in Collell’s work (Collell et al., 2014b). While fichtelite, a common tricyclic alkane in bitumen, is chosen as the representative saturated component, since most of long 6
ACCEPTED MANUSCRIPT aliphatic components disappear in the late catagenesis stage. Also, some methane molecules are included in the saturated components to represent residual lighter alkyl components. The specific molecular sizes and quantities of representative components are adjusted to match the experimental bitumen group separation results in the Cambrian OM. Lighter components. In this work, methane, ethane, carbon dioxide, nitrogen and water are defined as lighter components. These lighter components are assumed to be generated from kerogen during the maturation evolution. By adopting the atomic balance assumption (Loucks and Ruppel, 2007), namely there is no atomic loss during the maturation process, we determine the
CR IP T
initial total lighter components under reservoir conditions. Detailed analytical data of the immature type II kerogen in Duvernay shale (Kelemen et al., 2007) are selected for atomic balance analysis. However, it should be noted that most of these lighter components are inclined to be expelled from the OM, especially for gas shale in the late catagenesis stage. In this work, we first determine the compositions of lighter components for deduction by gas chromatography of produced gas ( Wu et al., 2014), and
US
then determine the quantity for deduction by restricting the deducted methane molecules . The deducted quantity of methane molecules from the initial total lighter components is determined by adjusting the residual methane quantity to match the
AN
experimental weight fraction of saturated hydrocarbons in bitumen. Note that we include methane in both the bitumen components and the initial lighter components. On one hand, we include methane in the saturated components of bitumen to represent resid ual
M
lighter alkyl components. This is also intended to match the experimental results of bitumen group separation in the Cambrian OM. On the other hand, we include methane in the initial total lighter components so as to meet the atomic balance assumption.
ED
However, it should be noted that methane molecules in the initial lighter components are deducted later. Thus in the final OM
PT
model, there are only methane molecules included in the bitumen components. OM. Kerogen molecules, bitumen fractions (asphaltene, resin, aromatic hydrocarbons and saturated hydrocarbons) and
CE
remaining lighter components (ethane, carbon dioxide, nitrogen and water) are integrated to generate the bulk molecular model of Cambrian OM. A total of 1700 atoms and a box size of 25 Å of bulk molecular model are considered as the minimized size
AC
requirements (Ungerer et al., 2014). To meet the size requirements, 10 kerogen molecules are selected to form the skeleton of OM model. Then a suite of bitumen fractions are introduced in the system according to the experimental bitumen group separation results. The total lighter components derived from the atomic balance assumption are introduced to construct the initial version of OM model. Thereafter, a large amount of lighter components with their compositions constrained by the gas chromatography results are deducted from the initial simulation system to generate the realistic molecular model of Cambrian OM under typical reservoir conditions.
2.4 Simulation details 2.4.1 Structure optimization In this work, a combination of geometry optimization and MD simulation is utilized to generate the bulk OM model using the Materials Studio software with pcff force field. This force field describes the dispersion-repulsion interactions using 7
ACCEPTED MANUSCRIPT Lennard-Jones (LJ) 6-9 function, and describe the Coulombic interactions using force field assigned point charges (Sun et al., 1994). This ab initio all-atom force field has been parameterized, tested and validated for application to various organic materials (Sun, 1995), and has performed well for predicting the volumetric and thermodynamic properties of kerogen (Yiannourakou et al., 2013; Ungerer et al., 2014). The simulation accuracy generally increases with the cutoff distance for short-range pair interaction. However, it is well known that the cutoff distance cannot exceed the half length of the minimum dimension of the simulation box so as to meet the minimum image convention. Since the cell length of our organic models ranges from 32.04 to 34.5 Å, we adopt a
CR IP T
fine cutoff distance of 15.5 Å to truncate the short-range LJ interactions. The Ewald summation method is applied to handle the long-range electrostatic interactions. The OM fractions are first geometry optimized by the smart minimization algorithm with a fine convergence criterion. The methodology of geometry optimization refers to the work of Ru et al. (Ru et al., 2012). A max iteration step of 5000 is adopted to guarantee the energy convergence. Thereafter, these optimized OM fractions are enclosed in an
US
amorphous cell with periodic boundary conditions to generate the initial c onfigurations of OM model. To investigate the uncertainties caused by different initial configurations, 10 configurations are generated in this stage for further discussion. The cell 3
AN
density is first targeted at 0.1 g/cm (Collell et al., 2014b; Ungerer et al., 2014; Huang et al., 2017b, 2018a, 2018b ), with the system temperature being 300 K.
M
Subsequently, these configurations are optimized with the annealing strategy to obtain the global lowest energy configurations (Ru et al., 2012). A total of 10 annealing cycles of dynamic runs are performed for each configuration. Each complete annealing
ED
cycle consists of 5 heating ramps and 5 cooling ramps, for which the target temperature is gradually increased from the initial
PT
temperature (300 K) to a mid-cycle temperature (800 K) and then decreased again to the initial temperature. The dynamic time step is 1 fs, and a total simulation time of 2000 ps is performed for the annealing dynamics run. The simulation is carried out using
CE
the NVT canonical ensemble (fixed molecular number, box volumes and system temperature) with the temperature controlled by the Nosé-Hoover thermostat (Hoover, 1985). After the annealing treatment, the output configurations are further geometry
AC
optimized.
These annealed configurations are then relaxed by a succession of MD simulations, as adopted in previous molecular simulation studies (Collell et al., 2014b; Huang et al., 2017b, 2018a, 2018b). NVT ensemble is first adopted to relax these configurations at 800 K. Then succeeding NPT ensemble (fixed molecular number, system pressure and temperature) simulations are performed at 33.6 MPa and gradually decreasing temperature of 800, 600, 400 and 300 K. The optimized configurations at each temperature are output as the initial structures for the next NPT run at a lower temperature. The pressure is imposed using the Berendsen barostat (Berendsen et al., 1984), while the temperature is imposed with the Nosé-Hoover thermostat (Hoover, 1985). The simulation time is 400 ps for the NVT simulation, and 500 ps for the NPT simulations. Collell et al. (Collell et al., 2014b) adopted 12.5~250 ps for each NVT simulation and 200~400 ps for each NPT simulation. Ungerer et al. (Ungerer et al., 2014) adopted 100 ps for each NVT simulation and 100~200 ps for each NPT simulation. The simulation time in this work is much 8
ACCEPTED MANUSCRIPT longer than these in previous work (Collell et al., 2014b; Ungerer et al., 2014). The simulation time of each run has also been checked to be long enough to ensure the density convergence. At last, these configurations are further relaxed at typical reservoir conditions (381 K and 33.6 MPa) for 1000 ps to obtain the equilibrium structures.
2.4.2 Pore structure characterization In this work, free pore volumes and specific surface areas are computed with the probe insertion method (Connolly, 1983). Probe molecules with fixed Connolly radius are randomly inserted into the system. These probes roll over the atomic van der
CR IP T
Waals surfaces to determine the solid skeleton surfaces, and the regions wrapped by the skeleton surfaces are identified as free pore volumes. Similarly, a series of spherical probes of stepwise increasing Connolly radius are adopted to detect the corresponding free pore volumes and specific surface areas. Accordingly, the pore volume or specific surface area distributions of organic models are computed by differentiating the free pore volumes or specific surface areas with respect to different probe
US
diameters.
2.4.3 Adsorption behavior simulation
AN
The GCMC simulations were performed to characterize the adsorption behaviors in the OM skeleton at lighter gas species (methane, carbon dioxide and nitrogen) and heavier hydrocarbons species (ethane, propane and n-butane), with the force field,
M
pair-wise interaction, cutoff distance, thermostat method consistent with that for MD simulations. These lighter gas species are dominate species in Cambrian shale gas according to the gas chromatography results (Wu et al., 2014). The adsorption of these
ED
heavier hydrocarbons species is investigated to analyze the possible effect of adsorbate size. Adsorbate molecules are adopted as
PT
full-atom rigid models. A total of 4×107 Monte Carlo steps are performed for each pressure point, wherein the first 2×10 7 equilibrium steps are utilized as relaxation, while another 2×10 7 production steps are used as statistics. The GCMC simulations
CE
outputs the absolute adsorption amount, which can be further converted to excess adsorption capacity so as to compared with experimental data by,
(1)
AC
ne na
Where na is the absolute adsorption amount, ne is the excess adsorption amount,
is the gas density calculated by the
Peng-Robinson equation, and is the free pore volume. The isosteric adsorption heat of a component is computed to reveal the pore filling mechanism, which is defined as the partial molar enthalpy of the sorbate component in the reservoir minus that in the framework, which is determined by,
Qst RT intra N E
(2)
Where Qst is the isosteric adsorption heat, R is universal gas constant, T intramolecular chemical potential,
N
is the system temperature, intra
is the ensemble averaged molecular number of sorbate, 9
E
is the
is the ensemble averaged
ACCEPTED MANUSCRIPT total energy for the sorbate-sorbent configuration.
3 Results and discussion 3.1 Chemical structure characterization 3.1.1 FTIR Spectroscopy The TOC content for shale sample Y1 and Y2 is 2.3 wt.% and 2.1 wt.%, respectively, indicating a good organic richness for Cambrian shale (Chalmers and Bustin, 2007). The equivalent vitrinite reflectance is 1.58 % and 2.01 % for sample Y1 and Y2,
CR IP T
respectively. The thermal maturity of studied Cambrian shale can be mainly classified as high mature level. The Van Krevelen diagram is shown in Fig. S1. The sample points of Cambrian shale are much closer to the type II curve, which indicates the source rocks mainly contain type II kerogen.
In this work, shale sample Y1 and Y2 and their isolated kerogen sample K1 and K2 were chosen for FTIR spectrum analysis
US
so as to characterize the chemical functional groups. Standard FTIR absorption spectra of these samples are displayed in Fig. 1. The band assignments of investigated functional groups are listed in Table S1 (Ojima, 2003; Chen et al., 2012).
K2
1620 C=C 1455 CH2,CH3
AN
1700 C=O
1286 C-OH 429 FeS2
M
Absorbance
1375 CH3
3680 K+I
ED
K1 Y2
3500
3000
AC
4000
613 A 1150 Q+A 707 CA 512 Q 912 K
CE
PT
Y1
1059 -O- 900-730 CH x
2500 2000 1500 -1 Wavenumber (cm )
1000
500
Fig. 1. FTIR spectra of Cambrian shale and kerogen samples. Symbol description: K denotes kaolinite, Q denotes quartz, CA denotes calcite, I denotes illite, A denotes anhydrite. As seen in Fig. 1, kerogen sample K1 and K2 show basically consistent absorption spectrum, which facilitates the generation of an average molecular structure of Cambrian kerogen. Compared with the kerogen spectra, the shale spectra show similar trend, but with different peak intensity. Also, a shift of peak sites toward the lower wavenumber can be observed for the shale spectra. This peak shift can be attributed to the mineral components in the shale samples. As seen on the shale spectra, a large number of peaks related to minerals such as kaolinite, quartz and calcite can be distinguished. The shale spectrum results from the overlay of absorption peaks for both OMs and minerals. The absorption peaks for organic groups become insignificant because of the 10
ACCEPTED MANUSCRIPT influence of various minerals. Most of the absorption peaks of minerals disappear for the kerogen spectra, indicating that acid treatment can remove most of the minerals in shale. However, we can still observe the absorption peaks for pyrites on the kerogen spectra. These adsorption peaks qualitatively shows the occurrence of pyrite minerals in the kerogen samples. This observation has also been previously reported, since pyrites cannot been removed by HCl/HF acid treatment (Kelemen et al., 2007). The kerogen spectra show significant absorption peaks of organic groups with fairly correct peak sites. The peaks at 1700, 1286 and 1059 cm–1 correspond to the oxygenated groups, the peaks at 1455 and 1375 cm –1 correspond to the saturated alkyl –1
CR IP T
groups, while the peaks at 1620 and 900-730 cm correspond to aromatic groups. These peaks indicate that the chemical structure of Cambrian kerogen contains aromatic components, saturated aliphatic components and oxygen-bearing groups. The kerogen –1
spectra show remarkable strong peaks of aromatic C=C ring stretching at 1620 cm and aromatic CHx out-of-plane deformation at 900-730 cm–1, suggesting that aromatic carbon is the major carbon type in Cambrian kerogen. Also, the peak intensities of the CH3
US
asymmetric bending vibration and the CH2 shear vibration at 1455 cm–1 and the CH3 symmetric bending vibration at 1375 cm–1 are fairly strong, indicating the existence of CH3 and CH2 groups in the kerogen structure. However, the kerogen spectra show an -1
AN
absence of absorption peak of vibration of long-chain alkane ((CH2 )n , n≥4) at 720 cm , which indicates that the aliphatic carbons are mainly in the forms of short side chains or linkages. The relative peak intensities between these oxygenated groups indicate
M
that ether groups and carbonyl groups are the main oxygen-containing groups, and there is nearly no alcohol group in the kerogen structure.
13
C NMR spectroscopy was performed on kerogen sample K1 and K2 to characterize the carbon
PT
In this work, solid-state
ED
3.1.2 13 C NMR Spectroscopy
skeleton structure of Cambrian kerogen. The measured kerogen spectra are presented in Fig. 2. The kerogen spectra are composed
CE
of aliphatic carbon region (0-90 ppm), aromatic carbon region (90-165 ppm) and carbonyl carbon region (165-220 ppm). Compared with the aliphatic carbon region and carbonyl carbon region, the aromatic carbon region has the remarkably largest peak
AC
area. We can derive that aromatic carbon is the major structural fraction of Cambrian kerogen, while aliphatic and carbonyl carbon mainly serve as linkages and short side chains, which is in good agreement with the FTIR results.
11
ACCEPTED MANUSCRIPT H
125 far 135 far 145 far 156,166 far
C M
D
40,48 fal 20 f A al
O
99,105 far
A
15 fal
K2
O O
50-90 fal
H
33 fal
B
25 fal
CR IP T
175,183 fa
B
K1
250
200
150
100
50
0
US
Chemical Shift (ppm)
-50
Fig. 2. 13 C NMR spectra of Cambrian kerogen samples. Symbol description and peak assignment can refer to Table 1.
AN
To quantitatively specify the carbon groups in the chemical structure, Lorentz peak-fitting of was conducted on the 13 C NMR spectrum of kerogen sample K2. The fitting curves are shown in Fig. S2. The structural parameters and chemical shifts have been
M
reported in the work of Trewhella et al. (Trewhella et al., 1986). The fitting area fractions are listed in Table 1. The important
ED
lattice parameters have been previously defined (Kelemen et al., 2007) and are shown in Table 2. Table 1
Structural
PT
Carbon structural parameters from 13 C NMR spectrum of kerogen sample K2. Chemical
Area fraction
Carbon type
Location
shift (ppm)
f al M
15
Aliphatic methyl
Cal -C*H3
5.29
f al A
20
Aromatic methyl
Car-C*H3
1.68
f al B
25
Aliphatic C(2) carbon
C*H2 -CH3
0.43
f al H
33
Methylene
C*H2
0.71
f al D
40, 48
C*, C*H
5.97
AC
CE
parameter
(%)
Methine, Quaternary sp 3 carbon
55
Oxy-methylene
C*H2 -O
4.45
62, 69
Oxy-methine
C*H-O
10.49
83, 89
Oxy-aliphatic ring
Cal *-O
1.85
f ar A
99, 105
Oxy-aromatic ring
Car*-Car-O
7.18
f ar H
125
Protonated aromatic carbon
Car*-H
29.50
f al O
12
ACCEPTED MANUSCRIPT f ar B
135
Bridged aromatic carbon
Car*-Car
8.49
f ar C
145
Branched aromatic carbon
Car*-Cal
13.74
f ar O
156, 166
Oxy-aromatic carbon
Car*-O
4.04
faO
175, 183
Carboxyl, carbonyl
RC*=OR
6.17
Table 2 Lattice parameters derived from 13 C NMR spectrum of kerogen sample K2 Lattice
Description
Value
f al
f al M + f al A + f al B + f al H + f al D + f al O
Aliphaticity
30.87 %
f ar
f ar A + f ar H + f ar B + f ar C + f ar O
Aromaticity
62.95 %
f al B + f al H )/ f ar C
Aliphatic chain length
0.07
f ar B /( f ar A + f ar H + f ar C + f ar O )
Degree of condensation
f al D +C*H-O)/ f al
Aliphatic branch
0.53
f ar H / f ar
Aromatic substitution
0.53
BI
(
1-
US
AN
X BP
(
0.16
M
Cn
CR IP T
Definition
parameter
ED
As seen in Table 2, the aromaticity of kerogen is up to 62.95 %, 10 benzenes or 6 naphthalenes or other similar structures can be derived for each 100 carbons in Cambrian kerogen. This indicates aromatic carbons constitute the main skeleton of Cambrian
PT
kerogen structure. In addition to aromatic carbons and aliphatic carbons, the kerogen structure also contains some carboxyl or carbonyl groups, which can serve as branched side chains or short linkages between aromatic clusters. The aliphatic chain length is
CE
0.07, indicating that the aliphatic carbons mainly exist in the form of short branched chains. This conclusion can be further
AC
validated by the negligible FTIR absorption peak at 720 cm -1. The degree of condensation ( X BP ) is closely related to aromatic cluster size. The derived X BP at 0.16 for Cambrian kerogen is smaller than that of naphthalene (0.25), which indicates that the condensation degree for aromatic cluster can be represented by 1-2 aromatic rings for Cambrian kerogen. The lattice parameter for aliphatic branch ( BI ) is 0.53, suggesting that the aliphatic carbons are highly branched for Cambrian kerogen. The aromatic substitution ( =0.53) indicates that 3-4 carbons in each aromatic ring are substituted by neighboring groups.
3.1.3 XPS results In this work, the elementary forms for O, N and S in Cambrian kerogen were characterized using the XPS technique. The XPS spectra of kerogen sample K1 and K2 are illustrated in Fig. 3. The two spectra are fairly similar to each other, and the XPS spectrum of K2 is utilized to quantify the elementary forms in Cambrian kerogen. The derived element concentrations for kerogen
13
ACCEPTED MANUSCRIPT K2 from XPS spectrum are listed in Table 3. In addition to C, O, N and S, Cambrian kerogen also contains a small trace of Cl and Fe. C is the major element in the kerogen structure, with its area concentration being 89.17 %. The absorption area concentration for O is 7.32 %, much larger than that for N and S. The observed Cl derives from the remains of acid treatment, while Fe is the main component of pyrite which has also been identified in the FTIR spectra. However, since the area concentrations for both Cl
US
Cl 2s S 2s Cl 2p S 2p
800
600 400 Binding Energy (ev)
200
0
M
1000
AN
K1
O 2s
N 1s
K2
Fe 2p
O KLL
O 1s
CR IP T
C 1s
and Fe are lower than 0.2 %, their effect on the kerogen chemical structure can be neglected.
ED
Fig. 3. XPS spectra of Cambrian kerogen samples. Table 3
PT
Element concentration derived from XPS spectrum of kerogen sample K2 Binding energy(eV)
Area concentration %
Fe 2p
714.20
0.11
O 1s
531.80
N 1s
401.90
C 1s
284.50
89.17
Cl 2p
200.30
0.13
S 2p
168.80
1.61
AC
CE
Element
7.32 1.66
To quantitatively specify the elementary forms for heteroatom including O, N and S in Cambrian kerogen, the XPS spectrum for K2 sample was fitted by the XPS peak-fitting software, and the results are shown in Table 4 and Fig. S3. The binding energies corresponding to functional groups have been documented in previous work (Tong et al., 2011; Li et al., 2014). Due to the effect of N and S, the peak-fitting of O 1s cannot accurately determine the binding modes between C and O. Therefore, the spectrum of C 1s signal was fitted to quantify the elemental forms of O associated with C atom. The results indicate that ether or alcohol groups are the main oxygen-containing groups in the kerogen structure. Also, the kerogen contains some carbonyl and carboxyl groups.
14
ACCEPTED MANUSCRIPT The analysis on N 1s signal indicates that amine groups are the dominant nitrogen-containing groups in the kerogen structure. Besides, the Cambrian kerogen contains a small amount of stable pyridine and pyrrole groups. The nitrogen oxide groups introduced from the acid treatment are observed to be well removed. The fitted results on S 2p signal indicate that S mainly exists in the organic components, and sulfoxide group is the dominant sulfur-containing group, with its area ratio being up to 86.43 %, followed by the thiophene group at 13.57 %. Table 4
Binding energy (eV)
Peak area ratio (%)
C 1s
C-C, C-H
284.51
79.70
C-O, C-OH
285.98
11. 79
C=O
287.05
3.13
COOH, COOR
288.45
2.37
Pyridine
398.09
15.80
Pyrrole
400.05
7.46
Amine
401.45
73.38
NOx
403.01
3.35
Sulfoether
162.13
0.00
Pyrite
163.10
Thiophene
164.09
Sulfoxide
168.97
0.00
13.57 86.43
PT
3.2 Molecular models
AN
S 2p
ED
N 1s
CR IP T
Functional group
US
Element
M
Elementary forms for heteroatom derived from XPS spectrum for K2 sample
3.2.1 Molecular structure of Cambrian kerogen
CE
Based on the analytical structural parameters, an initial 2D molecular structure of Cambrian kerogen was generated. The
13
C
NMR spectrum was predicted using the ACD/CNMR Predictor software, and compared with the experimental spectrum. The
AC
aliphatic carbons and heteroatom groups on the kerogen structure were continuously adjusted until good agreement was observed for the predicted and experimental spectrum (Fig. 4). The final representative Cambrian kerogen structure (C210 H184 O20 N4 S4) is illustrated in Fig. 5, with its compositions and structural parameters presented in Table 5. The kerogen atomic compositions are in good agreement with the elemental analys is results. The kerogen skeleton consists of 9 naphthalenes and 6 benzenes linked by branched aliphatic carbons and heteroatom groups. The lattice parameters of the kerogen structure such as degree of condensation, aliphatic branch and substitution are fairly close to the experimental
13
C NMR results. The oxygenated groups are mainly
composed of ether, carbonyl and carboxyl groups, consistent with the XPS and FTIR results. Moderate deviations are observed for N and S groups due to their low contents. The 4 N atoms are distributed in 3 amine groups and 1 pyridine group, while the 4 S atoms are involved in 3 sulfoxide groups and 1 thiophene group. In this work, 10 kerogen units are adopted to build the skeleton of
15
ACCEPTED MANUSCRIPT the Cambrian OM model, meeting the reported minimized size requirement for 1700 atoms per unit cell (Collell et al., 2014b).
006 Simulation Experiment
004
CR IP T
002
000 180
160
140
120 100 80 Chemical Shift (ppm) 13
40
C NMR spectrum of Cambrian kerogen K2.
PT
ED
M
AN
Fig. 4. Predicted and experimental
60
US
200
Table 5
CE
Fig. 5. Molecular structure of Cambrian kerogen.
Property
AC
Compositional and structural parameters of Cambrian kerogen
a
Structural parameter
Experimental
Modelling
H/C
0.87
0.88
Atomic
O/C
0.093
0.095
compositions
N/C
0.019
0.019
S/C
0.018
0.019
fal
30.87 %
36.19 %
f ar
62.95 %
61.43 %
Cn
0.07
0.18
X BP
0.16
0.16
C groups
16
20
0
ACCEPTED MANUSCRIPT
N groups
S groups
a
0.53
0.53
0.53
0.50
Carbonyl (mol % C~O)
15.89
17.64
Ether (mol % C~O)
59.99
58.82
Carboxyl (mol % C~O)
24.12
23.53
Pyridine (mol % N)
15.80
25.00
Pyrrole (mol % N)
7.46
0.00
Amine (mol % N)
73.38
75.00
NOx
3.35
0.00
Sulfoether (mol % organic S)
0.00
0.00
Thiophene (mol % organic S)
13.57
25.00
Sulfoxide (mol % organic S)
86.43
75.00
CR IP T
O groups
BI
Experimental ato mic co mpositions are derived fro m elemental analysis. Experimental C groups are derived fro m 13 C NMR spectrum.
Experimental O, N and S groups are derived from XPS spectrum.
US
3.2.2 Structures and compositions of bitumen components
The bitumen components are composed of representative components, with their compositions listed in Table 6. Note that
AN
resin is a derivative of higher plant. The OM in shale from the Lower Cambrian Niuititang Formation is mainly derived from
M
pelagic algae, fungi and zooplankton. Thus the content of resin is relatively low in our investigated OM. However, we integrate the resin in our OM system for the purpose of generating a complete bulk molecular model of OM. The molecular structures of resin
ED
unit and fichtelite unit are shown in Fig. 6. The molecular size (C33 H33 ON) of our resin unit refers to the documented one (C33 H41 O1.6 N0.4 S0.1) in the work of Tissot and Welte (Tissot and Welte, 1984). The topology including aromatic rings and alkane
PT
tails of our resin unit refers to the reported resin structure in the work of Collell et al. ( Collell et al., 2014b). The molecular
CE
structures and quantities of bitumen components are adjusted to match the experimental weight fraction. Good agreement between the simulated and experimental weight fraction distribution can be observed. The converted total bitumen content (0.0547 wt. %)
AC
is also fairly close to the experimental result (0.0511 wt. %). In this work, asphaltene components are excluded because of their low content, and saturated hydrocarbons and resin components are the dominant components in Cambrian bitumen.
17
ACCEPTED MANUSCRIPT (b)
CR IP T
(a)
Fig. 6. Molecular units of resin and fichtelite in Cambrian shale. (a) resin unit (C33 H33 ON); (b) fichtelite unit (C14 H24 ). Atom representation: black for C atom, red for O atom, blue for N atom and grey for H atom.
Weight composition of bitumen components in Cambrian OM
US
Table 6
Experimental data a
Representative
Molecule
Weight fraction
components
molecule
numbers
(%)
(%)
resin
resin
1
42.94
40.34
asphaltene
asphaltene
0
0
6.45
aromatics hydrocarbons
toluene
1
8.61
5.75
fichtelite
1
48.46
47.46
methane
20
a
M
ED
PT
saturated hydrocarbons
AN
Bitumen
Experimental results from bitumen group separation on Cambrian sample Y2.
CE
3.2.3 Molecular model of Cambrian OM
Associated with kerogen units and bitumen components, lighter components are also included to generate the initial Cambrian
AC
OM model adopting the assumption of atomic balance. The component distribution in the initial organic model is presented in Table S2. The validation of atomic balance is shown in Table 7. Reasonable consistencies for atomic compositions are observed between the initial organic model and an immature type II kerogen of Duvernary shale (Kelemen et al., 2007). The final Cambrian OM model is obtained by deducting a suite of lighter components prone to release from the organic skeleton. The compositions of deducted lighter components are restricted by gas chromatography results on Cambrian shale gas (Wu et al., 2014) (Table 8). The deducted quantity of methane is determined by adjusting the residual methane quantity to match the experimental weight fraction of saturated hydrocarbons in bitumen. The component distribution and atomic distribution of the final Cambrian OM model are shown in Table 9 and Table S3, respectively. Kerogen molecules dominate the compositions in the OM model. The final OM model corresponds to type II organic type (Fig. S1).
18
ACCEPTED MANUSCRIPT Table 7 Atomic balance between initial Cambrian OM model and immature kerogen phase Initial OM model Immature kerogen a
a
H/C
O/C
N/C
S/C
1.17
0.097
0.029
0.017
1.17
0.097
0.029
0.014
Analytic data of an immature kerogen of Duvernary shale documented by Kelemen et al. (Kelemen et al., 2007).
Table 8 Experimental and deducted lighter components distribution Molecular number
Mole fraction (%)
Experimental data a (%)
methane
194
96.52
96.52
ethane
1
0.50
0.35
nitrogen
2
1.00
1.24
carbon dioxide
4
1.99
1.75
Analytic data of gas chromatography on produced gas from Cambrian shale documented by Wu et al. (Wu et al., 2014).
US
a
CR IP T
Representative molecule
Components distribution in the final Cambrian OM model
AN
Table 9
Molecular number
Weight fraction (%)
kerogen
10
94.56
resin
1
1.35
toluene
1
0.27
fichtelite
1
methane
20
nitrogen
12
carbon dioxide
6
ED 0.58
PT
0.94
CE
water
M
Representative molecule
10
0.99 0.78 0.53
AC
3.3 Microstructure characterization 3.3.1 Properties of steric structure
In this section, the steric molecular model of Cambrian OM (Fig. 7) was constructed based on the compositions in Table 9 and the simulation details for structure optimization in section 2.4.1. As seen in Fig. 7, the OM structure is highly unordered and no remarkable stacking of aromatic units can be observed. This observation can be attributed to the chemical structure of Cambrian kerogen, as shown in Fig. 5. Cambrian kerogen contains a small size of polyaromatic clusters, which can be represented by 1 or 2 benzene rings. These benzene rings are connected by bendable aliphatic branches and heteroatom groups, which result in a strong steric hindrance on the stacking of aromatic units. Interestingly, we observe that lighter components such as water, methane, nitrogen and carbon dioxide tend to appear in pairs around heteroatom groups, showing competitive adsorption behaviors due to surface chemistry heterogeneity. Instead of aggregating into clusters, water molecules are scattered around hydrophilic groups, 19
ACCEPTED MANUSCRIPT which can be attributed to the abundant polar groups and poor inter-pore connectivity. On one hand, water molecules are preferentially adsorbed on hydrophilic groups on the OM skeleton. The OM skeleton contains sufficient polar groups to accommodate these water molecules considering the low moisture content. Thus, water molecules are prone to locate around thes e polar groups. On the other hand, the OM skeleton is highly unordered, and the effective inter-pores are poorly connected, which restricts the transport of water molecules to aggregate into clusters. The random distribution of bitumen components and lighter
M
AN
US
CR IP T
components aggravates the disorder of pore network in the OM system.
ED
Fig. 7. View of simulation box of the OM models at 381 K and 33.6 MPa. Molecule representation: cyan for resin, grey for fichtelite, red for toluene, green for methane, purple for water, blue for nitrogen, yellow for carbon dioxide.
PT
3.3.2 Properties of porous structure
As shown in Table 10, the averaged physical densities of our organic models are in the range of 1.22±0.01-1.30±0.03 g/cm3,
CE
fairly close to the reported experimental data (1.18-1.25 g/cm3) of high mature type II kerogen (Okiongbo et al., 2005). The helium probed porosities of these organic models range between 10.95±0.3 and 21.21±0.4 %, which are consistent with the calculated OM
AC
porosities (4.45-22.50 %) in Barnett mudstone based on SEM images (Loucks et al., 2009). In addition to the original annealing temperature cycle (300~800~300 K), a new temperature cycle (300~900~300 K) is adopted to relax the OM structure to account for possible temperature cycle effect. The mid-temperature is increased to be 900 K, consistent with the documented highest temperature for OM structure relaxations (Collell et al., 2014b; Ungerer et al., 2014). The OM structure under this temperature cycle is referred as OMe in Table 10. The comparison results show that the fluctuations of OM structure properties under the two different temperature cycles are tiny and acceptable. The unordered combination of kerogen molecular units constitutes the kerogen matrix, which is also defined as the OM skeleton. These lighter species and heavier species are included to represent the residual components in the porous network of OM skeleton. In Table 10, we present the effect of the inclusion of these species on the OM structure properties. The result shows that the densities increase while the porosities decrease for our organic models with the inclusion of these species. This observation suggests that the porous network of OM is mainly constituted by the kerogen skeleton, and these inclusions occupy some pore volumes, contributing to the increasing condensation of the OM system. In Table 10, we also present the effect of probe size on the structure properties of organic models. The porosities and specific surface areas probed by helium are far smaller than these of innate porous structures with a probe diameter of zero. Note that the specific surface 20
ACCEPTED MANUSCRIPT 2
areas exceed 4000 m /g, and the porosities can be up to 40 % under a probe diameter of zero. These high porosities and specific surface areas are mainly contributed by tiny pores, which cannot be accessed by gas molecules. These tiny pores cannot be measured in kerogen samples in laboratory. However, these tiny pores should not be neglected since they may be expanded due to kerogen swelling. When the probe size increases to that of helium, the porosities and specific surface areas show a significant decrease. The porosities decline to comparable values with experimental results, while the specific surface areas (1387-1524 m2 /g) are still higher than measured results. The higher specific surface area can be attributed to the narrow pore size distributions of our organic models. The pore size in actual kerogen matrix ranges from several angstroms up to a few micrometers. These pores wit h nanometer and micrometer sizes greatly reduces the specific surface areas. In our Cambrian OM system, the ultra-micropore is the dominant pore type. The small differences of innate porous properties between the externally accessible regions and the total
CR IP T
regions indicate that most of pores in the organic structures are externally accessible. The density, porosity and specific surface area diverge slightly under different initial structure configurations, with the uncertainties less than 5 %, which shows a good reproductivity.
Properties of different organic models at 300 K and 0.1 MPa a
Porosity (%)
3
Density (g/cm )
VDW b
VDW c
He
VDW b
VDW c
AN
He
Specific surface area (m2 /g)
1.30 ± 0.03
10.95 ± 0.3
36.24 ± 1.6
36.36 ±1.3
1524 ± 64
4838 ± 203
4844 ± 205
OM e
1.28 ± 0.02
11.21 ± 0.3
36.62 ± 1.5
36.73 ±1.2
1497 ± 59
4789 ± 177
4812 ± 192
Kerogen+lighter fractions
1.27 ± 0.02
13.84 ± 0.3
36.93 ± 1.6
37.07 ± 1.4
1389 ± 55
4379 ± 182
4414 ± 189
Kerogen+heavier fractions
1.25 ± 0.01
17.78 ± 0.4
38.08 ± 1.7
38.28 ± 1.6
1423 ± 62
4221 ± 191
4251 ± 177
Kerogen in OM
1.22 ±0.01
21.21 ± 0.4
40.12 ± 1.8
40.32 ± 1.9
1641 ± 77
4217 ± 178
4236 ± 192
PT
a
M
OM d
ED
Structure
US
Table 10
Probe diameter of He molecule is 2.6 Å. VDW b refers to a probe diameter of zero defined over the externally accessible regions, while
CE
VDW c is equivalent to a probe diameter of zero defined over the whole pore regions. OM d refers to the OM structure relaxed under the annealing cycle (300~800~300 K), while OM d refers to the OM structure under the annealing cycle (300~900~300 K). The average deviations are determined on 10 structure configurations of organic models.
AC
To further analyze the effect of lighter fractions (methane, water, nitrogen and carbon dioxide) and heavier fractions (resin, fichtelite and toluene), the pore volume and specific surface area distributions of various organic models are compared in Fig. 8. Note that the pore size distribution of our organic models are narrow due to the finite sizes of the molecular models. This k ind of narrow pore size distribution of organic models has also been reported in previous work of molecular simulations ( Collell et al., 2014b; Ho et al., 2018). These pore size distributions are close to the experimental micropore size distributions on kerogen calculated from carbon dioxide adsorption (Rexer et al., 2014). Therefore, the pore network in our OM model can reasonably represent the ultra-microporous and limited microporous network of OM in shale formations. In this work, most of the pore volumes and specific surface areas are contributed by the ultra-micropores with pore width ranging between 1 and 3 Å. With the introduction of inclusions in the OM system, the pores larger than 3.3 Å decrease while the smaller pores increase. This means that
21
ACCEPTED MANUSCRIPT the introduced fractions enter these larger pores and divide them into inaccessible pores. It is interesting to note that compared with lighter fractions, the heavier fractions have a much smaller effect on the kerogen porous structure. The pore distributions of different organic models are illustrated in Fig. 9. The pore distribution of kerogen containing heavier frac tions (Fig. 9b) is fairly similar to that of innate kerogen structure (Fig. 9d). By contrast, a large number of pores in the kerogen structure are occupied by lighter fractions (Fig. 9c). The heavier fractions are more inclined to be part of the porous skeleton, while the lighter fractions tend to distribute in enterable pores attracted by heteroatom groups on pore surfaces. Moreover, these helium probed pores in the OM
CR IP T
system are observed to be poorly connected. We conclude previously that most of pores in the organic structures are externally accessible. Note that this conclusion applies to the probe with a diameter of zero, which means that throats at any size in t he OM skeleton are included to analyze the connectivity. However, when helium is used as probe molecules, the connectivity is analyzed based on the throats that have bigger sizes than the helium molecule. Therefore, we can conclude that the throats in the OM
US
skeleton are highly connected, but most of them have smaller sizes than the helium molecule. However, these smaller throats cannot be neglected because they may become accessible for gas transport when considering the swelling effect. (b) OM Kerogen+lighter frac. Kerogen+heavier frac. Kerogen in OM
120
M
-2
80
ED
60 40
0
1
2
3 4 Pore width, Å
PT
20 0
AN
160 -2
100
dV/dr, 10 Å
2
120
200
dS/dr, 10 Å
(a) OM Kerogen+lighter frac. Kerogen+heavier frac. Kerogen in OM
140
5
6
80 40 0
0
1
2
3 4 Pore width, Å
5
6
AC
CE
Fig. 8. Pore volume and specific surface area distribution. (a) Pore volume distribution. (b) Specific surface area distribution.
22
ACCEPTED MANUSCRIPT (b)
(c)
(d)
AN
US
CR IP T
(a)
Fig. 9. Pore distribution of various organic models probed by He. Skeleton atoms are wrapped by the purple pore surfaces. (a) The
M
OM model. (b) The kerogen+heavier fractions model. (c) The kerogen+lighter fractions model. (d) The kerogen model.
3.4 Gases adsorption behaviors in OM skeleton
ED
In this section, the adsorption characteristics of lighter gas species (methane, carbon dioxide and nitrogen) and heavier hydrocarbons species (ethane, propane and n-butane) in the porous OM skeleton were investigated by the GCMC simulations. As
PT
shown in Table 8, methane, carbon dioxide and nitrogen are dominate species based on the analytic data of gas chromatography on produced gas from Cambrian shale (Wu et al., 2014). We first check the adsorption behaviors of these dominated gas species in
CE
our OM skeleton. The absolute adsorption isotherms present a trend of type I isotherm according to the classification of Gregg and Sing (Gregg and Sing, 1982) (Fig.10a). The absolute adsorption amounts increase with pressure, and the adsorption capacity is carbon dioxide > methane > nitrogen. To compare with experimental adsorption results, the excess adsorption amounts are
AC
computed for these dominate gas species (Fig.10b). The experimental excess adsorption data for methane, carbon dioxide and nitrogen are performed on the high-mature kerogen sample in Barnett shale (Hu, 2014), Ruhr coal sample (Gensterblum et al., 2013) and Sulcis coal sample (Ottiger et al., 2008), respectively. Note that we aim to validate the reasonability of our OM model to represent the adsorption characteristics of porous OM skeleton. It is well known that the porous structures of OM skeletons in coal and shale are similar at molecular level when neglecting the effect of micro-fractures. Therefore, the experimental data on coal can be used to be compared with our simulated results.
It is well known that some minerals cannot be removed from the
experimental OM samples. The experimental data are normalized by the TOC so as to eliminate the effect of residual minerals before the comparison. The simulated excess adsorption isotherms of methane, carbon dioxide and nitrogen on our OM skeleton match fairly well with experimental results on authentic organic pores in both trend and magnitude. Note that it is extremely difficult or impossible to represent the whole range of pore sizes of OM under molecular and atomistic scale. As discussed in section 3.3.2, the OM skeleton in this work can represent the ultra-microporous and limited microporous network of OM in shale samples. Thus, this observation in Fig.10b indicates that ultra-micropores and limited micropores dominate the pore types for 23
ACCEPTED MANUSCRIPT adsorption of methane, carbon dioxide and nitrogen in realistic OM skeleton. Therefore, the generated OM skeleton in this work can be used to study the adsorption behaviors of methane, carbon dioxide and nitrogen. The simulated maximum excess adsorption amount of carbon dioxide (1.80 mmol/g) is fairly close to the reported one (1.84 mmol/g). The excess adsorption amount of car bon dioxide initially increases and then decreases with the increase of pressure. Similar with the carbon dioxide isotherm, the simulated methane and nitrogen excess adsorption isotherms are comparable with the experimental results. Compared with the documented experimental data, the simulated methane isotherm presents a higher adsorption rate, and a slight decrease when the press ure exceeds ~9 MPa can be observed. Notably, the later drop of excess adsorption amount is not observed for the nitrogen isotherm . The excess adsorption is closely correlated with the density difference between the adsorbed layer and the bulk. The excess adsorption amount reaches the maximum when there is a maximal density difference. The decrease in excess isotherms is because
CR IP T
the density difference between the adsorbed layer and the bulk reduces at high pressures. It can be concluded that the decrease in excess isotherms can be affected by the gas adsorption capacity, and a higher adsorption capacity corresponds to a lower pressure when the decrease happens. To investigate the thermodynamic property of gas adsorption on the OM skeleton, the isosteric heat of gases adsorption is presented in Fig.10c. The isosteric heat of gases adsorption reflects the combined effect of adsorbent-adsorbate interaction and adsorbate-adsorbate interaction. The isosteric heat of carbon dioxide adsorption initially decreases and then gently
US
increases with further rising pressure, while that of methane and nitrogen adsorption shows continually slight increase as the pressure rises. The early decrease of isosteric heat of carbon dioxide adsorption reflects the chemically heterogeneity of OM
AN
skeleton, since carbon dioxide preferentially fills these high-energy adsorption sites and then occupies these low-energy sites. The subsequent increase of isosteric heat of gases adsorption can be attributed to the increasing interaction between gas molecules. The limit isosteric heat, defined as the isosteric heat of gas adsorption at zero coverage (gas adsorption approaches zero), reflects
M
quantitatively the adsorption enthalpy between adsorbent and adsorbate without the influence from the adsorbate-adsorbate interaction. The limit isosteric heat on the OM skeleton is carbon dioxide (27.33 KJ/mol) > methane (20.76 KJ/mol)> nitrogen
ED
(17.51 KJ/mol), consistent with the order of gas adsorption capacity.
PT
Although there are nearly no heavier hydrocarbon species in the Cambrian shale gas (Table 8), the adsorption isotherms of ethane, propane and n-butane are simulated and compared with the experimental results in the work of Zhao et al. (Zhao et al., 2018) so as to check the possible effect of pore size or adsorbate size. As shown in Fig.10d-f, the simulation results are in poor
CE
agreement with the experimental results. The simulated adsorption amounts reach the plateaus at lower pressures, and the simulated maximum excess adsorption amounts are much lower than the measured results. These observations indicate that the
AC
adsorption behaviors of heavier hydrocarbon species cannot be represented in our OM skeleton with ultra-micropores and limited micropores, and the contribution of some bigger pores in realistic OM cannot be neglected. This conclusion shows an effect of adsorbate size on the dominated pore size for adsorption. This effect of adsorbate size can also be reflected by the order of absolute adsorption amounts of these heavier hydrocarbon species (ethane > propane > n-butane). The accessible pores in the OM skeleton decrease as the molecular sizes of these hydrocarbon species increase.
24
Excess adsorption amount, mmol/g
2.5 2.0 1.5 CH4
1.0
CO2 N2
0.5
0
4
8 12 Pressure, MPa
16
30
1.5
1.0 0.5
0.0 0
20
CO2
N2
20
4
8 12 Pressure, MPa
16
20
1.5 1.0 0.5 0.0
-0.5
0.0 0.0
ED
0.5 1.0 Pressure, MPa
1.5
Adsorption amount, mmol/g
CE
0.5
PT
1.0
CO2_exp.
N2_sim.
N2_exp.
16
20
(d)
0
4
8 Pressure, MPa
12
16
1.2
(e)
Absolute adsorption_sim. Excess adsorption_sim. Excess adsorption_exp.
AC
Adsorption amount, mmol/g
1.5
CH4_exp.
CO2_sim.
8 12 Pressure, MPa
M
0
CH4_sim.
Absolute adsorption_sim. Excess adsorption_sim. Excess adsorption_exp.
AN
CH4
4
2.0
(c)
25
15
(b)
Adsorption amount, mmol/g
0.0
2.0
CR IP T
(a)
US
3.0
Isosteric heat of adsorption, KJ/mol
Absolute adsorption amount, mmol/g
ACCEPTED MANUSCRIPT
0.9
(f)
Absolute adsorption_sim. Excess adsorption_sim. Excess adsorption_exp.
0.6
0.3
0.0 0.0
0.1
0.2 0.3 Pressure, MPa
0.4
0.5
Fig. 10. Adsorption characteristics of gases in OM skeleton. (a) Absolute adsorption amount of methane, carbon dioxide and nitrogen (338 K). (b) Excess adsorption amount of methane (338 K), carbon dioxide (318 K) and nitrogen (343 K). (c) Isosteric heat of gases adsorption (338 K). (d) Adsorption amount of ethane (333 K). (e) Adsorption amount of propane (333 K). (f) Adsorption amount of n-butane (333 K).
4 Conclusions In this work, we proposed a novel strategy for constructing hybrid shale OM model at molecular scale. A realistic molecular model of the Cambrian OM was constructed based on analytical structural and compositional results. Microstructure
25
ACCEPTED MANUSCRIPT characterization of the OM model and gases adsorption characteristics in the OM skeleton were further discussed. The main conclusions are summarized below. The Cambrian OM model is composed of kerogen, bitumen components and residual lighter components. The OM model shows a reasonable representation of the Cambrian OM in terms of structural parameters, generic compositions, physical density and porosity. Thus, it can be utilized to study the Cambrian OM at molecular scale. The steric OM structure is highly unordered and shows a remarkable folding of aromatic units. Lighter fractions tend to
CR IP T
appear in pairs around heteroatom groups, showing competitive adsorption behaviors. Instead of aggregating into clusters, wat er molecules are scattered around hydrophilic groups. Externally accessible ultra-micropore is the dominant pore type in the Cambrian OM system, while the pores accessible for gas molecules are poorly connected. The heavier bitumen fractions are inclined to be part of the porous skeleton, exerting a small effect on the porous structure. By contrast, the lighter fractions tend to
US
distribute in enterable pores due to the attraction of heteroatom groups, thus dividing these pores into massive inaccessible pores. The GCMC simulations show that the adsorption capacity in the OM skeleton is carbon dioxide > methane > nitrogen, as
AN
validated by the adsorption amount and limit isosteric heat. The adsorption behaviors of methane, carbon dioxide and nitrogen can be well represented by the OM skeleton. The gas adsorption capacity affects the corresponding pressure when the excess isotherm
M
reaches the maximum, and a higher adsorption capacity corresponds to a lower pressure at the maximum. There is a significant effect of adsorbate size on the adsorption behaviors of heavier hydrocarbon species (ethane, propane and n-butane). To represent
ED
the adsorption behaviors of heavier hydrocarbon species reasonably, the OM system should be augmented to create a wider pore
PT
size distribution.
The reported novel construction process for shale OM model demonstrates potential application for the representative
CE
structure reproduce of realistic shale OM under reservoir conditions at molecular level. The generated OM model can serve as a starting point for further theoretical investigations on the fundamental mechanisms governing the storage, transport, and rec overy
AC
of shale gas in heterogeneous organic nanopores.
Appendix. A. Supplementary data Supplementary data contain Table S1-S3 and Fig. S1-S3 associated with this manuscript.
Notes The authors declare no competing financial interest.
Acknowledgments This work was supported by the National Natural Science Foundation of China (Grant No. 51774298 and 51504265) and the Science Foundation for the Excellent Youth Scholars of China University of Petroleum (Beijing) (Grant No.2462015YQ0223). Computer time for this study was provided by the HP High Performance Computing Cluster of the State Key Laboratory of Heavy
26
ACCEPTED MANUSCRIPT Oil Processing at China University of Petroleum (Beijing).
References Bajc, S., Cvetkovic, O., Ambles, A., Vitorović, D., 2008. Characterization of type III kerogen from Tyrolean shale (Hahntennjoch, Austria) based on its oxidation products. J. Serb. Chem. Soc. 73, 463-478. Berendsen, H.J.C., Postma, J.P.M., van Gunsteren, W.F., DiNola, A., Haak, J.R., 1984. Molecular dynamics with coupling to an external bath. J. Chem. Phys. 81, 3684-3690. Bous ige, C., Ghimbeu, C.M., Vix-Guterl, C., Pomerantz, A.E., Suleimenova, A., Vaughan, G., Garbarino, G., Feygenson, M.,
Mater. 15, 576-582.
CR IP T
Wildgruber, C., Ulm, F. J., Pellenq, R.J.M., Coasne, B., 2016. Realistic molecular model of kerogen/'s nanostructure. Nat.
Budinova, T., Razvigorova, M., Tsyntsarski, B., Petrova, B., Ekinci, E., Yardim, M.F., 2009. Characterization of Bulgarian oi l shale kerogen revealed by oxidative degradation. Chem. Erde-Geochem. 69, 235-245.
Chalmers, G.R.L., Bustin, R.M., 2007. The organic matter distribution and methane capacity of the Lower Cretaceous strata of Northeastern British Columbia, Canada. Int. J. Coal Geol. 70, 223–239.
US
Chen, S., Zhu, Y., Wang, H., Liu, H., Wei, W., Fang, J., 2011. Shale gas reservoir characterization: A typical case in the southern Sichuan Basin of China. Energy 36, 6609-6616.
Chen, Y., Mastalerz, M., Schimmelmann, A., 2012. Characterization of chemical functional groups in macerals across different
AN
coal ranks via micro-FTIR spectroscopy. Int. J. Coal Geol. 104, 22-33.
Collell, J., Galliero, G., Gouth, F., Montel, F., Pujol, M., Ungerer, P., Yiannourakou, M., 2014a. Molecular simulation and
conditions. Micropor. Mesopor. Mat. 197, 194-203.
M
modelisation of methane/ethane mixtures adsorption onto a microporous molecular model of kerogen under typical reservoir
ED
Collell, J., Ungerer, P., Galliero, G., Yiannourakou, M., Montel, F., Pujol, M., 2014b. Molecular simulation of bulk organic matter in Type II shales in the middle of the oil formation window. Energ. Fuel. 28, 7457-7466.
PT
Connolly, M.L., 1983. Solvent-accessible surfaces of proteins and nucleic acids. Science 221, 709-713. Curtis J.B., 2002. Fractured shale-gas systems. Am. Assoc. Pet. Geol. Bull. 86, 1921-1938. Feng, F., Akkutlu, I.Y., 2015. Flow of hydrocarbons in nanocapillary: A non-equilibrium molecular dynamics study. SPE Asia
CE
Pacific Unconventional Resources Conference and Exhibition, Brisbane, Australia, 9-11 November, SPE-177005. Gensterblum, Y., Merkel, A., Busch, A., Krooss, B.M., 2013. High-pressure CH4 and CO2 sorption isotherms as a function of coal
AC
maturity and the influence of moisture. Int. J. Coal Geol. 118, 45–57.
Gregg, S.J., Sing, K.S., 1982. Adsorption, surface area and porosity. Academic Press, New York, p.4. Ho, T.A., Wang, Y., Criscenti, L.J., 2018. Chemo-mechanical coupling in kerogen gas adsorption/desorption. Phys. Chem. Chem. Phys. 20, 12390-12395. Hoover, W.G., 1985. Canonical dynamics: Equilibrium phase-space distribution. Phys Rev A Gen Phys 31, 1695-1697. Hu, H., 2014. Methane adsorption comparison of different thermal maturity kerogens in shale gas system. Chin. J. Geochem. 33, 425-30.
Huang, L., Ning, Z.F., Wang, Q., Ye, H.T., Qin, H.B., 2017a. Molecular Simulation of CO2 sequestration and enhanced gas recovery in gas rich shale: An insight based on realistic kerogen model. SPE Abu Dhabi International Petroleum Exhibition & Conference, Abu Dhabi, UAE, 13-16 November, SPE-188216. Huang, L., Ning, Z.F., Wang, Q., Qi, R.R., Li, J., Zeng, Y., Ye, H.T., Qin, H.B., 2017b. Thermodynamic and structural characterization of bulk organic matter in Chinese Silurian Shale: Experimental and molecular modeling studies. Energ. Fuel. 31, 4851-4865.
27
ACCEPTED MANUSCRIPT Huang, L., Ning, Z.F., Wang, Q., Zhang, W.T., Cheng, Z.L., Wu, X.J., Qin, H.B., 2018a. Effect of organic type and moisture on CO2 /CH4 competitive adsorption in kerogen with implications for CO 2 sequestration and enhanced CH4 recovery, Appl. Energy 210, 28–43. Huang, L., Ning, Z.F., Wang, Q., Qi, R.R., Zeng, Y., Qin, H. B., Ye, H.T., Zhang, W.T., 2018b. Molecular simulation of adsorption behaviors of methane, carbon dioxide and their mixtures on kerogen: Effect of kerogen maturity and moisture content. Fuel 211, 159-172. Kazemi, M., Takbiri-Borujeni, A., 2016. Molecular dynamics study of carbon dioxide storage in carbon-based organic nanopores. SPE Annual Technical Conference and Exhibition, Dubai, UAE, 26-28 September, SPE-181705. Kelemen, S.R., Afeworki, M., Gorbaty, M.L., Sansone, M., Kwiatek, P.J., Walters, C.C., Freund, H., Siskin, M., 2007. Direct
CR IP T
characterization of kerogen by X-ray and solid-state 13 C nuclear magnetic resonance methods. Energ. Fuel. 21, 1548-1561. Kim, T.H., Cho, J., Lee, K.S., 2017. Evaluation of CO 2 injection in shale gas reservoirs with multi-component transport and geomechanical effects. Appl. Energy 190, 1195-1206.
Kurniawan, Y., Bhatia, S.K., Rudolph, V., 2006. Simulation of binary mixture adsorption of methane and CO 2 at supercritical conditions in carbons. AIChE J. 52, 957-967.
US
Li, J., Li, X.F., Wang, X.Z., Li, Y.Y., Wu, K.L., Shi, J.T., Yang, L., Feng, D., Zhang, T., Yu, P.L., 2016a. Water distribution characteristic and effect on methane adsorption capacity in shale clay. Int. J. Coal Geol. 159, 135-154.
AN
Li, J., Li X.F., Wu, K.L., Wang, X.Z., Shi, J.T., Yang, L., Zhang, H., Sun, Z., Wang, R., Feng, D., 2016b. Water sorption and distribution characteristics in shale clay: Effect of surface force. Energ. Fuel. 30, 8863−8874. Li, J., Li, X.F., Wu, K.L., Feng, D., Zhang, T., Zhang, Y.F., 2017. Thickness and stability of water film confined inside nanoslits
M
and nanocapillaries of shale and clay. Int. J. Coal Geol. 179, 253-268.
Li, Q., Han, X., Liu, Q. Jiang, X., 2014. Thermal decomposition of Huadian oil shale. Part 1. Critical organic intermediates, Fuel
ED
121, 109-116.
Li, S., Yue, C., 2003. Study of pyrolysis kinetics of oil shale. Fuel 82, 337-342.
PT
Liu, S., Deng, B., Zhong, Y., Ran, B., Yong, Z., Sun, W., Yang, D., Jiang, L., Ye Y., 2016a. Unique geological features of burial and superimposition of the Lower Paleozoic shale gas across the Sichuan Basin and its periphery. Earth Sci. Front. 23, 11-28.
CE
Liu, X.Q., He, X., Qiu, N.X., Yang, X., Tian, Z.Y., Li, M.J., Xue, Y., 2016b. Molecular simulation of CH4 , CO2, H2O and N2 molecules adsorption on heterogeneous surface models of coal. Appl. Surf. Sci. 389, 894-905. Loucks, R.G., Ruppel, S.C., 2007. Mississippian Barnett Shale: Lithofacies and depositional setting of a deep-water shale-gas
AC
succession in the Fort Worth Basin, Texas. Am. Assoc. Pet. Geol. Bull. 91, 579−601. Loucks, R.G., Reed, R.M., Ruppel, S.C., Jarvie, D.M., 2009. Morphology, genesis and distribution of nanometer -scale pores in siliceous mudstones of the Mississippian Barnett Shale. J. Sediment. Res. 79, 848-861. Lu, X.Q., Jin, D.L., Wei, S.X., Zhang, M.M., Zhu, Q., Shi, X.F., Deng, Z.G., Guo, W.Y., Shen, W.Z., 2015. Competitive adsorpt ion of a binary CO2 -CH4 mixture in nanoporous carbons: effects of edge-functionalization. Nanoscale 7, 1002-1012. Middleton, R.S., Carey, J.W., Currier, R.P., Hyman, J.D., Kang, Q., Karra, S., Martínez, J.J., Porter, M.L., Viswanathan, H.S ., 2015. Shale gas and non-aqueous fracturing fluids: Opportunities and challenges for supercritical CO 2 . Appl. Energy 147, 500-509. Ojima, J., 2003. Determining of crystalline silica in respirable dust samples by infrared spectrophotometry in the presence o f interferences. J. Occup. Health 45, 94-103. Okiongbo, K.S., Aplin, A.C., Larter, S.R., 2005. Changes in type II kerogen density as a function of maturity: Evidence from the Kimmeridge clay formation. Energ. Fuel. 19, 2495−2499. Orendt, A., Pimienta, I., Badu, S., 2013. Three-dimensional structure of the Siskin Green River Oil Shale Kerogen Model: A 28
ACCEPTED MANUSCRIPT comparison between calculated and observed properties. Energ. Fuel. 27, 702−710. Ottiger, S., Pini, R., Storti, G., Mazzotti, M., 2008. Measuring and modeling the competitive adsorption of CO 2 , CH4 , and N2 on a dry coal. Langmuir 24, 9531-9540. Paylor, A., 2015. The social–economic impact of shale gas extraction: A global perspective. Third World Q. 38, 340-355. Rexer, T.F., Mathia, E.J., Aplin, A.C., Thomas, K.M., 2014. High-pressure methane adsorption and characterization of pores in Posidonia shales and isolated kerogens. Energ. Fuel. 28, 2886-2901. Ross, D.J.K., Bustin, R.M., 2009. The importance of shale composition and pore structure upon gas storage potential of shale gas reservoirs. Mar. Pet. Geol. 26, 916-927. Ru, X., Cheng, Z., Song, L., Wang, H., Li, J., 2012. Experimental and computational studies on the average molecular structur e of
CR IP T
Chinese Huadian oil shale kerogen. J. Mol. Struct. 1030, 10-18.
Shinn, J.H., 2013. From coal to single-stage and two-stage products: A reactive model of coal structure. Fuel 63, 1187-1196. Sun, H., Mumby, S.J., Maple, J.R., Hagler, A.T., 1994. An ab initio CFF93 all-atom forcefield for polycarbonates. J. Am. Chem. Soc. 116, 2978-2987.
Sun, H., 1995. Ab initio calculations and forcefield development for computer simulation of polysilanes. Macromolecules 28,
US
701-712.
Sun, H.Y., Zhao, H., Qi, N., Li, Y., 2017. Molecular insights into the enhanced shale gas recovery by carbon dioxide in kerogen slit
AN
nanopores. J. Phys. Chem. C 121, 10233-10241.
Tissot, B.P., Welte, D.H., 1984. Petroleum Formation and Occurrence. Springer: Berlin, Germany, p720. Tong, J., Han, X., Wang S., 2011. Evaluation of structural characteristics of Huadian oil shale kerogen using direct techniques
M
(solid-state 13 C NMR, XPS, FTIR, and XRD). Energ. Fuel. 25, 4006-4013.
molecular modeling. Fuel 181, 330-339.
ED
Tong, J., Jiang, X., Han, X., Wang, X., 2016. Evaluation of the macromolecular structure of Huadian oil shale kerogen using
spectroscopy. Fuel 65, 541-546.
13
C nmr
PT
Trewhella, M.J., Poplett, I.J.F., Grint, A., 1986. Structure of Green River oil shale kerogen: Determination using solid state
Ungerer, P., Collell, J., Yiannourakou, M., 2014. Molecular modeling of the volumetric and thermodynamic properties of kerogen:
CE
Influence of organic type and maturity. Energ. Fuel. 29, 91-105. Vandenbroucke, M., 2003. Kerogen: from types to models of chemical structure. Oil Gas Sci. Technol. 58, 243−269. Vandenbroucke, M., Largeau, C., 2007. Kerogen origin, evolution and structure. Org. Geochem. 38, 719-833.
AC
Wang, J., Du, J., Chang, L., Xie, K., 2010. Study on the structure and pyrolysis characteristics of Chinese western coals. Fuel Process Technol. 91, 430-433.
Wang, X.Q., Zhai, Z.Q., Jin, X., Wu, S.T., Li, J.M., Sun, L., Liu, X.D., 2016. Molecular simulation of CO 2 /CH4 competitive adsorption in organic matter pores in shale under certain geological conditions. Petrol. Explor. Dev. 43, 841-848. Wei, Z., Gao, X., Zhang, D., Da, J., 2005. Assessment of thermal evolution of kerogen geopolymers with their structural parameters measured by solid-state 13 C NMR spectroscopy. Energ. Fuel. 19, 240-250. Weijermars, R., 2013. Economic appraisal of shale gas plays in Continental Europe. Appl. Energy 106, 100–115. Weijermars, R., 2014. US shale gas production outlook based on well roll-out rate scenarios. Appl. Energy 124, 283–297. Wu, K.L., Li, X.F., Wang, C.C., Yu, W., Chen, Z.X., 2015. Model for surface diffusion of adsorbed gas in nanopores of shale g as reservoirs. Ind. Eng. Chem. Res. 54, 3225-3236. Wu, K.L., Chen, Z.X., Li, X.F., Dong, X.H., 2016a. Methane storage in nanoporous material at supercritical temperature over a wide range of pressures. SCI REP-UK 6, 33461. 29
ACCEPTED MANUSCRIPT Wu, K.L., Li, X.F., Guo, C.H., Wang, C.C., Chen, Z.X., 2016b. A unified model for gas transfer in nanopores of shale-gas reservoirs: coupling pore diffusion and surface diffusion. SPE J. 21, 1583-1611. Wu, W., Huang, S., Hu, G., Gong, D., 2014. A comparison between shale gas and conventional gas on geochemical characteristics in Weiyuan area. Nat. Gas Geosci. 25, 1994-2002. Yang, F., Ning, Z.F., Wang, Q., Liao, H.M., 2013. Integrated study of reservoir characteristics of a shale gas reservoir: A c ase study from Sichuan Basin of China. SPE Asia Pacific Oil and Gas Conference and Exhibition, Jakarta, Indonesia, 22-24 October, SPE-165870. Yang, F., Ning, Z.F., Wang, Q., Liu, H.Q., 2016. Pore structure of Cambrian shales from the Sichuan Basin in China and implications to gas storage. Mar. Petrol. Geol. 70, 14-26.
CR IP T
Yen, T., 1979. Structural difference between petroleum and coal-derived asphaltenes. Am. Chem. Soc., Div. Pet. Chem., Prepr. 24, 901-909.
Yiannourakou, M., Ungerer, P., Leblanc, B., Rozanska, X., Saxe, P., Vidal-Gilbert, S., Gouth, F., Montel, F., 2013. Molecular simulation of adsorption in microporous materials. Oil Gas Sci. Technol. 68, 977−994.
Energy 148, 49–65.
US
Yuan, J., Luo, D., Feng, L., 2015a. A review of the technical and economic evaluation techniques for shale gas development. Appl.
Yuan, Q., Zhu, X., Lin, K., Zhao, Y.P., 2015b. Molecular dynamics simulations of the enhanced recovery of confined methane wit h
AN
carbon dioxide, Phys. Chem. Chem. Phys. 17, 31887−31893.
Zeng, Y., Wu, C., 2007. Raman and infrared spectroscopic study of kerogen treated at elevated temperatures and pressures. Fuel 86, 1192-1200.
M
Zhang, J., Liu, K., Clennell, M. B., Dewhurst, D.N., Pervukhina, M., 2015. Molecular simulation of CO 2 –CH4 competitive adsorption and induced coal swelling. Fuel 160, 309-317.
ED
Zhang, Q., Han, K., Li, S., Li, M., Li, J., Ren, K., 2017. Synthesis of garlic skin-derived 3D hierarchical porous carbon for high-performance supercapacitors. Nanoscale 10, 2427.
PT
Zhao, H., Wu, T., Firoozabadi, A., 2018. High pressure sorption of various hydrocarbons and carbon dioxide in Kimmeridge Blackstone and isolated kerogen. Fuel 224, 412-423
pores. Fuel 180, 718-726.
CE
Zhou, B., Xu, R., Jiang, P., 2016. Novel molecular simulation process design of adsorption in realistic shale kerogen spheric al
Zhu, X., Zhao, Y.P., 2014. Atomic mechanisms and equation of state of methane adsorption in carbon nanopores. J. Phys. Chem. C
AC
118, 17737-17744.
Graphical abstract Highlights: (1) Novel construction process for OM model is elaborated. (2) Direct characterization experiments are used to study kerogen chemical structure. (3) Representative molecular model of Cambrian OM is generated.
30
ACCEPTED MANUSCRIPT (4) Microstructure characterization is performed on the OM model.
AC
CE
PT
ED
M
AN
US
CR IP T
(5) Adsorption behaviors of gas components in the OM skeleton are studied.
31
Liang Huang, Zhengfu Ning, Qing Wang, Hongtao Ye, Zizheng Wang, Zheng Sun, Huibo Qin PII: DOI: Reference:
S0166-5162(18)30551-2 doi:10.1016/j.coal.2018.09.001 COGEL 3075
To appear in:
International Journal of Coal Geology
Received date: Revised date: Accepted date:
7 June 2018 31 August 2018 1 September 2018
Please cite this article as: Liang Huang, Zhengfu Ning, Qing Wang, Hongtao Ye, Zizheng Wang, Zheng Sun, Huibo Qin , Microstructure and adsorption properties of organic matter in Chinese Cambrian gas shale: Experimental characterization, molecular modeling and molecular simulation. Cogel (2018), doi:10.1016/j.coal.2018.09.001
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
ACCEPTED MANUSCRIPT Microstructure and adsorption properties of organic matter in Chinese Cambrian gas shale: Experimental characterization, molecular modeling and molecular simulation Liang Huang a,b , Zhengfu Ning a,c* , Qing Wang a
a,c
, Hongtao Ye
a,c
, Zizheng Wang c, Zheng Sun a,c, Huibo Qin d
State Key Laboratory of Petroleum Resources and Prospecting, China University of Petroleum (Beijing), Beijing
b
CR IP T
102249, China. Department of Chemical and Biomolecular Engineering, University of California, Berkeley, Berk eley, CA
94720-1462, United States.
Department of Petroleum Engineering, China University of Petroleum (Beijing), Beijing 102249, China.
d
State Key Laboratory of Heavy Oil Processing, China University of Petroleum (Beijing), Beijing, 102249, China.
AN
*Corresponding author: e-mail address: [email protected] .
US
c
ABSTRACT
M
A representative molecular model of shale organic matter (OM) is prerequisite to further theoretic investigation on the
ED
fundamental mechanisms governing storage, transport and recovery of shale gas. In this work, a systematic strategy to prepare structural and compositional properties of OM is reported first, and then a realistic molecular model of Chinese Cambrian OM is
PT
generated based on analytical experimental data. Microstructure characterization and adsorption simulation are further performed using molecular dynamics simulation and grand canonical Monte Carlo simulations, respectively. The OM model, composed of
CE
kerogen macromolecules, bitumen components and residual lighter components, shows a reasonable representation of realistic Cambrian OM with respect to structural parameters, generic compositions, physical density and porosity. The OM porous network
AC
consists of dominant ultra-micropores and limited micropores. Compared with heavier components, lighter components are more inclined to occupy accessible pores. Interestingly, lighter components are observed to appear in pairs due to competitive adsorption around heteroatom groups. Water molecules are scattered in the system because of abundant functional groups and poor pore connectivity. The OM skeleton represents the adsorption behaviors of methane, carbon dioxide and nitrogen well. The adsorption capacity is carbon dioxide > methane > nitrogen. A higher adsorption capacity corresponds to a lower pressure when the excess isotherm reaches the maximum. The adsorption behaviors of heavier hydrocarbon species (ethane, propane and n-butane) cannot be represented in the OM skeleton with ultra-micropores and limited micropores, and the effect of molecular sizes of these species cannot be neglected. This work reports a systematic construction process for realistic molecular model of shale OM, and the representative OM model can serve as a starting point to explore gas adsorption and transport mechanism in shale organic pores at
1
ACCEPTED MANUSCRIPT microscopic scale. Keywords: Organic matter; molecular model; kerogen; experimental characterization; shale; adsorption
1 Introduction Shale gas has been considered as one of the most important substitutes for conventional hydrocarbon resources due to its considerable resource abundance and high utilization efficiency (Weijermars, 2014; Yuan et al., 2015a; Paylor, 2017). However, the recovery of most shale gas reservoirs remains below 30 % because of complex gas transport mechanisms and tight reservoir
CR IP T
characteristics, as well as depleted development regime (Weijermars, 2013; Kim et al., 2017). With the rising demand for clean fossil fuels, potential technologies for enhanced shale gas recovery like carbon dioxide sequestration have attracted enormous attentions worldwide (Middleton et al., 2015; Huang et al., 2017a, 2018a, 2018b).
Shale gas in subsurface rock can be divided into adsorbed gas on various pore surfaces, free gas in fractures and pores and dissolved gas in organic matter (OM). (Wu et al., 2016a). A large proportion of shale gas is located in nanoporous nodules as a
US
form of adsorbed gas, which, for example, can take up 60-85 % in Lewis Shale (Curtis, 2002). Adsorbed gas exploitation is the key
AN
issue of enhanced shale gas recovery (Wu et al., 2015, 2016b). Adsorbed gas is traditionally believed to be distributed in organic matter and clay mineral (Curtis, 2002). However, recent investigations indicate that the contribution of clay minerals to adsorbed
M
gas content is much limited because of their hydrophilic nature, which greatly hinders the gas adsorption (Li et al., 2016a, 2016b, 2017). Accordingly, most of adsorbed gas in shale is stored in the OM nanopores, which have large specific surface areas and
ED
strong adsorption potentials (Ross and Bustin, 2009). A better understanding of gas adsorption and transport behaviors in shale OM nanopores can help improve the performance of enhanced shale gas recovery technologies. Compared with experimental
PT
measurements, molecular simulations including molecular dynamics (MD) and grand canonical Monte Carlo (GCMC) simulations
CE
are more powerful tools to gain insights into the microbehaviors of fluid mixtures in complex systems at high pressure and temperature conditions (Zhang et al., 2015). To further study the microscopic mechanisms of enhanced shale gas recovery with
AC
molecular simulation method, it is of fundamental significance to c haracterize the structure of shale OM and construct representative molecular model.
To date, various organic models related to shale OM fractions have been reported. However, most of these organic models are simplified carbon nanopore models such as carbon nanotubes (Zhu and Zhao, 2014; Yuan et al., 2015b) and graphene slits (Kurniawan et al., 2006; Kazemi and Takbiri-Borujen, 2016). Some attempts have been performed to improve the complexity of these simplified organic models. Wang et al. (Wang et al., 2016) constructed a composite pillar-layer model to represent the OM pore structure in Chinese Lower Silurian Longmaxi shale based on experimental characterization. Feng and Akkutlu (Feng and Akkutlu, 2015) built an improved carbon nanotube model with roughness by setting bumps and trenches in the nanotube surfaces. These simplified models fail to capture the heterogeneity of OM structure. Some investigations have also been conducted to consider the surface heterogeneity of organic models. Lu et al. (Lu et al., 2015) and Sun et al. (Sun et al., 2017) added various 2
ACCEPTED MANUSCRIPT heteroatom functional groups in the edge of carbon clusters to represent the chemical surface heterogeneity. Liu et al. (Liu et al., 2016b) generated physically heterogeneous surface model by adsorbing carbon adatoms on graphene. In addition to surface heterogeneity, realistic shale OM also has heterogeneous porous structure. Ungerer et al. (Ungerer et al., 2014) developed realistic kerogen molecular models with different organic types and thermal maturities based on documented experimental data by Kelemen and coworkers (Kelemen et al., 2007). The heterogeneous chemical structures of these kerogen models are constrained by analytical results, while the porous structures are generated by the phase equilibria using the MD simulations. To optimize the
CR IP T
porous network, Collell et al. (Collell et al. 2014a) inserted dummy particles into a mature type II kerogen matrix model, generating heterogeneous porous network with both artificial pores and kerogen matrix pores. Zhou et al. (Zhou et al., 2016) inserted cutter clusters instead of dummy particles into kerogen matrix models to further improve the pore shape of artificial pores. However, it should be noted that authentic shale OM is one kind of mixture of multi-components, containing bitumen components
US
and residual lighter components, as well as kerogen macromolecules. These kerogen models fail to consider the effect of OM fragments on the bulk OM skeleton structure. To our best knowledge, except in our recent work which reported the construction of
AN
a realistic OM molecular model of Chinese Silurian shale (Huang et al., 2017b), there is only one published report on bulk molecular model of shale OM. Collell et al. (Collell et al., 2014b) constructed a bulk type II shale OM model in the middle of the
M
oil formation window constrained by reported experimental data. They adopted the assumption that hydrocarbons produced in the maturation process were unexpelled, which cannot be applied to OM in gas shale with high maturity, because produced
ED
hydrocarbons are inclined to expel during the maturation process. Also, inclusion of bitumen components in their OM model was
PT
not constrained by analytical data. In our recent work (Huang et al., 2017b), we generated a bulk OM molecular model representative of Chinese Silurian shale. Nevertheless, most of structural parameters of the OM skeleton were derived from
structural parameters.
CE
experiment-based correlations rather than direct characterization results, which should be further improved with detailed analytical
AC
Laboratory experiments serve as fundamental approaches to characterize the compositional and structural properties of OM components, which are essential for the generation of molecular models of OM model (Vandenbroucke and Largeau, 2007; Bous ige et al., 2016). So far, various experimental methods including chemical approaches and physical approaches have been performed to characterize the structural parameters of OM components. Chemical approaches such as chemical degradation (Bajc et al., 2008; Budinova et al., 2009), thermal analys is (Li and Yue, 2003; Wang et al., 2010) and pyrolysis–gas chromatography– mass spectroscopy (Ru et al., 2012) require the decomposition of experimental sample. By contrast, physical approaches such as Fourier transform infrared (FTIR) spectroscopy (Zeng and Wu, 2007), Raman spectroscopy (Zeng and Wu, 2007; Zhang et al., 2017), solid-state
13
13
C nuclear magnetic resonance ( C NMR) spectroscopy (Wei et al., 2005; Tong et al., 2016) and X-ray
photoelectron spectroscopy (XPS) (Kelemen et al., 2007; Zhang et al., 2017) are nondestructive direct methods for structural characterization. Compared with destructive chemical approaches, direct physical characterization experiments can provide more 3
ACCEPTED MANUSCRIPT preserving structural information and less associated uncertainties. However, only averaged structural information can be obtained through direct characterization, and specific components and elementary compositions still need to be quantified by destructive methods (Kelemen et al., 2007). Therefore, to develop representative molecular models of complex organic samples, it is better to integrate experimental results from both indirect and direct characterization methods. In this work, the OM of Cambrian shale in Sichuan Basin, currently the most significant target zone for shale gas exploration and exploitation in China, was targeted for analysis. Based on our systematic construction strategy, we aimed to generate a realistic
CR IP T
bulk OM model containing kerogen macromolecules, bitumen components and residual lighter components by a combination of experimental characterization and molecular simulation. The chemical structure of Cambrian kerogen was characterized by elementary analysis, FTIR spectroscopy,
13
C NMR spectroscopy and XPS technique. The bitumen components were constrained
by the experimental group separation results. The residual lighter components resulted from the atomic balance assumption,
US
bitumen group separation and gas chromatography of produced gas. These OM components were integrated to generate the bulk OM molecular model using the molecular simulation method. Microstructure characterizations including steric structure and
AN
porous structure of the OM model were then conducted with the MD simulations. The gases adsorption characteristics including absolute and excess adsorption amount and isosteric heat of gas adsorption in the OM skeleton were further explored with the
M
GCMC simulations. This work provides a systematic construction process of bulk OM molecular model, which can facilitate the representative reproduce of realistic OM structure at molecular level. The proposed OM molecular model, along with its
ED
adsorption characteristics, is expected to serve as a starting point for further theoretical investigations on gas adsorption and
2 Materials and methodology
CE
2.1 Materials
PT
transport mechanism in organic pores at microscopic scale.
Two black shale outcrops were sampled from the Lower Cambrian Niuititang Formation in Sichuan Basin, Southwest China.
AC
The two samples were collected from different sampling sites located in the same district. The Lower Cambrian Niuititang Formation is deposited in a deep-water shelf environment. Pelagic algae, fungi and zooplankton are the main sources of Cambrian OM (Yang et al., 2016). A typical reservoir condition (381 K and 33.6 MPa) can be determined for the Cambrian OM according to the reported temperature and pressure gradients (Wu et al., 2014; Liu et al., 2016a; Wang et al., 2016). Detailed information on the geology, stratigraphy and resource potential for Cambrian shale has been previously reported (Yang et al., 2013).
2.2 Characterization experiments Large chunks of shale samples were crushed first, and then grounded by a ball mill. The milled powder samples were dry sieved by a 200 mesh sieve to collect particles smaller than 75 µm for subsequent measurement. A LECO CS230 carbon/sulfur analyzer was utilized to measure the total organic carbon (TOC) of investigated samples. The milled shale powders were treated with HCl to dissolve carbonates first, and then pyrolyzed up to 813.15 K. Vitrinite-like macerals 4
ACCEPTED MANUSCRIPT were targeted for thermal maturity evaluation because of the absence of vitrinite in shale of Niutitang Formation. We used a MPV-SP microscope equipped with oil-immersion objective lens and photometer to measure the reflectance of vitrinite-like materials. A reported mathematical conversion equation was then utilized to convert the measured reflectance into equivalent vitrinite reflectance (Chen et al., 2011). Kerogen samples analyzed below were isolated from the shale powders by a combination of acid treatment and Soxhlet extraction process. The shale powders were demineralized by a succession of HCL, and HF treatments. Carbonates were removed
CR IP T
by 20 wt% HCl acid at 343.15 K. Silicates were removed by a mixture of HCl and HF acids (20 wt% HCl + 40 wt% HF) at 343.15 K. Note that pyrites may not be removed from the samples, since there was no HNO 3 treatment in the acid treatment process. Therefore, it is important to check the possible effect of pyrites on kerogen characterization. The remains after each acid treatment were filtered and washed with deionized water and oven-dried at 353.15 K. Thereafter, Soxhlet extraction with toluene was utilized
US
to extract soluble bitumen components. The extracted bitumen components were then treated with chromatography column for group separation. The obtained kerogen samples were oven-dried under vacuum first, and then treated by a Thermo Scientific
In this work, FTIR spectroscopy,
13
AN
FLASH 2000 CHNS-O Analyzer for elemental analysis.
C NMR spectroscopy and XPS techniques were utilized to characterize the chemical
M
structure of the studied kerogen samples. The FTIR spectra were measured by a MAGNA-IR56 FTIR spectrometer. Sample powders of 0.5 mg were mixed with 200 mg dry KBr powders, then pressed under vacuum. FTIR spectra measurement was
13
C NMR spectra were measured at 100.64 MHz by a Bruker AVANCEIII TM-600 MHz spectrometer with the MAS magic
PT
The
ED
performed at ambient temperature, and spectra were recorded between 4000 and 400 cm -1 with 50 scans at a resolution of 4 cm -1.
angle spinning and using a ZrO2 rotor (4 mm diameter). A recycle delay time of 3 s and a contact time of 4 ms were adopted in the
CE
cross-polarization process. The XPS measurement on kerogen powders was conducted on a Thermo Scientific ESCALAB 250Xi 2
spectrometer using monochromatic Al Kα radiation (1486.6 ev). Ar+ scanning gun was used on an area of 8×8 mm . The sputtering
AC
rate was 4 nm/min, and the transmit current was 10 mA. To specify and quantify the functional groups and elemental forms, Lorentz peak-fitting was performed on the measured spectra curves.
2.3 Molecular modeling method In this section, the construction process of average molecular model for Cambrian OM is discussed. The OM model is integrated from individual fractions, including kerogen, bitumen fractions and lighter components. Kerogen. Kerogen has complex macromolecule structure. It is composed of polyaromatic units with cross-linkages such as alkyl, ether and sulfur bridge. Kerogen constitutes the basic skeleton of OM, and can represent most of the physical and chemical 13
characteristics of OM. In this work, we adopt a combination of elementary analysis, FTIR, C NMR and XPS techniques to obtain detailed analytical data of elemental compositions, structural parameters and functional groups for Cambrian kerogen. These analytic data are utilized to derive the kerogen chemical structure. The building procedure of the kerogen model follows the 5
ACCEPTED MANUSCRIPT process in our previous work (Huang et al., 2017b). In this study, carbon structural parameters of Cambrian kerogen are derived from the
13
C NMR spectrum. Only organic
carbon forms including aromatic carbons and aliphatic carbons are considered. The functional group forms of O, N and S are analyzed based on the XPS results. For oxygenated groups, we exclude these components in residual minerals and water. In particular, carbonyl groups, alcoholic groups and ether groups are considered. For nitrogen-bearing groups, pyridine, pyrrole, amine, as well as nitrogen oxide are considered. In the case of sulfur-bearing groups, we not only focus on the organic sulfur
CR IP T
groups including thiophene, sulfoxide and sulfoether, but also check the possible influence of pyrite minerals. The FTIR spectra are utilized to qualitatively characterize the functional groups of Cambrian kerogen, which can serve as a self-consistent check 13
with the results of C NMR and XPS techniques.
The kerogen structure topology is largely inspired from previous kerogen fractions in the literature (Shinn, 1984). A unit of
US
100-260 carbon atoms for each kerogen molecule has been reported as the minimized size requirement to capture the major structural features of kerogen (Orendt et al., 2013; Ungerer et al., 2014). A total number of 210 carbon atoms for single Cambrian
AN
kerogen molecule is adopted to meet the size requirement. In addition to linear alkyl fragments, branched alkyl, oxygen-, sulfurand nitrogen-bridge segments are also introduced as cross-links between aromatic units. Moreover, apart from aromatic ring units,
M
monocyclic and polycyclic saturated rings are also added to match the H/C ratio of Cambrian kerogen. The initial version of kerogen molecule is designed to be a flat-shape structure with minimized extension so as to approach the equilibrium conformation.
ED
Thereafter, adjustments on the functional group locations and polyaromatic unit shapes are performed to improve the consistency
PT
between the simulated and experimental 13 C NMR spectra. The pcff force field is adopted throughout this work. We adopt the pcff force field to assign atomic types and charges for kerogen atoms so as to prepare the kerogen unit for further structure relaxation.
CE
The reason that the pcff force field is used will be discussed later. Bitumen. Bitumen is defined as the component of OM that can be soluble in dichloromethane (Vandenbroucke, 2003). It can
AC
be divided into asphaltene, resin, aromatic hydrocarbons, and saturated hydrocarbons. The asphaltene and resin molecules are composed of polyaromatic units with polar heteroatoms. Their structures are similar to that of kerogen, but with a lower molecular weight. Asphaltene is soluble in polar solvent but insoluble in n-alkanes, while resin is soluble in liquid n-alkanes. Compared with asphaltene, resin has a lower aromaticity with more polar fractions and alkane tails. The structure sizes of asphaltene and resin are determined to ensure an economic and reasonable simulation system, as reported in our latest work (Huang et al., 2017b). To be consistent with the conclusion of Yen (Yen, 1979), the nitrogen group is represented by quinoline unit, while the oxygen group is in the form of ether unit. In the case of aromatic hydrocarbons and saturated hydrocarbons, a limited number of representative components are targeted in this study. Toluene is chosen as the representative aromatic component, as adopted in Collell’s work (Collell et al., 2014b). While fichtelite, a common tricyclic alkane in bitumen, is chosen as the representative saturated component, since most of long 6
ACCEPTED MANUSCRIPT aliphatic components disappear in the late catagenesis stage. Also, some methane molecules are included in the saturated components to represent residual lighter alkyl components. The specific molecular sizes and quantities of representative components are adjusted to match the experimental bitumen group separation results in the Cambrian OM. Lighter components. In this work, methane, ethane, carbon dioxide, nitrogen and water are defined as lighter components. These lighter components are assumed to be generated from kerogen during the maturation evolution. By adopting the atomic balance assumption (Loucks and Ruppel, 2007), namely there is no atomic loss during the maturation process, we determine the
CR IP T
initial total lighter components under reservoir conditions. Detailed analytical data of the immature type II kerogen in Duvernay shale (Kelemen et al., 2007) are selected for atomic balance analysis. However, it should be noted that most of these lighter components are inclined to be expelled from the OM, especially for gas shale in the late catagenesis stage. In this work, we first determine the compositions of lighter components for deduction by gas chromatography of produced gas ( Wu et al., 2014), and
US
then determine the quantity for deduction by restricting the deducted methane molecules . The deducted quantity of methane molecules from the initial total lighter components is determined by adjusting the residual methane quantity to match the
AN
experimental weight fraction of saturated hydrocarbons in bitumen. Note that we include methane in both the bitumen components and the initial lighter components. On one hand, we include methane in the saturated components of bitumen to represent resid ual
M
lighter alkyl components. This is also intended to match the experimental results of bitumen group separation in the Cambrian OM. On the other hand, we include methane in the initial total lighter components so as to meet the atomic balance assumption.
ED
However, it should be noted that methane molecules in the initial lighter components are deducted later. Thus in the final OM
PT
model, there are only methane molecules included in the bitumen components. OM. Kerogen molecules, bitumen fractions (asphaltene, resin, aromatic hydrocarbons and saturated hydrocarbons) and
CE
remaining lighter components (ethane, carbon dioxide, nitrogen and water) are integrated to generate the bulk molecular model of Cambrian OM. A total of 1700 atoms and a box size of 25 Å of bulk molecular model are considered as the minimized size
AC
requirements (Ungerer et al., 2014). To meet the size requirements, 10 kerogen molecules are selected to form the skeleton of OM model. Then a suite of bitumen fractions are introduced in the system according to the experimental bitumen group separation results. The total lighter components derived from the atomic balance assumption are introduced to construct the initial version of OM model. Thereafter, a large amount of lighter components with their compositions constrained by the gas chromatography results are deducted from the initial simulation system to generate the realistic molecular model of Cambrian OM under typical reservoir conditions.
2.4 Simulation details 2.4.1 Structure optimization In this work, a combination of geometry optimization and MD simulation is utilized to generate the bulk OM model using the Materials Studio software with pcff force field. This force field describes the dispersion-repulsion interactions using 7
ACCEPTED MANUSCRIPT Lennard-Jones (LJ) 6-9 function, and describe the Coulombic interactions using force field assigned point charges (Sun et al., 1994). This ab initio all-atom force field has been parameterized, tested and validated for application to various organic materials (Sun, 1995), and has performed well for predicting the volumetric and thermodynamic properties of kerogen (Yiannourakou et al., 2013; Ungerer et al., 2014). The simulation accuracy generally increases with the cutoff distance for short-range pair interaction. However, it is well known that the cutoff distance cannot exceed the half length of the minimum dimension of the simulation box so as to meet the minimum image convention. Since the cell length of our organic models ranges from 32.04 to 34.5 Å, we adopt a
CR IP T
fine cutoff distance of 15.5 Å to truncate the short-range LJ interactions. The Ewald summation method is applied to handle the long-range electrostatic interactions. The OM fractions are first geometry optimized by the smart minimization algorithm with a fine convergence criterion. The methodology of geometry optimization refers to the work of Ru et al. (Ru et al., 2012). A max iteration step of 5000 is adopted to guarantee the energy convergence. Thereafter, these optimized OM fractions are enclosed in an
US
amorphous cell with periodic boundary conditions to generate the initial c onfigurations of OM model. To investigate the uncertainties caused by different initial configurations, 10 configurations are generated in this stage for further discussion. The cell 3
AN
density is first targeted at 0.1 g/cm (Collell et al., 2014b; Ungerer et al., 2014; Huang et al., 2017b, 2018a, 2018b ), with the system temperature being 300 K.
M
Subsequently, these configurations are optimized with the annealing strategy to obtain the global lowest energy configurations (Ru et al., 2012). A total of 10 annealing cycles of dynamic runs are performed for each configuration. Each complete annealing
ED
cycle consists of 5 heating ramps and 5 cooling ramps, for which the target temperature is gradually increased from the initial
PT
temperature (300 K) to a mid-cycle temperature (800 K) and then decreased again to the initial temperature. The dynamic time step is 1 fs, and a total simulation time of 2000 ps is performed for the annealing dynamics run. The simulation is carried out using
CE
the NVT canonical ensemble (fixed molecular number, box volumes and system temperature) with the temperature controlled by the Nosé-Hoover thermostat (Hoover, 1985). After the annealing treatment, the output configurations are further geometry
AC
optimized.
These annealed configurations are then relaxed by a succession of MD simulations, as adopted in previous molecular simulation studies (Collell et al., 2014b; Huang et al., 2017b, 2018a, 2018b). NVT ensemble is first adopted to relax these configurations at 800 K. Then succeeding NPT ensemble (fixed molecular number, system pressure and temperature) simulations are performed at 33.6 MPa and gradually decreasing temperature of 800, 600, 400 and 300 K. The optimized configurations at each temperature are output as the initial structures for the next NPT run at a lower temperature. The pressure is imposed using the Berendsen barostat (Berendsen et al., 1984), while the temperature is imposed with the Nosé-Hoover thermostat (Hoover, 1985). The simulation time is 400 ps for the NVT simulation, and 500 ps for the NPT simulations. Collell et al. (Collell et al., 2014b) adopted 12.5~250 ps for each NVT simulation and 200~400 ps for each NPT simulation. Ungerer et al. (Ungerer et al., 2014) adopted 100 ps for each NVT simulation and 100~200 ps for each NPT simulation. The simulation time in this work is much 8
ACCEPTED MANUSCRIPT longer than these in previous work (Collell et al., 2014b; Ungerer et al., 2014). The simulation time of each run has also been checked to be long enough to ensure the density convergence. At last, these configurations are further relaxed at typical reservoir conditions (381 K and 33.6 MPa) for 1000 ps to obtain the equilibrium structures.
2.4.2 Pore structure characterization In this work, free pore volumes and specific surface areas are computed with the probe insertion method (Connolly, 1983). Probe molecules with fixed Connolly radius are randomly inserted into the system. These probes roll over the atomic van der
CR IP T
Waals surfaces to determine the solid skeleton surfaces, and the regions wrapped by the skeleton surfaces are identified as free pore volumes. Similarly, a series of spherical probes of stepwise increasing Connolly radius are adopted to detect the corresponding free pore volumes and specific surface areas. Accordingly, the pore volume or specific surface area distributions of organic models are computed by differentiating the free pore volumes or specific surface areas with respect to different probe
US
diameters.
2.4.3 Adsorption behavior simulation
AN
The GCMC simulations were performed to characterize the adsorption behaviors in the OM skeleton at lighter gas species (methane, carbon dioxide and nitrogen) and heavier hydrocarbons species (ethane, propane and n-butane), with the force field,
M
pair-wise interaction, cutoff distance, thermostat method consistent with that for MD simulations. These lighter gas species are dominate species in Cambrian shale gas according to the gas chromatography results (Wu et al., 2014). The adsorption of these
ED
heavier hydrocarbons species is investigated to analyze the possible effect of adsorbate size. Adsorbate molecules are adopted as
PT
full-atom rigid models. A total of 4×107 Monte Carlo steps are performed for each pressure point, wherein the first 2×10 7 equilibrium steps are utilized as relaxation, while another 2×10 7 production steps are used as statistics. The GCMC simulations
CE
outputs the absolute adsorption amount, which can be further converted to excess adsorption capacity so as to compared with experimental data by,
(1)
AC
ne na
Where na is the absolute adsorption amount, ne is the excess adsorption amount,
is the gas density calculated by the
Peng-Robinson equation, and is the free pore volume. The isosteric adsorption heat of a component is computed to reveal the pore filling mechanism, which is defined as the partial molar enthalpy of the sorbate component in the reservoir minus that in the framework, which is determined by,
Qst RT intra N E
(2)
Where Qst is the isosteric adsorption heat, R is universal gas constant, T intramolecular chemical potential,
N
is the system temperature, intra
is the ensemble averaged molecular number of sorbate, 9
E
is the
is the ensemble averaged
ACCEPTED MANUSCRIPT total energy for the sorbate-sorbent configuration.
3 Results and discussion 3.1 Chemical structure characterization 3.1.1 FTIR Spectroscopy The TOC content for shale sample Y1 and Y2 is 2.3 wt.% and 2.1 wt.%, respectively, indicating a good organic richness for Cambrian shale (Chalmers and Bustin, 2007). The equivalent vitrinite reflectance is 1.58 % and 2.01 % for sample Y1 and Y2,
CR IP T
respectively. The thermal maturity of studied Cambrian shale can be mainly classified as high mature level. The Van Krevelen diagram is shown in Fig. S1. The sample points of Cambrian shale are much closer to the type II curve, which indicates the source rocks mainly contain type II kerogen.
In this work, shale sample Y1 and Y2 and their isolated kerogen sample K1 and K2 were chosen for FTIR spectrum analysis
US
so as to characterize the chemical functional groups. Standard FTIR absorption spectra of these samples are displayed in Fig. 1. The band assignments of investigated functional groups are listed in Table S1 (Ojima, 2003; Chen et al., 2012).
K2
1620 C=C 1455 CH2,CH3
AN
1700 C=O
1286 C-OH 429 FeS2
M
Absorbance
1375 CH3
3680 K+I
ED
K1 Y2
3500
3000
AC
4000
613 A 1150 Q+A 707 CA 512 Q 912 K
CE
PT
Y1
1059 -O- 900-730 CH x
2500 2000 1500 -1 Wavenumber (cm )
1000
500
Fig. 1. FTIR spectra of Cambrian shale and kerogen samples. Symbol description: K denotes kaolinite, Q denotes quartz, CA denotes calcite, I denotes illite, A denotes anhydrite. As seen in Fig. 1, kerogen sample K1 and K2 show basically consistent absorption spectrum, which facilitates the generation of an average molecular structure of Cambrian kerogen. Compared with the kerogen spectra, the shale spectra show similar trend, but with different peak intensity. Also, a shift of peak sites toward the lower wavenumber can be observed for the shale spectra. This peak shift can be attributed to the mineral components in the shale samples. As seen on the shale spectra, a large number of peaks related to minerals such as kaolinite, quartz and calcite can be distinguished. The shale spectrum results from the overlay of absorption peaks for both OMs and minerals. The absorption peaks for organic groups become insignificant because of the 10
ACCEPTED MANUSCRIPT influence of various minerals. Most of the absorption peaks of minerals disappear for the kerogen spectra, indicating that acid treatment can remove most of the minerals in shale. However, we can still observe the absorption peaks for pyrites on the kerogen spectra. These adsorption peaks qualitatively shows the occurrence of pyrite minerals in the kerogen samples. This observation has also been previously reported, since pyrites cannot been removed by HCl/HF acid treatment (Kelemen et al., 2007). The kerogen spectra show significant absorption peaks of organic groups with fairly correct peak sites. The peaks at 1700, 1286 and 1059 cm–1 correspond to the oxygenated groups, the peaks at 1455 and 1375 cm –1 correspond to the saturated alkyl –1
CR IP T
groups, while the peaks at 1620 and 900-730 cm correspond to aromatic groups. These peaks indicate that the chemical structure of Cambrian kerogen contains aromatic components, saturated aliphatic components and oxygen-bearing groups. The kerogen –1
spectra show remarkable strong peaks of aromatic C=C ring stretching at 1620 cm and aromatic CHx out-of-plane deformation at 900-730 cm–1, suggesting that aromatic carbon is the major carbon type in Cambrian kerogen. Also, the peak intensities of the CH3
US
asymmetric bending vibration and the CH2 shear vibration at 1455 cm–1 and the CH3 symmetric bending vibration at 1375 cm–1 are fairly strong, indicating the existence of CH3 and CH2 groups in the kerogen structure. However, the kerogen spectra show an -1
AN
absence of absorption peak of vibration of long-chain alkane ((CH2 )n , n≥4) at 720 cm , which indicates that the aliphatic carbons are mainly in the forms of short side chains or linkages. The relative peak intensities between these oxygenated groups indicate
M
that ether groups and carbonyl groups are the main oxygen-containing groups, and there is nearly no alcohol group in the kerogen structure.
13
C NMR spectroscopy was performed on kerogen sample K1 and K2 to characterize the carbon
PT
In this work, solid-state
ED
3.1.2 13 C NMR Spectroscopy
skeleton structure of Cambrian kerogen. The measured kerogen spectra are presented in Fig. 2. The kerogen spectra are composed
CE
of aliphatic carbon region (0-90 ppm), aromatic carbon region (90-165 ppm) and carbonyl carbon region (165-220 ppm). Compared with the aliphatic carbon region and carbonyl carbon region, the aromatic carbon region has the remarkably largest peak
AC
area. We can derive that aromatic carbon is the major structural fraction of Cambrian kerogen, while aliphatic and carbonyl carbon mainly serve as linkages and short side chains, which is in good agreement with the FTIR results.
11
ACCEPTED MANUSCRIPT H
125 far 135 far 145 far 156,166 far
C M
D
40,48 fal 20 f A al
O
99,105 far
A
15 fal
K2
O O
50-90 fal
H
33 fal
B
25 fal
CR IP T
175,183 fa
B
K1
250
200
150
100
50
0
US
Chemical Shift (ppm)
-50
Fig. 2. 13 C NMR spectra of Cambrian kerogen samples. Symbol description and peak assignment can refer to Table 1.
AN
To quantitatively specify the carbon groups in the chemical structure, Lorentz peak-fitting of was conducted on the 13 C NMR spectrum of kerogen sample K2. The fitting curves are shown in Fig. S2. The structural parameters and chemical shifts have been
M
reported in the work of Trewhella et al. (Trewhella et al., 1986). The fitting area fractions are listed in Table 1. The important
ED
lattice parameters have been previously defined (Kelemen et al., 2007) and are shown in Table 2. Table 1
Structural
PT
Carbon structural parameters from 13 C NMR spectrum of kerogen sample K2. Chemical
Area fraction
Carbon type
Location
shift (ppm)
f al M
15
Aliphatic methyl
Cal -C*H3
5.29
f al A
20
Aromatic methyl
Car-C*H3
1.68
f al B
25
Aliphatic C(2) carbon
C*H2 -CH3
0.43
f al H
33
Methylene
C*H2
0.71
f al D
40, 48
C*, C*H
5.97
AC
CE
parameter
(%)
Methine, Quaternary sp 3 carbon
55
Oxy-methylene
C*H2 -O
4.45
62, 69
Oxy-methine
C*H-O
10.49
83, 89
Oxy-aliphatic ring
Cal *-O
1.85
f ar A
99, 105
Oxy-aromatic ring
Car*-Car-O
7.18
f ar H
125
Protonated aromatic carbon
Car*-H
29.50
f al O
12
ACCEPTED MANUSCRIPT f ar B
135
Bridged aromatic carbon
Car*-Car
8.49
f ar C
145
Branched aromatic carbon
Car*-Cal
13.74
f ar O
156, 166
Oxy-aromatic carbon
Car*-O
4.04
faO
175, 183
Carboxyl, carbonyl
RC*=OR
6.17
Table 2 Lattice parameters derived from 13 C NMR spectrum of kerogen sample K2 Lattice
Description
Value
f al
f al M + f al A + f al B + f al H + f al D + f al O
Aliphaticity
30.87 %
f ar
f ar A + f ar H + f ar B + f ar C + f ar O
Aromaticity
62.95 %
f al B + f al H )/ f ar C
Aliphatic chain length
0.07
f ar B /( f ar A + f ar H + f ar C + f ar O )
Degree of condensation
f al D +C*H-O)/ f al
Aliphatic branch
0.53
f ar H / f ar
Aromatic substitution
0.53
BI
(
1-
US
AN
X BP
(
0.16
M
Cn
CR IP T
Definition
parameter
ED
As seen in Table 2, the aromaticity of kerogen is up to 62.95 %, 10 benzenes or 6 naphthalenes or other similar structures can be derived for each 100 carbons in Cambrian kerogen. This indicates aromatic carbons constitute the main skeleton of Cambrian
PT
kerogen structure. In addition to aromatic carbons and aliphatic carbons, the kerogen structure also contains some carboxyl or carbonyl groups, which can serve as branched side chains or short linkages between aromatic clusters. The aliphatic chain length is
CE
0.07, indicating that the aliphatic carbons mainly exist in the form of short branched chains. This conclusion can be further
AC
validated by the negligible FTIR absorption peak at 720 cm -1. The degree of condensation ( X BP ) is closely related to aromatic cluster size. The derived X BP at 0.16 for Cambrian kerogen is smaller than that of naphthalene (0.25), which indicates that the condensation degree for aromatic cluster can be represented by 1-2 aromatic rings for Cambrian kerogen. The lattice parameter for aliphatic branch ( BI ) is 0.53, suggesting that the aliphatic carbons are highly branched for Cambrian kerogen. The aromatic substitution ( =0.53) indicates that 3-4 carbons in each aromatic ring are substituted by neighboring groups.
3.1.3 XPS results In this work, the elementary forms for O, N and S in Cambrian kerogen were characterized using the XPS technique. The XPS spectra of kerogen sample K1 and K2 are illustrated in Fig. 3. The two spectra are fairly similar to each other, and the XPS spectrum of K2 is utilized to quantify the elementary forms in Cambrian kerogen. The derived element concentrations for kerogen
13
ACCEPTED MANUSCRIPT K2 from XPS spectrum are listed in Table 3. In addition to C, O, N and S, Cambrian kerogen also contains a small trace of Cl and Fe. C is the major element in the kerogen structure, with its area concentration being 89.17 %. The absorption area concentration for O is 7.32 %, much larger than that for N and S. The observed Cl derives from the remains of acid treatment, while Fe is the main component of pyrite which has also been identified in the FTIR spectra. However, since the area concentrations for both Cl
US
Cl 2s S 2s Cl 2p S 2p
800
600 400 Binding Energy (ev)
200
0
M
1000
AN
K1
O 2s
N 1s
K2
Fe 2p
O KLL
O 1s
CR IP T
C 1s
and Fe are lower than 0.2 %, their effect on the kerogen chemical structure can be neglected.
ED
Fig. 3. XPS spectra of Cambrian kerogen samples. Table 3
PT
Element concentration derived from XPS spectrum of kerogen sample K2 Binding energy(eV)
Area concentration %
Fe 2p
714.20
0.11
O 1s
531.80
N 1s
401.90
C 1s
284.50
89.17
Cl 2p
200.30
0.13
S 2p
168.80
1.61
AC
CE
Element
7.32 1.66
To quantitatively specify the elementary forms for heteroatom including O, N and S in Cambrian kerogen, the XPS spectrum for K2 sample was fitted by the XPS peak-fitting software, and the results are shown in Table 4 and Fig. S3. The binding energies corresponding to functional groups have been documented in previous work (Tong et al., 2011; Li et al., 2014). Due to the effect of N and S, the peak-fitting of O 1s cannot accurately determine the binding modes between C and O. Therefore, the spectrum of C 1s signal was fitted to quantify the elemental forms of O associated with C atom. The results indicate that ether or alcohol groups are the main oxygen-containing groups in the kerogen structure. Also, the kerogen contains some carbonyl and carboxyl groups.
14
ACCEPTED MANUSCRIPT The analysis on N 1s signal indicates that amine groups are the dominant nitrogen-containing groups in the kerogen structure. Besides, the Cambrian kerogen contains a small amount of stable pyridine and pyrrole groups. The nitrogen oxide groups introduced from the acid treatment are observed to be well removed. The fitted results on S 2p signal indicate that S mainly exists in the organic components, and sulfoxide group is the dominant sulfur-containing group, with its area ratio being up to 86.43 %, followed by the thiophene group at 13.57 %. Table 4
Binding energy (eV)
Peak area ratio (%)
C 1s
C-C, C-H
284.51
79.70
C-O, C-OH
285.98
11. 79
C=O
287.05
3.13
COOH, COOR
288.45
2.37
Pyridine
398.09
15.80
Pyrrole
400.05
7.46
Amine
401.45
73.38
NOx
403.01
3.35
Sulfoether
162.13
0.00
Pyrite
163.10
Thiophene
164.09
Sulfoxide
168.97
0.00
13.57 86.43
PT
3.2 Molecular models
AN
S 2p
ED
N 1s
CR IP T
Functional group
US
Element
M
Elementary forms for heteroatom derived from XPS spectrum for K2 sample
3.2.1 Molecular structure of Cambrian kerogen
CE
Based on the analytical structural parameters, an initial 2D molecular structure of Cambrian kerogen was generated. The
13
C
NMR spectrum was predicted using the ACD/CNMR Predictor software, and compared with the experimental spectrum. The
AC
aliphatic carbons and heteroatom groups on the kerogen structure were continuously adjusted until good agreement was observed for the predicted and experimental spectrum (Fig. 4). The final representative Cambrian kerogen structure (C210 H184 O20 N4 S4) is illustrated in Fig. 5, with its compositions and structural parameters presented in Table 5. The kerogen atomic compositions are in good agreement with the elemental analys is results. The kerogen skeleton consists of 9 naphthalenes and 6 benzenes linked by branched aliphatic carbons and heteroatom groups. The lattice parameters of the kerogen structure such as degree of condensation, aliphatic branch and substitution are fairly close to the experimental
13
C NMR results. The oxygenated groups are mainly
composed of ether, carbonyl and carboxyl groups, consistent with the XPS and FTIR results. Moderate deviations are observed for N and S groups due to their low contents. The 4 N atoms are distributed in 3 amine groups and 1 pyridine group, while the 4 S atoms are involved in 3 sulfoxide groups and 1 thiophene group. In this work, 10 kerogen units are adopted to build the skeleton of
15
ACCEPTED MANUSCRIPT the Cambrian OM model, meeting the reported minimized size requirement for 1700 atoms per unit cell (Collell et al., 2014b).
006 Simulation Experiment
004
CR IP T
002
000 180
160
140
120 100 80 Chemical Shift (ppm) 13
40
C NMR spectrum of Cambrian kerogen K2.
PT
ED
M
AN
Fig. 4. Predicted and experimental
60
US
200
Table 5
CE
Fig. 5. Molecular structure of Cambrian kerogen.
Property
AC
Compositional and structural parameters of Cambrian kerogen
a
Structural parameter
Experimental
Modelling
H/C
0.87
0.88
Atomic
O/C
0.093
0.095
compositions
N/C
0.019
0.019
S/C
0.018
0.019
fal
30.87 %
36.19 %
f ar
62.95 %
61.43 %
Cn
0.07
0.18
X BP
0.16
0.16
C groups
16
20
0
ACCEPTED MANUSCRIPT
N groups
S groups
a
0.53
0.53
0.53
0.50
Carbonyl (mol % C~O)
15.89
17.64
Ether (mol % C~O)
59.99
58.82
Carboxyl (mol % C~O)
24.12
23.53
Pyridine (mol % N)
15.80
25.00
Pyrrole (mol % N)
7.46
0.00
Amine (mol % N)
73.38
75.00
NOx
3.35
0.00
Sulfoether (mol % organic S)
0.00
0.00
Thiophene (mol % organic S)
13.57
25.00
Sulfoxide (mol % organic S)
86.43
75.00
CR IP T
O groups
BI
Experimental ato mic co mpositions are derived fro m elemental analysis. Experimental C groups are derived fro m 13 C NMR spectrum.
Experimental O, N and S groups are derived from XPS spectrum.
US
3.2.2 Structures and compositions of bitumen components
The bitumen components are composed of representative components, with their compositions listed in Table 6. Note that
AN
resin is a derivative of higher plant. The OM in shale from the Lower Cambrian Niuititang Formation is mainly derived from
M
pelagic algae, fungi and zooplankton. Thus the content of resin is relatively low in our investigated OM. However, we integrate the resin in our OM system for the purpose of generating a complete bulk molecular model of OM. The molecular structures of resin
ED
unit and fichtelite unit are shown in Fig. 6. The molecular size (C33 H33 ON) of our resin unit refers to the documented one (C33 H41 O1.6 N0.4 S0.1) in the work of Tissot and Welte (Tissot and Welte, 1984). The topology including aromatic rings and alkane
PT
tails of our resin unit refers to the reported resin structure in the work of Collell et al. ( Collell et al., 2014b). The molecular
CE
structures and quantities of bitumen components are adjusted to match the experimental weight fraction. Good agreement between the simulated and experimental weight fraction distribution can be observed. The converted total bitumen content (0.0547 wt. %)
AC
is also fairly close to the experimental result (0.0511 wt. %). In this work, asphaltene components are excluded because of their low content, and saturated hydrocarbons and resin components are the dominant components in Cambrian bitumen.
17
ACCEPTED MANUSCRIPT (b)
CR IP T
(a)
Fig. 6. Molecular units of resin and fichtelite in Cambrian shale. (a) resin unit (C33 H33 ON); (b) fichtelite unit (C14 H24 ). Atom representation: black for C atom, red for O atom, blue for N atom and grey for H atom.
Weight composition of bitumen components in Cambrian OM
US
Table 6
Experimental data a
Representative
Molecule
Weight fraction
components
molecule
numbers
(%)
(%)
resin
resin
1
42.94
40.34
asphaltene
asphaltene
0
0
6.45
aromatics hydrocarbons
toluene
1
8.61
5.75
fichtelite
1
48.46
47.46
methane
20
a
M
ED
PT
saturated hydrocarbons
AN
Bitumen
Experimental results from bitumen group separation on Cambrian sample Y2.
CE
3.2.3 Molecular model of Cambrian OM
Associated with kerogen units and bitumen components, lighter components are also included to generate the initial Cambrian
AC
OM model adopting the assumption of atomic balance. The component distribution in the initial organic model is presented in Table S2. The validation of atomic balance is shown in Table 7. Reasonable consistencies for atomic compositions are observed between the initial organic model and an immature type II kerogen of Duvernary shale (Kelemen et al., 2007). The final Cambrian OM model is obtained by deducting a suite of lighter components prone to release from the organic skeleton. The compositions of deducted lighter components are restricted by gas chromatography results on Cambrian shale gas (Wu et al., 2014) (Table 8). The deducted quantity of methane is determined by adjusting the residual methane quantity to match the experimental weight fraction of saturated hydrocarbons in bitumen. The component distribution and atomic distribution of the final Cambrian OM model are shown in Table 9 and Table S3, respectively. Kerogen molecules dominate the compositions in the OM model. The final OM model corresponds to type II organic type (Fig. S1).
18
ACCEPTED MANUSCRIPT Table 7 Atomic balance between initial Cambrian OM model and immature kerogen phase Initial OM model Immature kerogen a
a
H/C
O/C
N/C
S/C
1.17
0.097
0.029
0.017
1.17
0.097
0.029
0.014
Analytic data of an immature kerogen of Duvernary shale documented by Kelemen et al. (Kelemen et al., 2007).
Table 8 Experimental and deducted lighter components distribution Molecular number
Mole fraction (%)
Experimental data a (%)
methane
194
96.52
96.52
ethane
1
0.50
0.35
nitrogen
2
1.00
1.24
carbon dioxide
4
1.99
1.75
Analytic data of gas chromatography on produced gas from Cambrian shale documented by Wu et al. (Wu et al., 2014).
US
a
CR IP T
Representative molecule
Components distribution in the final Cambrian OM model
AN
Table 9
Molecular number
Weight fraction (%)
kerogen
10
94.56
resin
1
1.35
toluene
1
0.27
fichtelite
1
methane
20
nitrogen
12
carbon dioxide
6
ED 0.58
PT
0.94
CE
water
M
Representative molecule
10
0.99 0.78 0.53
AC
3.3 Microstructure characterization 3.3.1 Properties of steric structure
In this section, the steric molecular model of Cambrian OM (Fig. 7) was constructed based on the compositions in Table 9 and the simulation details for structure optimization in section 2.4.1. As seen in Fig. 7, the OM structure is highly unordered and no remarkable stacking of aromatic units can be observed. This observation can be attributed to the chemical structure of Cambrian kerogen, as shown in Fig. 5. Cambrian kerogen contains a small size of polyaromatic clusters, which can be represented by 1 or 2 benzene rings. These benzene rings are connected by bendable aliphatic branches and heteroatom groups, which result in a strong steric hindrance on the stacking of aromatic units. Interestingly, we observe that lighter components such as water, methane, nitrogen and carbon dioxide tend to appear in pairs around heteroatom groups, showing competitive adsorption behaviors due to surface chemistry heterogeneity. Instead of aggregating into clusters, water molecules are scattered around hydrophilic groups, 19
ACCEPTED MANUSCRIPT which can be attributed to the abundant polar groups and poor inter-pore connectivity. On one hand, water molecules are preferentially adsorbed on hydrophilic groups on the OM skeleton. The OM skeleton contains sufficient polar groups to accommodate these water molecules considering the low moisture content. Thus, water molecules are prone to locate around thes e polar groups. On the other hand, the OM skeleton is highly unordered, and the effective inter-pores are poorly connected, which restricts the transport of water molecules to aggregate into clusters. The random distribution of bitumen components and lighter
M
AN
US
CR IP T
components aggravates the disorder of pore network in the OM system.
ED
Fig. 7. View of simulation box of the OM models at 381 K and 33.6 MPa. Molecule representation: cyan for resin, grey for fichtelite, red for toluene, green for methane, purple for water, blue for nitrogen, yellow for carbon dioxide.
PT
3.3.2 Properties of porous structure
As shown in Table 10, the averaged physical densities of our organic models are in the range of 1.22±0.01-1.30±0.03 g/cm3,
CE
fairly close to the reported experimental data (1.18-1.25 g/cm3) of high mature type II kerogen (Okiongbo et al., 2005). The helium probed porosities of these organic models range between 10.95±0.3 and 21.21±0.4 %, which are consistent with the calculated OM
AC
porosities (4.45-22.50 %) in Barnett mudstone based on SEM images (Loucks et al., 2009). In addition to the original annealing temperature cycle (300~800~300 K), a new temperature cycle (300~900~300 K) is adopted to relax the OM structure to account for possible temperature cycle effect. The mid-temperature is increased to be 900 K, consistent with the documented highest temperature for OM structure relaxations (Collell et al., 2014b; Ungerer et al., 2014). The OM structure under this temperature cycle is referred as OMe in Table 10. The comparison results show that the fluctuations of OM structure properties under the two different temperature cycles are tiny and acceptable. The unordered combination of kerogen molecular units constitutes the kerogen matrix, which is also defined as the OM skeleton. These lighter species and heavier species are included to represent the residual components in the porous network of OM skeleton. In Table 10, we present the effect of the inclusion of these species on the OM structure properties. The result shows that the densities increase while the porosities decrease for our organic models with the inclusion of these species. This observation suggests that the porous network of OM is mainly constituted by the kerogen skeleton, and these inclusions occupy some pore volumes, contributing to the increasing condensation of the OM system. In Table 10, we also present the effect of probe size on the structure properties of organic models. The porosities and specific surface areas probed by helium are far smaller than these of innate porous structures with a probe diameter of zero. Note that the specific surface 20
ACCEPTED MANUSCRIPT 2
areas exceed 4000 m /g, and the porosities can be up to 40 % under a probe diameter of zero. These high porosities and specific surface areas are mainly contributed by tiny pores, which cannot be accessed by gas molecules. These tiny pores cannot be measured in kerogen samples in laboratory. However, these tiny pores should not be neglected since they may be expanded due to kerogen swelling. When the probe size increases to that of helium, the porosities and specific surface areas show a significant decrease. The porosities decline to comparable values with experimental results, while the specific surface areas (1387-1524 m2 /g) are still higher than measured results. The higher specific surface area can be attributed to the narrow pore size distributions of our organic models. The pore size in actual kerogen matrix ranges from several angstroms up to a few micrometers. These pores wit h nanometer and micrometer sizes greatly reduces the specific surface areas. In our Cambrian OM system, the ultra-micropore is the dominant pore type. The small differences of innate porous properties between the externally accessible regions and the total
CR IP T
regions indicate that most of pores in the organic structures are externally accessible. The density, porosity and specific surface area diverge slightly under different initial structure configurations, with the uncertainties less than 5 %, which shows a good reproductivity.
Properties of different organic models at 300 K and 0.1 MPa a
Porosity (%)
3
Density (g/cm )
VDW b
VDW c
He
VDW b
VDW c
AN
He
Specific surface area (m2 /g)
1.30 ± 0.03
10.95 ± 0.3
36.24 ± 1.6
36.36 ±1.3
1524 ± 64
4838 ± 203
4844 ± 205
OM e
1.28 ± 0.02
11.21 ± 0.3
36.62 ± 1.5
36.73 ±1.2
1497 ± 59
4789 ± 177
4812 ± 192
Kerogen+lighter fractions
1.27 ± 0.02
13.84 ± 0.3
36.93 ± 1.6
37.07 ± 1.4
1389 ± 55
4379 ± 182
4414 ± 189
Kerogen+heavier fractions
1.25 ± 0.01
17.78 ± 0.4
38.08 ± 1.7
38.28 ± 1.6
1423 ± 62
4221 ± 191
4251 ± 177
Kerogen in OM
1.22 ±0.01
21.21 ± 0.4
40.12 ± 1.8
40.32 ± 1.9
1641 ± 77
4217 ± 178
4236 ± 192
PT
a
M
OM d
ED
Structure
US
Table 10
Probe diameter of He molecule is 2.6 Å. VDW b refers to a probe diameter of zero defined over the externally accessible regions, while
CE
VDW c is equivalent to a probe diameter of zero defined over the whole pore regions. OM d refers to the OM structure relaxed under the annealing cycle (300~800~300 K), while OM d refers to the OM structure under the annealing cycle (300~900~300 K). The average deviations are determined on 10 structure configurations of organic models.
AC
To further analyze the effect of lighter fractions (methane, water, nitrogen and carbon dioxide) and heavier fractions (resin, fichtelite and toluene), the pore volume and specific surface area distributions of various organic models are compared in Fig. 8. Note that the pore size distribution of our organic models are narrow due to the finite sizes of the molecular models. This k ind of narrow pore size distribution of organic models has also been reported in previous work of molecular simulations ( Collell et al., 2014b; Ho et al., 2018). These pore size distributions are close to the experimental micropore size distributions on kerogen calculated from carbon dioxide adsorption (Rexer et al., 2014). Therefore, the pore network in our OM model can reasonably represent the ultra-microporous and limited microporous network of OM in shale formations. In this work, most of the pore volumes and specific surface areas are contributed by the ultra-micropores with pore width ranging between 1 and 3 Å. With the introduction of inclusions in the OM system, the pores larger than 3.3 Å decrease while the smaller pores increase. This means that
21
ACCEPTED MANUSCRIPT the introduced fractions enter these larger pores and divide them into inaccessible pores. It is interesting to note that compared with lighter fractions, the heavier fractions have a much smaller effect on the kerogen porous structure. The pore distributions of different organic models are illustrated in Fig. 9. The pore distribution of kerogen containing heavier frac tions (Fig. 9b) is fairly similar to that of innate kerogen structure (Fig. 9d). By contrast, a large number of pores in the kerogen structure are occupied by lighter fractions (Fig. 9c). The heavier fractions are more inclined to be part of the porous skeleton, while the lighter fractions tend to distribute in enterable pores attracted by heteroatom groups on pore surfaces. Moreover, these helium probed pores in the OM
CR IP T
system are observed to be poorly connected. We conclude previously that most of pores in the organic structures are externally accessible. Note that this conclusion applies to the probe with a diameter of zero, which means that throats at any size in t he OM skeleton are included to analyze the connectivity. However, when helium is used as probe molecules, the connectivity is analyzed based on the throats that have bigger sizes than the helium molecule. Therefore, we can conclude that the throats in the OM
US
skeleton are highly connected, but most of them have smaller sizes than the helium molecule. However, these smaller throats cannot be neglected because they may become accessible for gas transport when considering the swelling effect. (b) OM Kerogen+lighter frac. Kerogen+heavier frac. Kerogen in OM
120
M
-2
80
ED
60 40
0
1
2
3 4 Pore width, Å
PT
20 0
AN
160 -2
100
dV/dr, 10 Å
2
120
200
dS/dr, 10 Å
(a) OM Kerogen+lighter frac. Kerogen+heavier frac. Kerogen in OM
140
5
6
80 40 0
0
1
2
3 4 Pore width, Å
5
6
AC
CE
Fig. 8. Pore volume and specific surface area distribution. (a) Pore volume distribution. (b) Specific surface area distribution.
22
ACCEPTED MANUSCRIPT (b)
(c)
(d)
AN
US
CR IP T
(a)
Fig. 9. Pore distribution of various organic models probed by He. Skeleton atoms are wrapped by the purple pore surfaces. (a) The
M
OM model. (b) The kerogen+heavier fractions model. (c) The kerogen+lighter fractions model. (d) The kerogen model.
3.4 Gases adsorption behaviors in OM skeleton
ED
In this section, the adsorption characteristics of lighter gas species (methane, carbon dioxide and nitrogen) and heavier hydrocarbons species (ethane, propane and n-butane) in the porous OM skeleton were investigated by the GCMC simulations. As
PT
shown in Table 8, methane, carbon dioxide and nitrogen are dominate species based on the analytic data of gas chromatography on produced gas from Cambrian shale (Wu et al., 2014). We first check the adsorption behaviors of these dominated gas species in
CE
our OM skeleton. The absolute adsorption isotherms present a trend of type I isotherm according to the classification of Gregg and Sing (Gregg and Sing, 1982) (Fig.10a). The absolute adsorption amounts increase with pressure, and the adsorption capacity is carbon dioxide > methane > nitrogen. To compare with experimental adsorption results, the excess adsorption amounts are
AC
computed for these dominate gas species (Fig.10b). The experimental excess adsorption data for methane, carbon dioxide and nitrogen are performed on the high-mature kerogen sample in Barnett shale (Hu, 2014), Ruhr coal sample (Gensterblum et al., 2013) and Sulcis coal sample (Ottiger et al., 2008), respectively. Note that we aim to validate the reasonability of our OM model to represent the adsorption characteristics of porous OM skeleton. It is well known that the porous structures of OM skeletons in coal and shale are similar at molecular level when neglecting the effect of micro-fractures. Therefore, the experimental data on coal can be used to be compared with our simulated results.
It is well known that some minerals cannot be removed from the
experimental OM samples. The experimental data are normalized by the TOC so as to eliminate the effect of residual minerals before the comparison. The simulated excess adsorption isotherms of methane, carbon dioxide and nitrogen on our OM skeleton match fairly well with experimental results on authentic organic pores in both trend and magnitude. Note that it is extremely difficult or impossible to represent the whole range of pore sizes of OM under molecular and atomistic scale. As discussed in section 3.3.2, the OM skeleton in this work can represent the ultra-microporous and limited microporous network of OM in shale samples. Thus, this observation in Fig.10b indicates that ultra-micropores and limited micropores dominate the pore types for 23
ACCEPTED MANUSCRIPT adsorption of methane, carbon dioxide and nitrogen in realistic OM skeleton. Therefore, the generated OM skeleton in this work can be used to study the adsorption behaviors of methane, carbon dioxide and nitrogen. The simulated maximum excess adsorption amount of carbon dioxide (1.80 mmol/g) is fairly close to the reported one (1.84 mmol/g). The excess adsorption amount of car bon dioxide initially increases and then decreases with the increase of pressure. Similar with the carbon dioxide isotherm, the simulated methane and nitrogen excess adsorption isotherms are comparable with the experimental results. Compared with the documented experimental data, the simulated methane isotherm presents a higher adsorption rate, and a slight decrease when the press ure exceeds ~9 MPa can be observed. Notably, the later drop of excess adsorption amount is not observed for the nitrogen isotherm . The excess adsorption is closely correlated with the density difference between the adsorbed layer and the bulk. The excess adsorption amount reaches the maximum when there is a maximal density difference. The decrease in excess isotherms is because
CR IP T
the density difference between the adsorbed layer and the bulk reduces at high pressures. It can be concluded that the decrease in excess isotherms can be affected by the gas adsorption capacity, and a higher adsorption capacity corresponds to a lower pressure when the decrease happens. To investigate the thermodynamic property of gas adsorption on the OM skeleton, the isosteric heat of gases adsorption is presented in Fig.10c. The isosteric heat of gases adsorption reflects the combined effect of adsorbent-adsorbate interaction and adsorbate-adsorbate interaction. The isosteric heat of carbon dioxide adsorption initially decreases and then gently
US
increases with further rising pressure, while that of methane and nitrogen adsorption shows continually slight increase as the pressure rises. The early decrease of isosteric heat of carbon dioxide adsorption reflects the chemically heterogeneity of OM
AN
skeleton, since carbon dioxide preferentially fills these high-energy adsorption sites and then occupies these low-energy sites. The subsequent increase of isosteric heat of gases adsorption can be attributed to the increasing interaction between gas molecules. The limit isosteric heat, defined as the isosteric heat of gas adsorption at zero coverage (gas adsorption approaches zero), reflects
M
quantitatively the adsorption enthalpy between adsorbent and adsorbate without the influence from the adsorbate-adsorbate interaction. The limit isosteric heat on the OM skeleton is carbon dioxide (27.33 KJ/mol) > methane (20.76 KJ/mol)> nitrogen
ED
(17.51 KJ/mol), consistent with the order of gas adsorption capacity.
PT
Although there are nearly no heavier hydrocarbon species in the Cambrian shale gas (Table 8), the adsorption isotherms of ethane, propane and n-butane are simulated and compared with the experimental results in the work of Zhao et al. (Zhao et al., 2018) so as to check the possible effect of pore size or adsorbate size. As shown in Fig.10d-f, the simulation results are in poor
CE
agreement with the experimental results. The simulated adsorption amounts reach the plateaus at lower pressures, and the simulated maximum excess adsorption amounts are much lower than the measured results. These observations indicate that the
AC
adsorption behaviors of heavier hydrocarbon species cannot be represented in our OM skeleton with ultra-micropores and limited micropores, and the contribution of some bigger pores in realistic OM cannot be neglected. This conclusion shows an effect of adsorbate size on the dominated pore size for adsorption. This effect of adsorbate size can also be reflected by the order of absolute adsorption amounts of these heavier hydrocarbon species (ethane > propane > n-butane). The accessible pores in the OM skeleton decrease as the molecular sizes of these hydrocarbon species increase.
24
Excess adsorption amount, mmol/g
2.5 2.0 1.5 CH4
1.0
CO2 N2
0.5
0
4
8 12 Pressure, MPa
16
30
1.5
1.0 0.5
0.0 0
20
CO2
N2
20
4
8 12 Pressure, MPa
16
20
1.5 1.0 0.5 0.0
-0.5
0.0 0.0
ED
0.5 1.0 Pressure, MPa
1.5
Adsorption amount, mmol/g
CE
0.5
PT
1.0
CO2_exp.
N2_sim.
N2_exp.
16
20
(d)
0
4
8 Pressure, MPa
12
16
1.2
(e)
Absolute adsorption_sim. Excess adsorption_sim. Excess adsorption_exp.
AC
Adsorption amount, mmol/g
1.5
CH4_exp.
CO2_sim.
8 12 Pressure, MPa
M
0
CH4_sim.
Absolute adsorption_sim. Excess adsorption_sim. Excess adsorption_exp.
AN
CH4
4
2.0
(c)
25
15
(b)
Adsorption amount, mmol/g
0.0
2.0
CR IP T
(a)
US
3.0
Isosteric heat of adsorption, KJ/mol
Absolute adsorption amount, mmol/g
ACCEPTED MANUSCRIPT
0.9
(f)
Absolute adsorption_sim. Excess adsorption_sim. Excess adsorption_exp.
0.6
0.3
0.0 0.0
0.1
0.2 0.3 Pressure, MPa
0.4
0.5
Fig. 10. Adsorption characteristics of gases in OM skeleton. (a) Absolute adsorption amount of methane, carbon dioxide and nitrogen (338 K). (b) Excess adsorption amount of methane (338 K), carbon dioxide (318 K) and nitrogen (343 K). (c) Isosteric heat of gases adsorption (338 K). (d) Adsorption amount of ethane (333 K). (e) Adsorption amount of propane (333 K). (f) Adsorption amount of n-butane (333 K).
4 Conclusions In this work, we proposed a novel strategy for constructing hybrid shale OM model at molecular scale. A realistic molecular model of the Cambrian OM was constructed based on analytical structural and compositional results. Microstructure
25
ACCEPTED MANUSCRIPT characterization of the OM model and gases adsorption characteristics in the OM skeleton were further discussed. The main conclusions are summarized below. The Cambrian OM model is composed of kerogen, bitumen components and residual lighter components. The OM model shows a reasonable representation of the Cambrian OM in terms of structural parameters, generic compositions, physical density and porosity. Thus, it can be utilized to study the Cambrian OM at molecular scale. The steric OM structure is highly unordered and shows a remarkable folding of aromatic units. Lighter fractions tend to
CR IP T
appear in pairs around heteroatom groups, showing competitive adsorption behaviors. Instead of aggregating into clusters, wat er molecules are scattered around hydrophilic groups. Externally accessible ultra-micropore is the dominant pore type in the Cambrian OM system, while the pores accessible for gas molecules are poorly connected. The heavier bitumen fractions are inclined to be part of the porous skeleton, exerting a small effect on the porous structure. By contrast, the lighter fractions tend to
US
distribute in enterable pores due to the attraction of heteroatom groups, thus dividing these pores into massive inaccessible pores. The GCMC simulations show that the adsorption capacity in the OM skeleton is carbon dioxide > methane > nitrogen, as
AN
validated by the adsorption amount and limit isosteric heat. The adsorption behaviors of methane, carbon dioxide and nitrogen can be well represented by the OM skeleton. The gas adsorption capacity affects the corresponding pressure when the excess isotherm
M
reaches the maximum, and a higher adsorption capacity corresponds to a lower pressure at the maximum. There is a significant effect of adsorbate size on the adsorption behaviors of heavier hydrocarbon species (ethane, propane and n-butane). To represent
ED
the adsorption behaviors of heavier hydrocarbon species reasonably, the OM system should be augmented to create a wider pore
PT
size distribution.
The reported novel construction process for shale OM model demonstrates potential application for the representative
CE
structure reproduce of realistic shale OM under reservoir conditions at molecular level. The generated OM model can serve as a starting point for further theoretical investigations on the fundamental mechanisms governing the storage, transport, and rec overy
AC
of shale gas in heterogeneous organic nanopores.
Appendix. A. Supplementary data Supplementary data contain Table S1-S3 and Fig. S1-S3 associated with this manuscript.
Notes The authors declare no competing financial interest.
Acknowledgments This work was supported by the National Natural Science Foundation of China (Grant No. 51774298 and 51504265) and the Science Foundation for the Excellent Youth Scholars of China University of Petroleum (Beijing) (Grant No.2462015YQ0223). Computer time for this study was provided by the HP High Performance Computing Cluster of the State Key Laboratory of Heavy
26
ACCEPTED MANUSCRIPT Oil Processing at China University of Petroleum (Beijing).
References Bajc, S., Cvetkovic, O., Ambles, A., Vitorović, D., 2008. Characterization of type III kerogen from Tyrolean shale (Hahntennjoch, Austria) based on its oxidation products. J. Serb. Chem. Soc. 73, 463-478. Berendsen, H.J.C., Postma, J.P.M., van Gunsteren, W.F., DiNola, A., Haak, J.R., 1984. Molecular dynamics with coupling to an external bath. J. Chem. Phys. 81, 3684-3690. Bous ige, C., Ghimbeu, C.M., Vix-Guterl, C., Pomerantz, A.E., Suleimenova, A., Vaughan, G., Garbarino, G., Feygenson, M.,
Mater. 15, 576-582.
CR IP T
Wildgruber, C., Ulm, F. J., Pellenq, R.J.M., Coasne, B., 2016. Realistic molecular model of kerogen/'s nanostructure. Nat.
Budinova, T., Razvigorova, M., Tsyntsarski, B., Petrova, B., Ekinci, E., Yardim, M.F., 2009. Characterization of Bulgarian oi l shale kerogen revealed by oxidative degradation. Chem. Erde-Geochem. 69, 235-245.
Chalmers, G.R.L., Bustin, R.M., 2007. The organic matter distribution and methane capacity of the Lower Cretaceous strata of Northeastern British Columbia, Canada. Int. J. Coal Geol. 70, 223–239.
US
Chen, S., Zhu, Y., Wang, H., Liu, H., Wei, W., Fang, J., 2011. Shale gas reservoir characterization: A typical case in the southern Sichuan Basin of China. Energy 36, 6609-6616.
Chen, Y., Mastalerz, M., Schimmelmann, A., 2012. Characterization of chemical functional groups in macerals across different
AN
coal ranks via micro-FTIR spectroscopy. Int. J. Coal Geol. 104, 22-33.
Collell, J., Galliero, G., Gouth, F., Montel, F., Pujol, M., Ungerer, P., Yiannourakou, M., 2014a. Molecular simulation and
conditions. Micropor. Mesopor. Mat. 197, 194-203.
M
modelisation of methane/ethane mixtures adsorption onto a microporous molecular model of kerogen under typical reservoir
ED
Collell, J., Ungerer, P., Galliero, G., Yiannourakou, M., Montel, F., Pujol, M., 2014b. Molecular simulation of bulk organic matter in Type II shales in the middle of the oil formation window. Energ. Fuel. 28, 7457-7466.
PT
Connolly, M.L., 1983. Solvent-accessible surfaces of proteins and nucleic acids. Science 221, 709-713. Curtis J.B., 2002. Fractured shale-gas systems. Am. Assoc. Pet. Geol. Bull. 86, 1921-1938. Feng, F., Akkutlu, I.Y., 2015. Flow of hydrocarbons in nanocapillary: A non-equilibrium molecular dynamics study. SPE Asia
CE
Pacific Unconventional Resources Conference and Exhibition, Brisbane, Australia, 9-11 November, SPE-177005. Gensterblum, Y., Merkel, A., Busch, A., Krooss, B.M., 2013. High-pressure CH4 and CO2 sorption isotherms as a function of coal
AC
maturity and the influence of moisture. Int. J. Coal Geol. 118, 45–57.
Gregg, S.J., Sing, K.S., 1982. Adsorption, surface area and porosity. Academic Press, New York, p.4. Ho, T.A., Wang, Y., Criscenti, L.J., 2018. Chemo-mechanical coupling in kerogen gas adsorption/desorption. Phys. Chem. Chem. Phys. 20, 12390-12395. Hoover, W.G., 1985. Canonical dynamics: Equilibrium phase-space distribution. Phys Rev A Gen Phys 31, 1695-1697. Hu, H., 2014. Methane adsorption comparison of different thermal maturity kerogens in shale gas system. Chin. J. Geochem. 33, 425-30.
Huang, L., Ning, Z.F., Wang, Q., Ye, H.T., Qin, H.B., 2017a. Molecular Simulation of CO2 sequestration and enhanced gas recovery in gas rich shale: An insight based on realistic kerogen model. SPE Abu Dhabi International Petroleum Exhibition & Conference, Abu Dhabi, UAE, 13-16 November, SPE-188216. Huang, L., Ning, Z.F., Wang, Q., Qi, R.R., Li, J., Zeng, Y., Ye, H.T., Qin, H.B., 2017b. Thermodynamic and structural characterization of bulk organic matter in Chinese Silurian Shale: Experimental and molecular modeling studies. Energ. Fuel. 31, 4851-4865.
27
ACCEPTED MANUSCRIPT Huang, L., Ning, Z.F., Wang, Q., Zhang, W.T., Cheng, Z.L., Wu, X.J., Qin, H.B., 2018a. Effect of organic type and moisture on CO2 /CH4 competitive adsorption in kerogen with implications for CO 2 sequestration and enhanced CH4 recovery, Appl. Energy 210, 28–43. Huang, L., Ning, Z.F., Wang, Q., Qi, R.R., Zeng, Y., Qin, H. B., Ye, H.T., Zhang, W.T., 2018b. Molecular simulation of adsorption behaviors of methane, carbon dioxide and their mixtures on kerogen: Effect of kerogen maturity and moisture content. Fuel 211, 159-172. Kazemi, M., Takbiri-Borujeni, A., 2016. Molecular dynamics study of carbon dioxide storage in carbon-based organic nanopores. SPE Annual Technical Conference and Exhibition, Dubai, UAE, 26-28 September, SPE-181705. Kelemen, S.R., Afeworki, M., Gorbaty, M.L., Sansone, M., Kwiatek, P.J., Walters, C.C., Freund, H., Siskin, M., 2007. Direct
CR IP T
characterization of kerogen by X-ray and solid-state 13 C nuclear magnetic resonance methods. Energ. Fuel. 21, 1548-1561. Kim, T.H., Cho, J., Lee, K.S., 2017. Evaluation of CO 2 injection in shale gas reservoirs with multi-component transport and geomechanical effects. Appl. Energy 190, 1195-1206.
Kurniawan, Y., Bhatia, S.K., Rudolph, V., 2006. Simulation of binary mixture adsorption of methane and CO 2 at supercritical conditions in carbons. AIChE J. 52, 957-967.
US
Li, J., Li, X.F., Wang, X.Z., Li, Y.Y., Wu, K.L., Shi, J.T., Yang, L., Feng, D., Zhang, T., Yu, P.L., 2016a. Water distribution characteristic and effect on methane adsorption capacity in shale clay. Int. J. Coal Geol. 159, 135-154.
AN
Li, J., Li X.F., Wu, K.L., Wang, X.Z., Shi, J.T., Yang, L., Zhang, H., Sun, Z., Wang, R., Feng, D., 2016b. Water sorption and distribution characteristics in shale clay: Effect of surface force. Energ. Fuel. 30, 8863−8874. Li, J., Li, X.F., Wu, K.L., Feng, D., Zhang, T., Zhang, Y.F., 2017. Thickness and stability of water film confined inside nanoslits
M
and nanocapillaries of shale and clay. Int. J. Coal Geol. 179, 253-268.
Li, Q., Han, X., Liu, Q. Jiang, X., 2014. Thermal decomposition of Huadian oil shale. Part 1. Critical organic intermediates, Fuel
ED
121, 109-116.
Li, S., Yue, C., 2003. Study of pyrolysis kinetics of oil shale. Fuel 82, 337-342.
PT
Liu, S., Deng, B., Zhong, Y., Ran, B., Yong, Z., Sun, W., Yang, D., Jiang, L., Ye Y., 2016a. Unique geological features of burial and superimposition of the Lower Paleozoic shale gas across the Sichuan Basin and its periphery. Earth Sci. Front. 23, 11-28.
CE
Liu, X.Q., He, X., Qiu, N.X., Yang, X., Tian, Z.Y., Li, M.J., Xue, Y., 2016b. Molecular simulation of CH4 , CO2, H2O and N2 molecules adsorption on heterogeneous surface models of coal. Appl. Surf. Sci. 389, 894-905. Loucks, R.G., Ruppel, S.C., 2007. Mississippian Barnett Shale: Lithofacies and depositional setting of a deep-water shale-gas
AC
succession in the Fort Worth Basin, Texas. Am. Assoc. Pet. Geol. Bull. 91, 579−601. Loucks, R.G., Reed, R.M., Ruppel, S.C., Jarvie, D.M., 2009. Morphology, genesis and distribution of nanometer -scale pores in siliceous mudstones of the Mississippian Barnett Shale. J. Sediment. Res. 79, 848-861. Lu, X.Q., Jin, D.L., Wei, S.X., Zhang, M.M., Zhu, Q., Shi, X.F., Deng, Z.G., Guo, W.Y., Shen, W.Z., 2015. Competitive adsorpt ion of a binary CO2 -CH4 mixture in nanoporous carbons: effects of edge-functionalization. Nanoscale 7, 1002-1012. Middleton, R.S., Carey, J.W., Currier, R.P., Hyman, J.D., Kang, Q., Karra, S., Martínez, J.J., Porter, M.L., Viswanathan, H.S ., 2015. Shale gas and non-aqueous fracturing fluids: Opportunities and challenges for supercritical CO 2 . Appl. Energy 147, 500-509. Ojima, J., 2003. Determining of crystalline silica in respirable dust samples by infrared spectrophotometry in the presence o f interferences. J. Occup. Health 45, 94-103. Okiongbo, K.S., Aplin, A.C., Larter, S.R., 2005. Changes in type II kerogen density as a function of maturity: Evidence from the Kimmeridge clay formation. Energ. Fuel. 19, 2495−2499. Orendt, A., Pimienta, I., Badu, S., 2013. Three-dimensional structure of the Siskin Green River Oil Shale Kerogen Model: A 28
ACCEPTED MANUSCRIPT comparison between calculated and observed properties. Energ. Fuel. 27, 702−710. Ottiger, S., Pini, R., Storti, G., Mazzotti, M., 2008. Measuring and modeling the competitive adsorption of CO 2 , CH4 , and N2 on a dry coal. Langmuir 24, 9531-9540. Paylor, A., 2015. The social–economic impact of shale gas extraction: A global perspective. Third World Q. 38, 340-355. Rexer, T.F., Mathia, E.J., Aplin, A.C., Thomas, K.M., 2014. High-pressure methane adsorption and characterization of pores in Posidonia shales and isolated kerogens. Energ. Fuel. 28, 2886-2901. Ross, D.J.K., Bustin, R.M., 2009. The importance of shale composition and pore structure upon gas storage potential of shale gas reservoirs. Mar. Pet. Geol. 26, 916-927. Ru, X., Cheng, Z., Song, L., Wang, H., Li, J., 2012. Experimental and computational studies on the average molecular structur e of
CR IP T
Chinese Huadian oil shale kerogen. J. Mol. Struct. 1030, 10-18.
Shinn, J.H., 2013. From coal to single-stage and two-stage products: A reactive model of coal structure. Fuel 63, 1187-1196. Sun, H., Mumby, S.J., Maple, J.R., Hagler, A.T., 1994. An ab initio CFF93 all-atom forcefield for polycarbonates. J. Am. Chem. Soc. 116, 2978-2987.
Sun, H., 1995. Ab initio calculations and forcefield development for computer simulation of polysilanes. Macromolecules 28,
US
701-712.
Sun, H.Y., Zhao, H., Qi, N., Li, Y., 2017. Molecular insights into the enhanced shale gas recovery by carbon dioxide in kerogen slit
AN
nanopores. J. Phys. Chem. C 121, 10233-10241.
Tissot, B.P., Welte, D.H., 1984. Petroleum Formation and Occurrence. Springer: Berlin, Germany, p720. Tong, J., Han, X., Wang S., 2011. Evaluation of structural characteristics of Huadian oil shale kerogen using direct techniques
M
(solid-state 13 C NMR, XPS, FTIR, and XRD). Energ. Fuel. 25, 4006-4013.
molecular modeling. Fuel 181, 330-339.
ED
Tong, J., Jiang, X., Han, X., Wang, X., 2016. Evaluation of the macromolecular structure of Huadian oil shale kerogen using
spectroscopy. Fuel 65, 541-546.
13
C nmr
PT
Trewhella, M.J., Poplett, I.J.F., Grint, A., 1986. Structure of Green River oil shale kerogen: Determination using solid state
Ungerer, P., Collell, J., Yiannourakou, M., 2014. Molecular modeling of the volumetric and thermodynamic properties of kerogen:
CE
Influence of organic type and maturity. Energ. Fuel. 29, 91-105. Vandenbroucke, M., 2003. Kerogen: from types to models of chemical structure. Oil Gas Sci. Technol. 58, 243−269. Vandenbroucke, M., Largeau, C., 2007. Kerogen origin, evolution and structure. Org. Geochem. 38, 719-833.
AC
Wang, J., Du, J., Chang, L., Xie, K., 2010. Study on the structure and pyrolysis characteristics of Chinese western coals. Fuel Process Technol. 91, 430-433.
Wang, X.Q., Zhai, Z.Q., Jin, X., Wu, S.T., Li, J.M., Sun, L., Liu, X.D., 2016. Molecular simulation of CO 2 /CH4 competitive adsorption in organic matter pores in shale under certain geological conditions. Petrol. Explor. Dev. 43, 841-848. Wei, Z., Gao, X., Zhang, D., Da, J., 2005. Assessment of thermal evolution of kerogen geopolymers with their structural parameters measured by solid-state 13 C NMR spectroscopy. Energ. Fuel. 19, 240-250. Weijermars, R., 2013. Economic appraisal of shale gas plays in Continental Europe. Appl. Energy 106, 100–115. Weijermars, R., 2014. US shale gas production outlook based on well roll-out rate scenarios. Appl. Energy 124, 283–297. Wu, K.L., Li, X.F., Wang, C.C., Yu, W., Chen, Z.X., 2015. Model for surface diffusion of adsorbed gas in nanopores of shale g as reservoirs. Ind. Eng. Chem. Res. 54, 3225-3236. Wu, K.L., Chen, Z.X., Li, X.F., Dong, X.H., 2016a. Methane storage in nanoporous material at supercritical temperature over a wide range of pressures. SCI REP-UK 6, 33461. 29
ACCEPTED MANUSCRIPT Wu, K.L., Li, X.F., Guo, C.H., Wang, C.C., Chen, Z.X., 2016b. A unified model for gas transfer in nanopores of shale-gas reservoirs: coupling pore diffusion and surface diffusion. SPE J. 21, 1583-1611. Wu, W., Huang, S., Hu, G., Gong, D., 2014. A comparison between shale gas and conventional gas on geochemical characteristics in Weiyuan area. Nat. Gas Geosci. 25, 1994-2002. Yang, F., Ning, Z.F., Wang, Q., Liao, H.M., 2013. Integrated study of reservoir characteristics of a shale gas reservoir: A c ase study from Sichuan Basin of China. SPE Asia Pacific Oil and Gas Conference and Exhibition, Jakarta, Indonesia, 22-24 October, SPE-165870. Yang, F., Ning, Z.F., Wang, Q., Liu, H.Q., 2016. Pore structure of Cambrian shales from the Sichuan Basin in China and implications to gas storage. Mar. Petrol. Geol. 70, 14-26.
CR IP T
Yen, T., 1979. Structural difference between petroleum and coal-derived asphaltenes. Am. Chem. Soc., Div. Pet. Chem., Prepr. 24, 901-909.
Yiannourakou, M., Ungerer, P., Leblanc, B., Rozanska, X., Saxe, P., Vidal-Gilbert, S., Gouth, F., Montel, F., 2013. Molecular simulation of adsorption in microporous materials. Oil Gas Sci. Technol. 68, 977−994.
Energy 148, 49–65.
US
Yuan, J., Luo, D., Feng, L., 2015a. A review of the technical and economic evaluation techniques for shale gas development. Appl.
Yuan, Q., Zhu, X., Lin, K., Zhao, Y.P., 2015b. Molecular dynamics simulations of the enhanced recovery of confined methane wit h
AN
carbon dioxide, Phys. Chem. Chem. Phys. 17, 31887−31893.
Zeng, Y., Wu, C., 2007. Raman and infrared spectroscopic study of kerogen treated at elevated temperatures and pressures. Fuel 86, 1192-1200.
M
Zhang, J., Liu, K., Clennell, M. B., Dewhurst, D.N., Pervukhina, M., 2015. Molecular simulation of CO 2 –CH4 competitive adsorption and induced coal swelling. Fuel 160, 309-317.
ED
Zhang, Q., Han, K., Li, S., Li, M., Li, J., Ren, K., 2017. Synthesis of garlic skin-derived 3D hierarchical porous carbon for high-performance supercapacitors. Nanoscale 10, 2427.
PT
Zhao, H., Wu, T., Firoozabadi, A., 2018. High pressure sorption of various hydrocarbons and carbon dioxide in Kimmeridge Blackstone and isolated kerogen. Fuel 224, 412-423
pores. Fuel 180, 718-726.
CE
Zhou, B., Xu, R., Jiang, P., 2016. Novel molecular simulation process design of adsorption in realistic shale kerogen spheric al
Zhu, X., Zhao, Y.P., 2014. Atomic mechanisms and equation of state of methane adsorption in carbon nanopores. J. Phys. Chem. C
AC
118, 17737-17744.
Graphical abstract Highlights: (1) Novel construction process for OM model is elaborated. (2) Direct characterization experiments are used to study kerogen chemical structure. (3) Representative molecular model of Cambrian OM is generated.
30
ACCEPTED MANUSCRIPT (4) Microstructure characterization is performed on the OM model.
AC
CE
PT
ED
M
AN
US
CR IP T
(5) Adsorption behaviors of gas components in the OM skeleton are studied.
31