Molecular Simulation of Bulk Organic Matter in Type II Shales in the

Nov 12, 2014 - A quantitative comparison to PSD determined on shales cannot be achieved, because experimental PSDs reported in the literature have bee...
0 downloads 8 Views 2MB Size
Article pubs.acs.org/EF

Molecular Simulation of Bulk Organic Matter in Type II Shales in the Middle of the Oil Formation Window Julien Collell,† Philippe Ungerer,‡ Guillaume Galliero,*,† Marianna Yiannourakou,‡ François Montel,§ and Magali Pujol§ †

Laboratoire des Fluides Complexes et leurs Réservoirs (LFC-R), UMR-5150 with CNRS and TOTAL, Université de Pau et des Pays de l’Adour, bP 1155, 64013 Pau, France ‡ Materials Design SARL, 18 rue de Saisset, 92120 Montrouge, France § Centre Scientifique et Technique Jean-Féger (CSTJF), TOTAL S.A., Avenue Larribau, 64018 Pau, France S Supporting Information *

ABSTRACT: In this work, we use molecular simulations to determine the structural and physical properties of the organic matter present in type II shales in the middle of the oil generation window. The construction of molecular models of organic matter constrained by experimental data is discussed. Using a realistic molecular model of organic matter, we generate, by molecular dynamics simulations, structures that mimic bulk organic matter under typical reservoir conditions. Consistent results on density, diffusion, and specific adsorption are found between simulated and experimental data. These structures enable us to provide information on the fluid distribution within the organic matter, the pore size distributions, the isothermal compressibility, and the dynamic of the fluids within the kerogen matrix. This study shows that a consistent description at the molecular level combined with molecular simulations can be useful, in complement of experiments, to investigate the organic matter present in shales.



III, humic) and their composition.10 The organic matter is mainly made of kerogen (generally the most abundant), asphaltenes, resins, hydrocarbons, and other fluids (such as carbon dioxide, water, and nitrogen).10 The experimental characterization of the physical and structural properties of the organic matter in reservoir conditions is either difficult or impossible. Indeed, it requires, most of the time, the extraction of the organic matter from the inorganic phase. This step is, in itself, quite challenging.11,12 This results in the production of a loosely consolidated sintered material, in which concerns remain about the textural and mechanical integrity of the samples. In this work, we propose to overcome some of these issues by the means of molecular simulations based on our previous works13,15,16 to generate molecular models of organic matter. The organic matter is considered as bulk, and the interactions with the inorganic matrix are not represented. To do so, we take into account the kerogen, asphaltenes, resins, carbon dioxide, water, and hydrocarbons. We aim to build a molecular model that mimics a type II organic matter of intermediate maturity in the middle of the oil formation window, whose hydrocarbon distribution corresponds to a wet gas or a condensate gas.17,18 Then, using molecular simulations, we produce molecular structures of bulk organic matter, under typical reservoir conditions. The species that composes the organic matter generates an intrinsic porous network, which results from the

INTRODUCTION Since the beginning of their industrial development 10 years ago, shales are a resource of growing interest to face the conventional natural gas resource rarefaction. According to the International Energy Agency (IEA), their production allowed for the division by a factor 2 of the price of natural gas in the United States since 2010.1 As a consequence, the oil and gas industry is now interested in producing shales that contain liquid-rich fluids, such as wet gases or oils, which are more valuable compared to dry gas from mature shales. Shales are multi-phase materials, made of organic nodules scattered within an inorganic matrix,2 mostly made of quartz, clays, and calcite.3,4 Although the organic nodules represent only a few percentages of the total volume in shales,5 a large part of the hydrocarbons contained in shales is stored in them. The organic matter that composes these nodules is a microporous material,6 with an important affinity for the hydrocarbons. Ross et al.4 reported the Brunauer−Emmett− Teller (BET) surface area in shale samples of approximately 30 m2/g, whereas Pribylov et al.7 reported the BET surface area on kerogen samples of 80 m2/g on average. Despite the fact that these values may strongly vary from one sample to another, these results indicate that most of the adsorption area is located in the organic matter. As a result, almost half of the total hydrocarbons is assumed to be adsorbed in the organic matter.8 Shale organic matter comes from the sedimentary organic matter maturation during the burial process.9 Dependent upon its composition and the history of the burying process, the properties of the organic matter may strongly differ. They can be sorted depending upon their type (I, algal; II, liptinic; and © 2014 American Chemical Society

Received: September 24, 2014 Revised: November 7, 2014 Published: November 12, 2014 7457

dx.doi.org/10.1021/ef5021632 | Energy Fuels 2014, 28, 7457−7466

Energy & Fuels

Article

hydrocarbons and soluble in liquid n-alkanes, while asphaltenes are insoluble in n-alkanes but soluble in polar solvent. Because of similarities in their compositions, the distinction between asphaltenes and resins is mostly based on a molecular weight basis. The typical molecular weight of resins is roughly 300 and 700 g/mol for asphaltenes in such shales.23,22 Compounds Released during Catagenesis. The carbon dioxide, water, and asphaltene/resin compounds are produced during the early stages of the maturation process of the sedimentary organic matter.9,10 At later stages of maturation, alkyl chains and cycles are generated together with a minor fraction of alkanes and aromatics, leading first to crude oils, and then a more important share of alkanes and aromatics is generated, leading to lighter gases. In proposing the values shown in Table 1, Tissot and Welte10 have assumed that the atomic balance is preserved during the maturation process. This is equivalent to consider that the compounds generated during the diagenesis and catagenesis processes have not been expelled and that carbon dioxide has been preserved from reaction with the inorganic matrix. The hydrocarbon distribution in such shales corresponds to a typical condensate gas.17,18,24 Size of the Three-Dimensional Structures. The size of the generated structures may have an impact on the properties of the structures because finite size effects may prevent them to reach consistent macroscopic values. Moreover, it is well-known that amorphous materials typically exhibit structural features ranging from a few angströms to a few nanometers.25,26 These kinds of features would only be reproduced if the molecular model is bigger than their typical length. Pikunic et al.27,28 and Jain et al.29 performed molecular simulations on reconstructed disordered nanoporous carbon structures of 25 Å long by reverse Monte Carlo simulations. The simulated and experimental nitrogen adsorption isotherms agree for the pressure range corresponding to adsorption in the ultra-microporous up to the microporous network.6 This indicates that such structures have a sufficient length to reproduce the porous network relevant at the nanometer scale. Orendt et al.30 generated molecular structures of immature Green River type I kerogen containing at least 1700 atoms. The simulated structural data on C solid-state NMR and particle distribution function were in agreement with experimental results. In previous works,13,15,16 we simulated structures of approximately 1300 atoms (more than 25 Å long), with kerogen molecular units of 100−260 carbon atoms. Qualitative and quantitative agreements were found on the density and adsorption capacity between the results from molecular simulations and experiments. On the basis of these considerations, a simulated structure of 25 Å (1300 atoms) has been considered to be sufficient to capture most of the structural features present in kerogen, at the nanometer scale. Molecular Model of Organic Matter. Using the characteristics of the kerogen for a typical type II shale in the oil formation window, detailed in the previous section, we propose in the following a generic molecular model of the organic matter in such shales, including the products of kerogen breakdown. The present section is organized as follows: First, the molecular model for the kerogen is introduced. Second, the molecular model used to represent the asphaltene/resin fraction is presented. Finally, the introduction of the hydrocarbons, carbon dioxide, and water in the model of organic matter is discussed. Molecular Model of Kerogen. Discussions on the molecular model of kerogen used in this work, with a detailed study of its volumetric and thermodynamic properties, are given in a companion paper.13 It aims to reproduce a generic type II kerogen of maturity corresponding to the middle of the oil formation window based on experimental data from Kelemen et al.14 The proposed model contains a total of 481 atoms, with 242 atoms of C, and its composition is given in Table 2. The molecular weight of MW = 3469.6 g/mol, a representation of the molecular model, is given in Figure 1. As seen in Table 2, the composition of the model is in really good agreement with the values reported by Kelemen et al.14 The strongest deviations are for the groupments involving nitrogen, whose fractions are too low in the atomic balance (5 atoms of nitrogen in the

phase equilibria of the different compounds. This enables us to reproduce structures of organic matter generated by the unexpelled gas and oil during the thermal cracking of the immature organic matter.19,20 This paper is organized as follows: In the first section, the properties of the organic matter present in type II organic matter of marine shales in the oil formation window are discussed and the molecular model of organic matter is detailed. In the second section, the molecular simulation techniques used to generate organic matter structures and determine their properties are presented. In the third section, the results of the simulations are discussed, and in the final section, some conclusions related to this work are given.



GENERATION OF ORGANIC MATTER MODELS

In this section, the development of the models of organic matter is detailed. The discussions include insights on the kerogen, asphaltenes, resins, hydrocarbons, carbon dioxide, and water. It should be noted that these discussions are limited to type II organic matter in the middle of the oil generation window and that the characteristics of the organic matter may strongly differ depending upon the organic matter type and maturity. Composition of the Organic Matter. In the following, the composition and characteristics of each fraction that composes the organic matter are detailed and discussed: the kerogen fraction, asphaltene/resin fraction, hydrocarbon fraction, and carbon dioxide/ water fraction. On the basis of the published data by Tissot and Welte,10 the composition of these fractions is given in Table 1.

Table 1. Atomic Balance of the Major Compounds of Organic Matter Present in Typical Type II Shales in the Middle of the Oil Formation Window, from Tissot and Welte10 phase

C

H

O

N

S

kerogen resin asphaltene CO2 H2O

806 33 22 5 0

876 41 23 0 102

40 1.6 1.6 10 51

18 0.4 0.4 0 0

5 0.1 0.1 0 0

Kerogen. The kerogen is defined as the fraction of sedimentary organic matter insoluble in dichloromethane.11 This corresponds to the heavier compounds found in the organic matter, made of polyaromatic units cross-linked by alkyl, naphthenic, ether, and sulfur bridges.21 The characteristics of a generic type II kerogen in the middle of the oil generation window have to be determined in terms of (1) elemental analyses (H, C, O, N, and S atoms), (2) functional groupments involved, and (3) molecular topology; kerogen is made of polyaromatic cluster,21 whose size and connectivity depend upon type and maturity. Detailed analytical data have been reported by Kelemen et al.14 by means of X-ray photoelectron spectroscopy (XPS) and solid-state 13C nuclear magnetic resonance (NMR) methods. Because kerogen is a cross-linked organic solid, its molecular weight cannot be measured reliably by experimental means. The molecular weight of current molecular models of kerogen is thus somewhat arbitrary. A molecular weight of 8000 g/mol11 has been initially proposed to provide a detailed representation of the chemical structure, but simplified models of 3000−4000 g/mol are sufficient to provide a fair account of the chemical structure and thermodynamic properties.14 Asphaltenes/Resins. The asphaltene and resin fractions correspond to the hydrocarbon compounds, made of polyaromatic units of lower weight compared to the kerogen. They include an important amount of heteroatoms (N, S, and O).10,22 Resins and asphaltenes present in shale of such type and maturity are similar to those observed in crude oils.10,22 Resins are defined as compounds more polar than 7458

dx.doi.org/10.1021/ef5021632 | Energy Fuels 2014, 28, 7457−7466

Energy & Fuels

Article

Table 2. Composition and Structural Parameters of a Generic Type II Kerogen in the Middle of the Oil Formation Windowa property composition

C groups

O groups

N groups

S groups

quantity H/C O/C N/C S/C aromatic C (%) C atoms per cluster attached aromatic C (%) O atoms in C−O per 100 C O in carboxylic groups per 100 C O in ketones (>CO) per 100 C pyrrolic (mol % of N) pyridinic (mol % of N) quaternary (mol % of N) percentage of aromatic S percentage of aliphatic S (sulfide and thiols)

target

this model

0.890 0.050 0.021 0.006 54 19 30 3.5 0.7 0.8 65 18 17 54 46

0.904 0.054 0.021 0.008 57.5 19.7 29.7 3.8 0.83 0.83 60 40 0 50 50

Figure 2. Molecular model of resin used to represent the asphaltene/ resin fraction. The formula is C26OH32, and MW = 360.5 g/mol. Atom representation: black for carbon atoms, white for hydrogen atoms, and red for oxygen atoms. given in Table 1. Two molecules of asphaltenes/resins are included in the model of organic matter. Inclusion of the Other Compounds. As assumed by Tissot and Welte,10 we suppose that the compounds produced during the maturation process are unexpelled. This scheme has also been suggested in other works20,21 and seems consistent with the high hydrocarbon content in the organic matter. Nevertheless, it cannot be generalized to all source rocks, because the hydrocarbons may be expelled during the maturation process. As a result of this assumption, the mass balance between the immature and mature organic matter is preserved. The immature organic matter is almost exclusively composed of immature kerogen,10 and the other species are neglected. Thereby, the elemental analysis of the total mature organic matter has to match the atomic balance of the immature kerogen, given in Table 2. Thus, 15 molecules of carbon dioxide and 27 molecules of water are introduced in the model of organic matter. Because of the very small amount of nitrogen and sulfur in other compounds than kerogen, they have been neglected in the atomic mass balance. Concerning the hydrocarbon distribution, we choose to use lumped fractions. Even if detailed analysis of the hydrocarbon distribution was available, the use of well-defined components would raise two problems: (1) Hundreds of individual components are present in natural gas and oil, but a detailed analysis of all of the isomers for heavy hydrocarbons is not always doable. (2) The representation of diluted compounds with a finite number of molecules would require huge systems, which is not tractable by molecular simulations. As a result, we simplified our detailed analysis to a limited number of lumped fractions. The hydrocarbon distribution generated includes 8 lumped compounds: 6 paraffins from 1−14 carbon (C) atoms and 2 representative aromatic compounds,10 toluene and dimethylnaphtalene. The content in each compound is adjusted to match the hydrocarbon distribution17,18,24,31 and the atomic mass balance of the immature kerogen, given in Table 2. The generated hydrocarbon distribution is given in Table 3. Table 4 compares the H, C, and O atomic balance of our model of organic matter and the reference immature kerogen.14 Both values are in really good agreement, and the distributions generated are consistent with results published on such shales.18 Our process allows us to build a molecular model of organic matter, which reproduces the most abundant compounds present in such shales. The organic matter molecular model contains a total of 2767 atoms. According to published data from Alberta Geological Survey (AGS),17,32,33 the carbon dioxide in such shales represents 3.2% in weight of the organic matter, which is in fair agreement with our model of organic matter. Moreover, the S1 peak from Rock-Eval analysis34 gives an estimation of the free hydrocarbon content in the organic matter. This analysis permits us to quantify the fraction of

a

The target corresponds to the reference data by Kelemen et al.14 The model corresponds to the molecular model of kerogen II C detailed in the companion paper by Ungerer et al.13

Figure 1. Molecular model of type II kerogen in the middle of the oil window. The formula is S2C242N5O13H219, and MW = 3469.6 g/mol. Atom representation: black for carbon atoms, white for hydrogen atoms, red for oxygen atoms, yellow for sulfur atoms, and purple for nitrogen atoms.

molecule) and whose quaternary nitrogen groupments are not introduced. As discussed in the previous section, a minimum of 1300 atoms is required to generate big enough representative structures. Thus, 4 molecules of kerogen are included in the molecular model of organic matter, which represents 1924 atoms for the kerogen units. Molecular Model for the Asphaltene/Resin Fraction. As seen from Table 1, in organic matter of such type and maturity, the remaining fraction for the asphaltenes and resins is small. Because of the similarities between these two fractions, we have chosen to model them with one molecule. A generic molecular model used to represent the asphaltene/resin molecules is proposed in Figure 2. The arrangement of the aromatic units is typical of the asphaltenes/resins encountered in crude oils.10,22,23 Nitrogen and sulfur groupments are neglected because their concentrations are too low. Oxygenic groupments and alkyl chains are included to match the atomic balance 7459

dx.doi.org/10.1021/ef5021632 | Energy Fuels 2014, 28, 7457−7466

Energy & Fuels

Article

interactions. Long distance corrections are used to estimate LJ energy beyond the cutoff,35,36 and the PPPM method is used for electrostatic interactions.37 The temperature is imposed during the simulations at each time step with a Nose/Hoover thermostat,41 and the pressure is imposed using a Nose/Hoover barostat.42 The time step used for the MD time integration is set to 1 fs, with a damping factor for the Nose/Hoover barostat and thermostat of 100 fs. The organic matter model is initially placed in the simulation box of 100 Å long, with the initial density being approximately 0.1 g/cm3. The first step of the MD simulations consists of a fast relaxation of the molecules at 900 K, with the volume of the simulation cell being kept constant (NVT). From the performance of succeeding NPT simulations with a stepwise decreasing temperature from 900 to 300 K, the organic matter is condensed and brought in reservoir conditions. The whole MD simulation procedure with the succeeding NVT/NPT steps is exposed in Table 5.

Table 3. Molar Composition of the the Hydrocarbon Distribution in Type II Organic Matter in the Middle of the Oil Generation Windowa

a

representative molecule

number of molecules

molecular fraction (%)

weight fraction (%)

methane ethane propane n-butane n-octane n-tetradecane toluene dimetylnaphtalene total

40 5 3 2 4 3 2 2 61

65.7 8.2 4.8 3.3 6.6 4.8 3.3 3.3 100

24.8 5.8 5.1 4.5 17.6 23.0 7.1 12.1 100

This distribution corresponds to a condensate-rich gas.18,24,31.

Table 4. Comparison of the Atomic Balance of H, C, and O between Our Molecular Model of Organic Matter and the Reference Data on Immature Type II Kerogen, Studied by Kelemen et al.14 phase

H/C

O/C

kerogen +hydrocarbons, CO2, and H2O +asphaltenes/resins =our model immature kerogen

0.91 2.53 1.23 1.18 1.17

0.054 0.294 0.039 0.091 0.097

Table 5. NVT/NPT relaxation of the molecular models of type II organic matter and kerogen in the middle of the oil formation window

linear and cyclic alkanes that contains between 1 and approximately 10 C atoms. However, this does not characterize the initial gas in place because an important fraction of lighter compounds (mostly methane) is outgassed during the extraction of the field samples. The data published by the AGS17,32,33 reports that the free hydrocarbon amount from the S1 peak represents 7.13% of the organic matter weight. The fraction of linear alkanes from 1 to 8 C, including toluene in our model, represents 9.8% of the total weight of the organic matter and 6.3% of the total weight if the methane is considered as totally outgassed. The value with 50% of the initial methane content is 8.1%, which is consistent with the data from the AGS.

ensemble

T (K)

t (ps)

P (MPa)

NVT NPT NPT NPT NPT

900 900 700 500 300

12.5−250 200 200 400 400

20 20 20 20

The length of each run has been set to be long enough to ensure the convergence of the density of the simulation cell to its equilibrium value; thus, the occurrence of fast and slow structural relaxation to generate equilibrated bulk structures of organic matter is allowed. The relative shortness of the runs is due to the fact that most of the kerogen and organic matter structural features are generated at the highest temperatures (with fast relaxation times), whereas the coldest runs only involve small volume relaxation. In addition, it should be noted that the most computationally expensive simulations are not related to the generation of the bulk structures but to the determination of their properties (isothermal compressibilities and self-diffusion coefficients), as mentioned in the following. A total of 20 different structures are generated by increasing the initial NVT stage from 12.5 to 250 ps. Results in the following correspond to averages on these 20 different structures. Equilibrium MD Simulations. Equilibrium MD simulations are used to determine the isothermal compressibility of the organic matter and the self-diffusion coefficients of the light compounds. Isothermal Compressibility. The isothermal compressibility χT of the structures of organic matter can be determined from volume fluctuations at equilibrium in the NPT ensemble43,44



MOLECULAR SIMULATIONS In this section, we detail the molecular simulations used in this study. In the first part, the molecular dynamics (MD)35,36 simulations used to generate molecular structures of organic matter are discussed. In the second part, the MD simulations conducted to determine the properties of the organic matter structures are introduced. In the third part, the Monte Carlo (MC) algorithms used to characterize the properties of the porous microporous structures present in the organic matter are discussed. MD. Structures of organic matter are generated by compressing the molecular model of organic matter, presented in the previous section, from high temperatures to typical reservoir conditions (300 K and 20 MPa). To do so, we use the MedeA-LAMMPS software.37 The allatoms pcff+ force field16 is used to represent the interatomic interactions. Dispersion/repulsion [Lennard−Jones (LJ) potential 6-9],38,39 electrostatic, and intramolecular interactions are taken into account. The Wäldman and Hagler combining rules40 are used for LJ interactions between unlike species. The force field enables the reproduction with high accuracy of the thermodynamic properties of organic compounds.16 A cutoff of 9.5 Å is selected in evaluating explicit LJ and Coulombic

χT =

⟨V 2⟩ − ⟨V ⟩2 kT ⟨V ⟩

(1)

where V is the volume of the porous structure, k is the Boltzmann constant, and ⟨···⟩ denotes the ensemble average. These coefficients have been determined from equilibrium NPT simulations over 50 ns. 7460

dx.doi.org/10.1021/ef5021632 | Energy Fuels 2014, 28, 7457−7466

Energy & Fuels

Article

Figure 3. Configuration and distribution of the rigid skeleton and fluid phase that compose the organic matter.

Self-Diffusion Coefficients. The self-diffusion coefficients Ds of each species are estimated from MD simulations via the calculation of the mean square displacement (MSD) of a group of particles with time. For a three-dimensional system, Ds is given by35,36

where nacc is the number of accepted insertions and ntot is the total number of achieved insertion tests. Pore Size Distribution (PSD). The PSDs are simulated by the method of sphere insertion proposed by Bhattacharya and Gubbins.46 It consists of picking points randomly distributed within the solid and evaluating the radius rs of the biggest sphere that contains this point without overlapping the solid, such as σ rij ≥ i + rs (5) 2

N

Ds =

⟨r (t + Δt ) − ri(t )⟩2 1 lim ∑ i 6N Δt →∞ i = 1 Δt

(2)

where ri(t + Δt) is the position of the particle i at time (t + Δt), ri(t) is the position of the particle i at time t, and N is the number of molecules of the species under study. Because the self-diffusion coefficients derive from the trajectory of the particles, they are usually determined in the microcanonical ensemble (NVE).35,36 However, because of the important size of our systems, the LJ and electrostatic interactions are cut at 9.5 Å to accelerate the simulations. Therefore, because of the slow convergence of the electrostatic interactions, the Hamiltonian of our system is not totally conserved during the simulation. It results in a constant rise of the temperature of the system during NVE simulations. To avoid this issue, kerogen molecules were weakly coupled to a Nose/Hoover thermostat, while the other compounds are not thermostated. The damping factor is set to 1000 fs. Thus, the thermostat has only a small impact on the kerogen motion, and it is supposed to negligibly modify the motion of the other hydrocarbons. Self-diffusion coefficients are determined from 100 ns runs. Characterization of the Porous Structure. The porous network properties are determined from Monte Carlo algorithms to measure the porous volume, the specific adsorption area, and the pore size distribution. Pore Volume. The porous volume has been determined using the algorithm proposed by Herrera et al.45 It consists of carrying out test insertions of a probe molecule within the simulation box based on a simple distance criteria. If the probe molecule does not overlap with the atoms from the solid, then the test insertion is accepted. The distance criteria used thus far is rij ≥

where rij is the distance between the sphere center and a particles from the solid, σi is the LJ diameter of the nearest particle from the solid, and rs is the sphere radius. This method provides an estimation of pore radius distribution. The PSD P(rs) is normalized, such as

∫0

P(rs) drs = 1

(6)

Direct comparison to experimental PSD is possible by the means of the pore diameter ds = 2rs. Specific Adsorption Area. The accessible adsorbent area is calculated using the methodology proposed by Düren et al.,47 where a probe molecule is “rolled” over the solid surface. The probe is randomly inserted around the atoms of the solid. The fraction of the probe inserted, which does not overlap with other atoms from the solid, contributes to the accessible pore surface. This method enables consistent prediction of the accessible BET surface in metal organic frameworks.47



RESULTS AND DISCUSSION In this section, some structural and physical properties of the organic matter structures generated by MD simulations are discussed. It is organized as follows: In the first part, the bulk density and isothermal compressibility of the structures are discussed. In the second part, the properties of the porous network of the rigid skeleton of the organic matter are discussed. Finally, in the third part, the self-diffusion of the fluid species present in the organic matter is discussed. Physical Properties of the Organic Matter. In this part of the discussion, we present results about the density and isothermal compressibility of our organic matter structures. The structures of organic matter have an average density of 1.26 ± 0.02 g/cm3, in agreement with the experimental data from Okiongbo et al.48 on isolated type II kerogen at the same maturity level. The weak dispersion of density of the final structures, generated from 20 uncorrelated initial structures,

σprobe + σi 2

+∞

(3)

where rij is the distance between the probe center and the center of the nearest particle i from the solid, σprobe is the LJ diameter of the probe molecule, and σi is the LJ diameter of the nearest particle i. The porosity is given by n Vporous (%) = acc ntot (4) 7461

dx.doi.org/10.1021/ef5021632 | Energy Fuels 2014, 28, 7457−7466

Energy & Fuels

Article

Table 6. Properties of the Organic Matter, Kerogen, and Bulk Fluid under Typical Reservoir Conditions (300 K and 20 MPa)a structure organic matter 2 × 2 × 2 model kerogen bulk fluid

χT (Pa−1)

density (g/cm3) 1.26 1.23 1.42 0.515

± ± ± ±

0.02 0.02 0.01 0.003

6.7 6.9 2.7 2.9

± ± ± ±

1.0 1.0 0.4 0.1

× × × ×

−10

10 10−10 10−10 10−8

porosity (%)

specific adsorption area (m2/g)

6.0 ± 1 7.3 0.15 ± 0.08

86 ± 10 93 4.3 ± 3

a

Porosities and specific adsorption area are determined using a methane probe.57 The uncertainties correspond to the average deviation determined on the 20 structures of organic matter.

sterical effects but also by the thermodynamic equilibrium of the mixture. The isothermal compressibility of the organic matter is presented in Table 6. The isothermal compressibilities of the bulk fluid phase are typical of a compressed gas, and the isothermal compressibilities of the organic matter are typical of a liquid,53 whereas pure kerogen is closer to that of a solid. These values are consistent with the bulk modulus K determined by atomic force microscopy (AFM) on organic matter by Ahmadov et al.54 (K−1 = 2−7.2 × 10−10). As stated before, because of the confinement of the fluid phase within the kerogen structure, the isothermal compressibility of the organic matter is not equal to the sum of the bulk fluid phase and pure kerogen isothermal compressibility, weighted by their fraction in the organic matter. In addition to the 20 structures of organic matter, a bigger structure of approximately 60 Å long has been generated to account for possible finite size effects, referred thereafter as the 2 × 2 × 2 model. The initial simulation cell consists of 2 × 2 × 2 models of organic matter with a total of 22 136 atoms. Because of the important computational cost required to generate such a structure (approximately 60 days on one processor) and to determine its properties, only one single structure is generated. Nevertheless, no finite size effects appear on the density of this structure (1.23 ± 002 g/cm3), which equals the density of the small organic matter structures within uncertainties. Properties of the Porous Structure. To compare the properties of the simulated porous structures to experiments, we need to remove the volatile compounds. Indeed, when performing experiments, the shale samples are priority-heated under vacuum. It results in the removal of carbon dioxide and the free hydrocarbon compounds. As seen from comparison to Rock-Eval analysis, the free hydrocarbon amount roughly corresponds to the amount of alkanes ranging from 1 to 8 C atoms. As a consequence, the structures, whose porous network properties are studied in the following, correspond to the organic matter structures in which alkanes from 1 to 8 C atoms and carbon dioxide have been removed. As a result, we assume that the remaining rigid skeleton of organic matter is composed of the kerogen, asphaltenes/resins, toluene, dimethylnaphtalene, and n-tetradecane. The average porosity and specific adsorption area are measured with the method described in the previous section. As mentioned in other works on microporous materials,47,51,55,56 the properties are dependent upon the size of the probing molecule used to estimate the properties of the porous structure. Thus, we choose to use a methane probe, represented by a single LJ site, whose LJ parameters are given by Müller et al.57 Average porosities and specific adsorption area are given in Table 6. The porosity and specific adsorption area of the two sizes of structures are equaled within uncertainties. No finite

indicates that the dense system has reached an equilibrium stage. It has been verified on five structures that the cutoff of the electrostatic interactions does not modify, within uncertainties, the density of the organic matter structures. A significant tendency of the polyaromatic units to stack is observed for structures of pure kerogen and organic matter, as seen in Figure 3. These observations are in agreement with recent molecular simulation studies13,15,16,30 and experimental work.10,12,49,50 Thereafter, the asphaltenes/resins, hydrocarbons, carbon dioxide, and water are referred as the fluid phase. Because of the inclusion of the fluid phase, the density of the organic matter drops from 1.42 ± 0.01 g/cm3 for the pure kerogen to 1.26 ± 0.02 g/cm3. The fluid phase fraction represents 24.3% of the total weight of the organic matter. The density of the bulk fluid phase has been measured by MD simulations at 0.515 ± 0.003 g/cm3. Considering the density of the bulk fluid phase and the pure kerogen, the density of the “unmixed” organic matter, composed of xf = 24.3% of bulk fluid phase and xk = 75.7% of pure kerogen, is 1.20 g/cm3. However, the density of the organic matter structures is higher. It indicates that, because of the interaction of the fluid phase with the kerogen molecules and the induced confinement, the fluid phase density is higher in the organic matter compared to the equivalent bulk fluid. Equivalently, this means that a negative excess volume (Vexcess = Vmixed − xfVf − xkVk) is found when mixing kerogen and its generated products (where Vexcess is the excess volume, Vmixed is the volume of the organic matter, Vf is the volume of the bulk fluid phase, and Vk is the volume of the pure kerogen). Panels a−c of Figure 3 show examples of fluid and kerogen distribution within the organic matter. It is interesting to note that, because of the global hydrophobic behavior of the organic matter, the water molecules form their own clusters. At the molecular level, the clustering of water molecules is well-explained by the formation of hydrogen bonds between water molecules, while they are marginal between water and kerogen. Furthermore, a separation of the hydrocarbon compounds occurs within the organic matter model. This segregation has two origins: (1) solubility/affinity (the hydrocarbons occupy different subdomains in the organic matter, with the preferential localization of light gas in nanometered pockets, whereas heavier hydrocarbons are incorporated within the kerogen molecules) and (2) adsorption selectivity inside the hydrocarbon pockets, whose composition varies close to the kerogen/hydrocarbon interface.16 It has been suggested51 that the kerogen acts as a molecular sieve52 during experiments; dependent upon the size of the probing molecule, some pores are not available, because they are smaller than the probing molecules. As seen from Figure 3, the scheme here is in fact more complicated because the affinity implies that the distribution of the fluid is not only driven by 7462

dx.doi.org/10.1021/ef5021632 | Energy Fuels 2014, 28, 7457−7466

Energy & Fuels

Article

unconnected network, in complement to experiments. It seems to confirm the presence of a porous network in the organic matter generated by the unexpelled oils and gas, as suggested in the literature.19,20,63 Diffusion of Light Compounds in the Organic Matter. The self-diffusion coefficients of light n-alkanes (methane to octane) and carbon dioxide, referred thereafter as light compounds, have been measured in the bulk fluid and in organic matter. A third case corresponds to the organic matter in which half of the light compounds have been removed (noted 1/2 HC). Then, a new NPT stage at 300 K and 20 MPa is achieved to relax the structure. This case is studied to investigate the evolution of the diffusion coefficients of the remaining compounds when fluids contained in the organic matter have been expelled, for example, during the recovery process. A total of 20 additional structures of organic matter with half of the light compounds have been generated, following the NVT/NPT process given in Table 5. The self-diffusion coefficients of the light compounds are given in Table 7. The self-diffusion coefficients (Ds) in the organic matter are in the order of 10−9 m2/s, which corresponds to self-diffusion coefficients in a liquid.64 These coefficients are comparable to the effective diffusivities of hydrocarbons in the organic matter reported by Thomas and Clouse65,66 and NMR experimental results.67 Theses values are also consistent with self-diffusion coefficients of gas in bulk semi-crystalline polymers.68−71 Contrary to static properties, which have low uncertainties because of a weak dispersion of the results from one structure of organic matter to another, the dynamic properties heave greater uncertainties. These are due to the strong variation of the topology and connectivity of the porous network from one structure to another. As a results, despite absolute average deviations on the density of the organic matter being lower than 2%, the absolute average deviation for self-diffusion coefficients is 29% for methane and even greater for the most diluted compounds, as butane. Nevertheless, the average deviation determined on each single structure is from 1 to 2 orders of magnitude lower (for example, the absolute average deviation of the methane self-diffusion coefficient on each structure is 1.3%). Thus, a quantitative comparison of methane self-diffusion can be achieved between the three cases presented here within a few percentages of uncertainty. In comparison to the bulk fluid, the diffusion of light compounds in the organic matter is slowed by an order of magnitude.67 Moreover, the decrease of the self-diffusion of nalkanes (proportional to 1/M0.5 in the gas phase64,72) with their molecular weight observed for bulk fluids, is flattened in the organic matter and almost insignificant for the structures of organic matter with only half of the light compounds. As the confinement of the light compounds increases (from the organic matter structures to the 1/2 HC structures), the diffusion of the light compounds is slowed, by roughly a third

size effects are noticeable on density, porosity, or specific adsorption area. The values of the specific adsorption area are consistent with experimental results reported on organic matter.7 The PSDs of the rigid skeleton of organic matter are given in Figure 4. If the PSDs measured on a single structure are

Figure 4. PSD of the molecular structures of organic matter. Black is for the PSD of the 19th structure of organic matter; red is for the average PSD on the 20 structures of organic matter; and blue is for the PSD of the big structure of approximately 60 Å long.

peculiar, the average PSD that is determined on the 20 structures exhibits a smoother shape. The average pore radius is 3.6 Å, and therefore, the pore diameter is 7.2 Å. This is approximately 1.9 times the methane molecule.57 The average pore size increases to 4.5 Å for the 2 × 2 × 2 model, which may indicate the occurrence of weak finite size effects, hardly seen on other properties. According to the International Union of Pure and Applied Chemistry (IUPAC)6 classification, the PSDs range from ultramicropores up to micropores. A quantitative comparison to PSD determined on shales cannot be achieved, because experimental PSDs reported in the literature have been determined on complete shale samples, which include both organic and inorganic phases. Nevertheless, this range corresponds to the most confined pores that have been observed on complete shale samples.58−60 Such ultra-microporous and microporous structures are similar to those obtained in disordered nanoporous carbon27,28 or in swollen bituminous coals.61,62 Clarkson et al.59 reported PSD of type II shales obtained from nitrogen and carbon dioxide adsorption analysis. The PSDs exhibit unimodal or multimodal distributions, with their lowest peak corresponding to a pore radius below 10 Å. Even if the porous volume corresponding to these pores is very low, most of the adsorption surface seems to arise from these ultramicropores and micropores.59 Thus, our models of organic matter yield microporous and ultra-microporous networks. This provides information on this

Table 7. Self-Diffusion Coefficients of Light Compounds in the Bulk Fluid and in the Organic Mattera

a

Ds (×10−9, m2/s)

C1

C2

C3

C4

C8

CO2

organic matter 1 /2 HC bulk fluid

4.8 ± 1.4 3.2 ± 1.8 24

3.3 ± 2.0 3.3 ± 2.0 23

2.4 ± 1.4 2.2 ± 1.6 16

3.2 ± 2.2 3.0 ± 2.3 18

2.5 ± 1.4 2.4 ± 1.7 5.0

2.4 ± 1.4 2.6 ± 1.6 6.0

The uncertainties correspond to the average deviation of the self-diffusion coefficients determined on the 20 structures of organic matter. 7463

dx.doi.org/10.1021/ef5021632 | Energy Fuels 2014, 28, 7457−7466

Energy & Fuels

Article

diffusion is slowed by an order of magnitude. A weak dependency of the diffusion coefficients to the weight of the diffusing compounds, typical of strong interactions of this compound with the porous skeleton, has been observed. This study shows that, using a consistent description at the molecular level and molecular simulations, it is possible to obtain insights on the bulk organic present in shales, whose experiments are quite challenging. Molecular simulations can, hence, be used as an alternative or in complement to experiments.

for methane molecules. Such slowing effects have been reported on mixture diffusion in zeolites and permeation membranes by Krishna et al.,73 where the motion of the fastest molecules is slowed because of the slowest molecules. In our case, the slowest molecules are the kerogen molecules, which are also the heaviest molecules, which can “damp” the diffusion of the lightest species.74 Another specificity of our system is that the “solid” skeleton (kerogen) is not considered as rigid and its porous network is not connected. As a result, the motion of the light compounds may be affected by the thermal motion of the kerogen molecules and the displacement of the fluid pockets during the simulation.



ASSOCIATED CONTENT

S Supporting Information *



Molecular models of type II kerogen. This material is available free of charge via the Internet at http://pubs.acs.org.

CONCLUSION In this work, on the basis of both experimental10,14 and simulation13,16,30 work, we developed a molecular model of generic type II organic matter in the middle of the oil formation window. Our model includes the most relevant phases present in the organic matter: kerogen, asphaltenes/resins, carbon dioxide, water, and hydrocarbons. From this approach, the schematic composition of the hydrocarbons and heavy products formed during the maturation process is in agreement with (1) the composition of the oils/gases and (2) the atomic balance of the immature organic matter. Furthermore, the content in free hydrocarbons and carbon dioxide has been validated with available experimental data.32,33 Using this realistic molecular model of organic matter, we generated structures by MD simulations, which mimic the bulk organic matter under typical reservoir conditions. The density of our structures of organic matter is 1.26 ± 0.02 g/cm3, in agreement with experimental data48 on such kerogen at the same maturity level. Because the density is related to the stacking of the polyaromatic cluster,10,12,50 our molecular model allows for a consistent modeling of the structuration of the organic matter. The specific adsorption area determined from molecular simulations is measured at 86 ± 10 m2/g. It is consistent with the BET specific adsorption area determined on such shales7 from adsorption isotherms. The important adsorption area compared to the inorganic matrix4 confirms that most of the hydrocarbons are adsorbed in the nodules of organic matter.8 The isothermal compressibility is consistent with the bulk modulus measurements performed by AFM54 and greater than typical isothermal compressibility of the inorganic matrix by an order of magnitude. In our structures, the porous network is generated by the unexpelled oils and gases during the thermal cracking.19,20 The water molecules are gathered in small clusters, indicating that the hydrophobic behavior is predicted as a consequence of its chemical structure, as in overmature kerogen.16 The spatial distribution of the hydrocarbons exhibits a local phase split between the heavy and light compounds, with a preferential localization of the heavy hydrocarbons within the kerogen molecules, whereas the light compounds generate nanometered pockets of gases. This is confirmed by the PSD, which indicates a porous network corresponding to ultra-micropores and micropores, comparable to those observed in nanoporous carbons27,28 and bituminous coals.61,62 The self-diffusion coefficients of light compounds are in agreement with permeability measurements65,66 and NMR experiments.67 In comparison to the bulk fluid, the self-



AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.



ACKNOWLEDGMENTS The authors gratefully acknowledge TOTAL S.A. for financial support and a grant for Julien Collell. Computer time for this study was provided by the computing facilities Mésocentre de Calcul Intensif Aquitain (MCIA) of the Université de Bordeaux and the Université de Pau et des Pays de l’Adour.



REFERENCES

(1) U.S. International Energy Agency (IEA). World Energy Outlook 2013; IEA: Washington, D.C., 2013. (2) Curtis, M. E.; Ambrose, R. J.; Sondergeld, C. H. Investigation of the relationship between organic porosity and thermal maturity in the Marcellus Shale. Proceedings of the North American Unconventional Gas Conference and Exhibition; The Woodlands, TX, June 14−16, 2012; SPE Paper 144370. (3) Gasparik, M.; Ghanizadeh, A.; Bertier, P.; Gensterblum, Y.; Bouw, S.; Krooss, B. M. High-pressure methane sorption isotherms of black shales from the Netherlands. Energy Fuels 2012, 26, 4995−5004. (4) Ross, D. J. K.; Bustin, R. M. The importance of shale composition and pore structure upon gas storage potential of shale gas reservoirs. Mar. Pet. Geol. 2009, 26, 916−927. (5) Curtis, M. E.; Ambrose, R. J.; Sondergeld, C. H. Structural characterization of gas shales on the micro- and nano-scales. Proceedings of the Canadian Unconventional Resources and International Petroleum Conference; Calgary, Alberta, Canada, Oct 19−21, 2010; SPE Paper 137693. (6) International Union of Pure and Applied Chemistry (IUPAC). Compendium of Chemical TherminologyGold Book; IUPAC: Zürich, Switzerland, 2012. (7) Pribylov, A. A.; Skibitskaya, N. A.; Zekel’, L. A. Sorption of methane, ethane, propane, butane, carbon dioxide, and nitrogen on kerogen. Russ. J. Phys. Chem. A 2014, 88, 1028−1036. (8) Wang, F. P.; Reed, R. M. Pore networks and fluid flow in gas shales. Proceedings of the SPE Annual Technical Conference and Exhibition; New Orleans, LA, Oct 4−7, 2009; SPE Paper 124253. (9) Durand, B. Kerogen: Insoluble Organic Matter from Sedimentary Rocks; Technip: Paris, France, 1980; p 519. (10) Tissot, B. P.; Welte, D. H. Petroleum Formation and Occurrence; Springer: Berlin, Germany, 1984; p 720. (11) Vandenbroucke, M. Kerogen: From types to models of chemical structure. Oil Gas Sci. Technol. 2003, 58, 243−269. (12) Vandenbroucke, M.; Largeau, C. Kerogen origin, evolution and structure. Org. Geochem. 2007, 38, 719−833.

7464

dx.doi.org/10.1021/ef5021632 | Energy Fuels 2014, 28, 7457−7466

Energy & Fuels

Article

(13) Ungerer, P.; Collell, J.; Yiannourakou, M. Molecular modeling of the volumetric and thermodynamic properties of kerogen: Influence of organic type and maturity. Energy Fuels 2015. (14) Kelemen, S. R.; Afeworki, M.; Gorbaty, M. L.; Sansone, M.; Kwiatek, P. J.; Walters, C. C.; Freund, H.; Siskin, M.; Bence, A. E.; Curry, D. J.; Solum, M.; Pugmire, R. J.; Vandenbroucke, M.; Leblond, M.; Behar, F. Direct characterization of kerogen by X-ray and solidstate 13C nuclear magnetic resonance methods. Energy Fuels 2007, 21, 1548−1561. (15) Collell, J.; Galliero, G.; Gouth, F.; Montel, F.; Pujol, M.; Ungerer, P.; Yiannourakou, M. Molecular simulation and modelisation of methane/ethane mixtures adsorption onto a microporous molecular model of kerogen under typical reservoir conditions. Microporous Mesoporous Mater. 2014, 197, 194−203. (16) Yiannourakou, M.; Ungerer, P.; Leblanc, B.; Rozanska, X.; Saxe, P.; Vidal-Gilbert, S.; Gouth, F.; Montel, F. Molecular simulation of adsorption in microporous materials. Oil Gas Sci. Technol. 2013, 68, 977−994. (17) Rokosh, C. D.; Lyster, S.; Anderson, S. D. A.; Beaton, A. P.; Berhane, H.; Brazzoni, T.; Chen, D.; Cheng, Y.; Mack, T.; Pana, C.; Pawlowicz, J. G. Open File Report: Summary of Alberta’s Shale- and Siltstone-Hosted Hydrocarbon Resource Potential; Alberta Geological Survey (AGS)/Energy Resources Conservation Board (ERCB): Edmonton, Alberta, Canada, 2012; p 339. (18) Smith Low, W.; Kennedy, J.; Fisher, G.; Greene, M. The Duvernay ShaleThe New Millennium Gold is Condensate; BMO Capital Makets: Calgary, Alberta, Canada, 2013; p 18. (19) Bernard, S.; Wirth, R.; Schreiber, A.; Schulz, H.-M.; Horsfield, B. Formation of nanoporous pyrobitumen residues during maturation of the Barnett Shale (Fort Worth Basin). Int. J. Coal Geol. 2012, 103, 3− 11. (20) Loucks, R. G.; Ruppel, S. C. Mississippian Barnett Shale: Lithofacies and depositional setting of a deep-water shale-gas succession in the Fort Worth Basin, Texas. Am. Assoc. Pet. Geol. Bull. 2007, 91, 57−601. (21) Behar, F.; Vandenbroucke, M. Représentation chimique de la structure des kérogènes et des asphaltènes en fonction de leur origine et de leur degré d’évolution. Oil Gas Sci. Technol. 1986, 41, 173−188. (22) Castex, H. Résines et asphaltènes: Évolution en fonction des types de matière organique et de leur enfouissement. Oil Gas Sci. Technol. 1985, 40, 169−189. (23) Mullins, O. C.; Sheu, E. Y.; Hammami, A.; Marshall, A. G. Asphaltenes, Heavy Oils, and Petroleomics; Springer: New York, 2007; p 677. (24) Lagache, M.; Ungerer, P.; Boutin, A. Prediction of thermodynamic derivative properties of natural condensate gases at high pressure by Monte Carlo simulation. Fluid Phase Equilib. 2004, 220, 211−223. (25) Palmer, J. C.; Gubbins, K. E. Atomistic models for disordered nanoporous carbons using reactive force fields. Microporous Mesoporous Mater. 2012, 154, 24−37. (26) Salazar, R.; Gelb, L. D. A computational study of the reconstruction of amorphous mesoporous materials from gas adsorption isotherms and structure factors via evolutionary optimization. Langmuir 2007, 23, 530−541. (27) Pikunic, J.; Clinard, C.; Cohaut, N.; Gubbins, K. Structural modeling of porous carbons: constrained reverse Monte Carlo method. Langmuir 2003, 19, 8565−8582. (28) Pikunic, J.; Llewellyn, P.; Pellenq, R.; Gubbins, K. Argon and nitrogen adsorption in disordered nanoporous carbons: Simulation and experiment. Langmuir 2005, 21, 4431−4440. (29) Jain, S.; Pellenq, R.; Pikunic, J.; Gubbins, K. Molecular modeling of porous carbons using the hybrid reverse Monte Carlo method. Langmuir 2005, 22, 9942−9948. (30) Orendt, A.; Pimienta, I.; Badu, S. Three-dimensional structure of the Siskin Green River oil shale kerogen model: A comparison between calculated and observed properties. Energy Fuels 2013, 27, 702−710.

(31) Arnaud, J. F. Caractérisation des propriétés physiques et thermodynamiques des fluides pétroliers à haute pression. Ph.D. Thesis, Institut Francais du Pétrole, Rueil-Malmaison, France, 1995; p 109. (32) Beaton, A. P.; Pawlowicz, J. G.; Anderson, S. D. A.; Berhane, H.; Rokosh, C. D. Rock Eval, Total Organic Carbon and Adsorption Isotherms of the Duvernay and Muskwa Formations in Alberta: Shale Gas Data Release; Alberta Geological Survey (AGS): Edmonton, Alberta, Canada, 2010; p 39. (33) Beaton, A. P.; Pawlowicz, J. G.; Anderson, S. D. A.; Berhane, H.; Rokosh, C. D. Mineralogy, Permeametry, Mercury Porosimetry, Pycnometry and Scanning Electron Microscope Imaging of Duvernay and Muskwa Formations in Alberta: Shale Gas Data Release; Alberta Geological Survey (AGS): Edmonton, Alberta, Canada, 2010; p 61. (34) Lafargue, E.; Marquis, F.; Pillot, D. Rock-Eval 6 applications in hydrocarbon exploration, production, and soil contamination studies. Oil Gas Sci. Technol. 1998, 53, 421−437. (35) Frenkel, D.; Smit, B. Understanding Molecular SimulationFrom Algorithms to Applications; Academic Press: San Diego, CA, 1996; p 658. (36) Allen, M.; Tildesley, D. J. Computer Simulations of Liquids, Computer Simulations of Liquids; Oxford University Press: Oxford, U.K., 1989; p 400. (37) Plimpton, S. Fast parallel algorithms for short-range molecular dynamics. J. Comput. Phys. 1995, 117, 1−19. (38) Jones, J. E. On the determination of molecular fields. II. From the equation of state of a gas. Proc. R. Soc. London 1924, 106, 463−477. (39) Sun, H. COMPASS: An ab initio force-field optimized for condensed-phase applications overview with details on alkane and benzene compounds. J. Phys. Chem. B 1998, 5647, 7338−7364. (40) Waldman, M.; Hagler, A. T. New combining rules for rare gas van der waals parameters. J. Comput. Chem. 1993, 14, 1077−1084. (41) Hoover, W. G. Canonical dynamics: Equilibrium phase-space distributions. Phys. Rev. A: At., Mol., Opt. Phys. 1985, 31, 1695−1697. (42) Hoover, W. G. Constant-pressure equations of motion. Phys. Rev. A: At., Mol., Opt. Phys. 1986, 34, 2499−2500. (43) Ungerer, P.; Tavitian, B.; Boutin, A. Applications of Molecular Simulation in the Oil and Gas Industry: Monte Carlo Methods; Editions Technip: Paris, France, 2005; p 295. (44) Lagache, M.; Ungerer, P.; Boutin, A.; Fuchs, A. H. Prediction of thermodynamic derivative properties of fluids by Monte Carlo simulation. Phys. Chem. Chem. Phys. 2001, 3, 4333−4339. (45) Herrera, L.; Do, D. D.; Nicholson, D. A. Monte Carlo integration method to determine accessible volume, accessible surface area and its fractal dimension. J. Colloid Interface Sci. 2010, 348, 529− 536. (46) Bhattacharya, S.; Gubbins, K. E. Fast method for computing pore size distributions of model materials. Langmuir 2006, 22, 7726− 7731. (47) Düren, T.; Millange, F.; Férey, G. Calculating geometric surface areas as a characterization tool for metal−organic frameworks. J. Phys. Chem. C 2007, 111, 15350−15356. (48) Okiongbo, K. S.; Aplin, A. C.; Larter, S. R. Changes in type II kerogen density as a function of maturity: Evidence from the Kimmeridge Clay Formation. Energy Fuels 2005, 19, 2495−2499. (49) Zhang, L.; Greenfield, M. L. Molecular orientation in model asphalts using molecular simulation. Energy Fuels 2007, 21, 1102− 1111. (50) Oberlin, A.; Boulmier, J. L.; Durand, B. Electron microscope investigation of the structure of naturally and artificially metamorphosed kerogen. Geochim. Cosmochim. Acta 1974, 38, 647−650. (51) Kang, S.; Fathi, E.; Ambrose, R.; Akkutlu, I.; Sigal, R. CO2 storage capacity of organic-rich shales. Proceedings of the SPE Annual Technical Conference and Exhibition; Florence, Italy, Sept 19−22, 2010; SPE Paper 134583. (52) Kärger, J.; Ruthven, D. M.; Theodorou, D. N. Diffusion in Nanoporous Materials; Wiley: Hoboken, NJ, 2012; p 872. (53) Motakabbir, K.; Berkowitz, M. Isothermal compressibility of SPC/E water. J. Phys. Chem. 1990, 94, 8359−8362. 7465

dx.doi.org/10.1021/ef5021632 | Energy Fuels 2014, 28, 7457−7466

Energy & Fuels

Article

(54) Ahmadov, R.; Vanorio, T.; Mavko, G. Confocal laser scanning and atomic-force microscopy in estimation of elastic properties of the organic-rich Bazhenov Formation. Leading Edge 2009, 28, 18−23. (55) Navarrete, R.; Llewellyn, P.; Rouquerol, F.; Denoyel, R.; Rouquerol, J. Calorimetry by immersion into liquid nitrogen and liquid argon: A better way to determine the internal surface area of micropores. J. Colloid Interface Sci. 2004, 277, 383−386. (56) Llewellyn, P.; Rodriquez-Reinoso, F.; Rouquerol, J.; Seaton, N. A. Characterization of Porous Solids VII; Elsevier: Amsterdam, Netherlands, 2007; pp 49−56. (57) Müller, D.; Ó przynski, J.; Müller, A.; Fischer, J. Prediction of thermodynamic properties of fluid mixtures by molecular dynamics simulations: Methane−ethane. Mol. Phys. 1992, 75, 363−378. (58) Bustin, R. M.; Bustin, A. M. M.; Cui, A.; Ross, D.; Murthy Pathi, V. Impact of shale properties on pore structure and storage characteristics. Proceedings of the SPE Shale Gas Production Conference; Fort Worth, TX, Nov 16−18, 2008; SPE Paper 119892. (59) Clarkson, C. R.; Solano, N.; Bustin, R. M.; Bustin, A. M. M.; Chalmers, G. R. L.; He, L.; Melnichenko, Y. B.; Radlinski, A. P.; Blach, T. P. Pore structure characterization of North American shale gas reservoirs using USANS/SANS, gas adsorption, and mercury intrusion. Fuel 2013, 103, 606−616. (60) Jun-yi, L.; Zheng-song, Q.; Wei-an, H. Nano-pore structure characterization of shales using gas adsorption and mercury intrusion techniques. J. Chem. Pharm. Res. 2014, 6, 850−857. (61) Van Niekerk, D.; Mathews, J. P. Molecular dynamic simulation of coal−solvent interactions in Permian-aged South African coals. Fuel Process. Technol. 2011, 92, 729−734. (62) Mathews, J. P.; van Duin, A. C. T.; Chaffee, A. L. The utility of coal molecular models. Fuel Process. Technol. 2011, 92, 718−728. (63) Ungerer, P.; Behar, E.; Discamps, D. Tentative calculation of the overall volume expansion of organic matter during hydrocarbon genesis from geochemistry data. Implications for primary migration. Adv. Org. Geochem. 1981, 30, 1105−1125. (64) Poling, B.; Prausnitz, J.; Paul, O. J.; Reid, R. The Properties of Gases and Liquids; McGraw-Hill: New York, 2001; p 803. (65) Thomas, M. M.; Clouse, J. A. Primary migration by diffusion through kerogen: I. Model experiments with organic-coated rocks. Geochim. Cosmochim. Acta 1990, 54, 2775−2779. (66) Thomas, M. M.; Clouse, J. A. Primary migration by diffusion through kerogen: II. Hydrocarbon diffusivities in kerogen. Geochim. Cosmochim. Acta 1990, 54, 2781−2792. (67) Viswanathan, R. K. K.; Minh, C. C.; Zielinski, L.; Vissapragada, B.; Akkurt, R.; Song, Y.-Q.; Liu, C.; Jones, S.; Blair, E. Characterization of gas dynamics in kerogen nanopores by NMR. Proceedings of the SPE Annual Technical Conference and Exhibition; Denver, CO, Oct 30−Nov 2, 2011; SPE Paper 147198. (68) Kucukpinar, E.; Doruker, P. Molecular simulations of small gas diffusion and solubility in copolymers of styrene. Polymer 2003, 44, 3607−3620. (69) Kucukpinar, E.; Doruker, P. Molecular simulations of gas transport in nitrile rubber and styrene butadiene rubber. Polymer 2006, 47, 7835−7845. (70) Memari, P.; Lachet, V.; Rousseau, B. Gas permeation in semicrystalline polyethylene as studied by molecular simulation and elastic model. Oil Gas Sci. Technol. 2013, 1965, 9. (71) Gusev, A. A.; Müller-Plathe, F.; van Gunsteren, W. F.; Suter, U. W. Dynamics of Small Molecules in Bulk Polymers; Springer: Berlin, Germany, 1994; pp 207−247. (72) Yu, Y.; Gao, G. Lennard−Jones chain model for self-diffusion of n-alkanes. Int. J. Thermophys. 2000, 21, 14. (73) Krishna, R.; van Baten, J. M. Mutual slowing-down effects in mixture diffusion in zeolites. J. Phys. Chem. C 2010, 114, 13154− 13156. (74) Facelli, J. C.; Pugmire, R. J.; Pimienta, I. S.; Badu, S.; Orendt, A. M. Atomistic Modeling of Oil Shale Kerogens and Asphaltenes along with Their Interactions with the Inorganic Mineral Matrix; Institute of Clean and Secure Energy, University of Utah: Salt Lake City, UT, 2011; Report DOE-FE0001243, p 39. 7466

dx.doi.org/10.1021/ef5021632 | Energy Fuels 2014, 28, 7457−7466