Kerogen Swelling and Confinement: Its implication on Fluid Thermodynamic Properties in Shales

Type I kerogen was isolated from Green River Shale and characterized using SEM, TGA, DSC, and nitrogen adsorption. The swelling behavior of this kerogen with decane was analyzed using traditional test-tube swelling experiment and Dynamic Light Scattering. The TGA and DSC were used to analyze the thermodynamic behavior of decane that was sorbed in the kerogen and show that kerogen suppresses the boiling point of decane due to the effect of confinement. However, the suppression is larger when oil (a multicomponent mixture) was used, possibly due to the combined effect of differential uptake of components by kerogen (kerogen prefers and sorbs polars and aromatics more than saturates, leading to splitting of oil into a sorbed and a free phase) and confinement in nano pores. Test-tube swelling, TGA, and DSC experiments were also performed on pyridine(polar-aromatic)-swelled kerogen. The combined and individual contributions from the two effects (the effect of confinement and differential uptake of hydrocarbon components) on properties of liquid in contact with kerogen, are studied in this work. Molecular Dynamics (MD) simulations revealed the variation in the swelling of type II kerogen in the presence of same amount of different liquids (differential swelling of kerogen).

The simulation matrix was selected in order to perform simulation of swelling of kerogen with fixed masses of a series of hydrocarbon liquids that ranged from being polar, non-polar, saturates, and aromatic compounds. Thus, a model of type II kerogen was simulated with 17 different liquids in multiple MD simulations. Each simulation consists of a fixed mass of kerogen simulated with a liquid (mass of liquid equals 20% of mass of the kerogen). The kerogen and liquid molecules were randomly inserted in the simulation box and the system was annealed together using the 'Simulated Annealing' procedure discussed later. Initial configuration consists of the liquid molecules placed outside of the kerogen in a simulation box. The liquid molecules diffused into the kerogen and were sorbed by it during the subsequent MD simulations. The properties of the final system at the end of a successful simulation annealing procedure are not a function of the time when the liquid molecules are inserted in the simulation box.

Workflow for Simulated Annealing and Molecular Dynamics Simulations
A method has been developed to construct a molecular structure of bulk kerogen from a single molecule. A simulated annealing procedure was performed on a group of several kerogen molecules using the GAFF force field, which was originally developed for molecular dynamics simulations of organic molecules 1 . In the current work, the application of GAFF was extended from simpler organic molecules to complex organic macromolecules like kerogen. The type II C kerogen model was chosen from Ungerer et al. (2015) 2 . The chemical formula of Model II C is C242H219O13N5S2 with 481 atoms in one molecule. The H/C and O/C ratios of this model unit are 0.905 and 0.054 respectively, and its location on a Van Krevelen diagram, near the end of the oil window is shown in Figure 1.
The force field parameters for GAFF were chosen from Wang et al, 2004 1 and the 'gaff.lt' file available in the Moletemplate 3 software distribution. The Gasteiger-Marsili 4 charges were assigned to all the atoms using Avogadro molecular modeling tool 5 .
In order to perform MD simulations, a workflow was developed comprising of series of steps as shown in figure S1. The steps are discussed below and in figure S1: 1. First, a molecule of Type II C kerogen was converted from (x,y,z) format to protein data bank (PDB) format. A single kerogen model II C is shown in figure 10. 2. The atomtype program of Ambertools 6 was used to convert the file into GAFF friendly PDB format 3. The Packmol (Martinez L, 2009) tool was used to randomly stuff 35 of these kerogen molecule models into a simulation box 4. The Visual Molecular Dynamics, VMD 7 tool was used to produce an input file for LAMMPS 8 5. The bond, angles and dihedral angle were manually entered for all the atoms of kerogen from gaff.lt 6. An instruction file for LAMMPS, which included the series of simulated annealing steps was written 7. The simulations for simulated annealing, followed by production runs were run in the LAMMPS in parallel with multiple processors of a super computer 8. The density of kerogen and probability distribution function of the C-C bond lengths were calculated at the end of simulated annealing (after relaxation steps). 9. Swelling simulations started with first annealing the bulk kerogen with organic liquid (of weight equal to 20% of that of kerogen) till density becomes constant. VMD was used to visualize the final relaxed structure of pure kerogen and kerogen with liquid ( Figure 15) at the end of MD simulation following the annealing steps. was used to calculate the non-bonded molecular forces. Here, is the depth of the potential well, σ is the finite distance at which the inter-particle potential is zero, r is the distance between the particles. The Lorentz-Berthelot classical mixing rule 10,11 was used to calculate the parameters for unlike atoms in the L-J potential. Here, q1 and q2 are the charges on interacting atoms, ∈ is the dielectric constant of the medium and C is a constant.

Simulated Annealing Procedure
In order to relax the system, the bulk kerogen structure made up of fixed number of kerogen molecules systematically went through series of energy minimization steps. The initial density of kerogen in the 100Å cube was 0.206 g/cc. First, the structure was simulated in a canonical (NVT) ensemble as the temperature was increased to 2000K, during a total time of 0.2 ns (in 2 million time steps) with density of 0.22 g/cc. The system was then simulated in a series of isothermal-isobaric (NPT) ensembles for a total time of 1.4 ns, consisting of 5 smaller steps of 0.28 ns each with step-wise reduction in temperatures and pressures. Each of these steps consist of MD simulations in NPT ensemble with pressure, temperature in each step as (1000K, 1000atm), (1000K-2000K, 1000atm), (500K, 500atm), (400K, 300atm), (300K, 300atm) respectively. These temperatures and pressures during each step is indicated in Figure 13 of the manuscript. The final pressure and temperature was close to reservoir pressure and temperature in order to replicate the realistic in-situ conditions for kerogen. A schematic of all the annealing steps are given in Figure S1. For any of the NPT or NVT cycle, the time steps varied from 0.0001 to 0.01 ns. The density of the final annealed structure of bulk kerogen at minimum energy state was calculated at the end of the simulated annealing procedure.
A Nose-Hoover temperature thermostat and Nose-Hoover pressure barostat were used in LAMMPS molecular simulator to generate positions and velocities sampled from the canonical (NVT) and isothermal-isobaric (NPT) ensembles. Since the pressure includes a kinetic component due to particle velocities, both thermostatting and barostatting in MD simulations require calculation of the temperature. A target temperature (T) and/or pressure (P) mentioned in former paragraph was specified in each run, and the thermostat or barostat attempted to equilibrate the system to the requested T and/or P. The thermostatting and barostatting is achieved by adding some dynamic variables which are coupled to the particle velocities (thermostatting) and simulation domain dimensions (barostatting).
The kerogen-liquid simulations were carried by randomly inserting kerogen and liquid molecules in a 10nm simulation cubic-box. Simulated annealing procedure discussed above was used to first minimize the energy of the system. MD simulations were then performed in NPT ensemble to measure the volume of the system after quasiequilibrium was reached. The initial configuration of kerogen would have some effect on density of kerogen 12,13 and therefore, on swelling ratios (Qv) calculated by MD simulations. However, this work focused more on the comparison of the Qv for each liquid than the absolute value of Qv for each liquid, i.e. the interest was in obtaining a trend in the swelling ratios of kerogen with different liquids.