Reactive Molecular Dynamics Simulation of Kerogen Thermal

Subscriber access provided by LAURENTIAN UNIV

Article

Reactive Molecular Dynamics Simulation of Kerogen Thermal Maturation and Crosslinking Pathways Gorakh Pawar, Paul Meakin, and Hai Huang Energy Fuels, Just Accepted Manuscript • DOI: 10.1021/acs.energyfuels.7b01555 • Publication Date (Web): 27 Sep 2017 Downloaded from http://pubs.acs.org on October 9, 2017

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

Energy & Fuels is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 47

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Energy & Fuels

Reactive Molecular Dynamics Simulation of Kerogen Thermal Maturation and Crosslinking Pathways GORAKH PAWAR1, PAUL MEAKIN2 and HAI HUANG1* 1

Idaho National Laboratory, 2525 Fremont Avenue, Idaho Falls, ID, United States 83402

2

Temple University, 1801 N Broad St, Philadelphia, PA, United States, 19122 * Corresponding author - Email: [email protected]; Tel: (208)526-4895

ACS Paragon Plus Environment

1

Energy & Fuels

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 47

ABSTRACT Molecular dynamics simulations were performed with a ReaxFF reactive force field to investigate bond breaking and bond formation mechanisms during the thermal maturation of kerogen and potential crosslinking pathways towards the formation of three-dimensional (3D) quasi-infinite molecular networks (crosslinked kerogen macromolecules). Starting with small ensembles of high molecular mass models for immature type I Green River Shale kerogen (kerogen 1-I), top of the oil window type II kerogen (kerogen 2-L), and low maturity type III kerogen (kerogen 3-L), low molecular mass species including H2O, C2H4, C3H6 were produced as the maturities of the remaining kerogens increased. Highly reactive fragments, which are not detected in pyrolysis experiments, were also produced. Further, the crosslinking mechanisms in the newly develop polymeric kerogen networks appear to be highly complex, and covalent —C—S—, —C—O— and —C—C— bonds were the primary crosslinks that structurally bound kerogen monomers together. The trends observed in the simulated thermochemical transformation of kerogen and the kerogen crosslinking pathways are consistent with theoretical and experimental studies reported in the scientific literature. Reactive force field molecular dynamics simulation provides a potentially valuable approach to the development of realistic three-dimensional models for kerogen molecular networks. However, the conversion of kerogen molecules into a 3D cross-linked molecular network was low (the density of the intermolecular crosslinks was low).


2

Page 3 of 47

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Energy & Fuels

I. Introduction Kerogen is the part of the solid organic matter in buried sediment and sedimentary rocks, that is insoluble in typical organic solvents1, and it is formed by the concerted decomposition and reassembly of the particulate and dissolved organic matter deposited together with inorganic materials in sedimentary systems2. The buried organic matter undergoes a continuous chemical transformation in the subsurface, but under progressive burial conditions, where the rates of organic matter loss and change in atomic composition decrease with increasing burial depth because of consumption of electron acceptors, increased recalcitrance of the remaining organic matter and because conditions become less favorable for microorganisms. Consequently, most of the loss of organic matter and changes in its atomic composition occur at depths that are typically less than a few meters, until thermal decomposition begins on geological time scales at much greater burial depths and times. However, the discovery of microbes at depths of several kilometers3-4, suggests that enzyme catalyzed organic diagenesis may proceed at much higher depths under some conditions, albeit at rates that continue to decrease with increasing depth. At all depths, part of the organic matter is insoluble, and there is no well-defined burial depth or burial time at which the insoluble organic matter becomes kerogen. However, kerogen can be considered to be the insoluble organic matter that remains when early organic diagenesis (low temperature diagenesis at depths up to a few hundred meters) is essentially complete. At much greater depths and higher temperatures, hydrocarbons produced from kerogen may disproportionate into hydrogen rich low molecular mass hydrocarbons and an insoluble high carbon solid (pyrobitumen), which may be difficult to distinguish from high maturity kerogen (kerogen that has been transformed into a low H/C and low O/C ratio form accompanied by the production of hydrocarbon fluids), and there may be a continuum of high carbon insoluble solids


3

Energy & Fuels

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 4 of 47

ranging from high maturity kerogen (formed directly from kerogen of lower maturity) to pyrobitumen (formed only from hydrocarbon fluids produced from kerogen). It is well accepted that when kerogen and other organic matter is buried to depths at which the temperature exceed 50 C, thermally driven chemical transformation (catagenesis) commences on geological time scales. At first only low activation energy reactions such as breaking of the weakest chemical bonds and reaction of the free radicals produced by bond breaking occur. These reactions lead to the formation of oil, CO2, and H2O, a decrease in the H/C and O/C ratios in the remaining insoluble organic matter (kerogen) and an increase in the kerogen crosslink density. As progressive burial continues to increase the temperature, a wider range of chemical transformations occur, and the average molecular mass of the hydrocarbon reaction products decreases until only methane, with very small amounts of ethane and other, higher molecular mass, hydrocarbons are produced at temperatures above 150 C. Most of the world’s oil and gas resources have been produced by the thermal maturation of kerogen (the soluble bitumen associated with kerogen may also play a significant role in some cases) and the expulsion of hydrocarbons from maturing kerogen. Further, the retention of hydrocarbons in kerogen is important in primary migration, and it is one of the phenomena that limits shale oil and shale gas recovery factors. If kerogens had well determined chemical structures, molecular modeling could play an important role in developing a better understanding of and information about kerogen-fluid molecule systems, processes such as the transport of fluid molecules in nanoporous kerogen matrices and phenomena such as the retention of fluid molecules, including hydrocarbons, within kerogen and their expulsion from it. This knowledge and understanding could provide new insights and a theoretical foundation for the development of better technologies for increased


4

Page 5 of 47

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Energy & Fuels

recovery of hydrocarbons from organic-rich tight rocks (a.k.a. shales), CO2 sequestration in depleted shale and other geotechnical applications. However, kerogen is formed primarily from detrital organic matter particles with a variety of origins and chemical natures, and this complex organic material is not completely homogenized by the depolymerization, polymerization and other processes that occur during organic diagenesis. In some depositional environments, the detrital organic matter may be dominated by a particular source such as the remains of phytoplankton, but part of the organic matter may have a quite different source such as dissolved or particulate organic matter of terrestrial origin. As a result, kerogen particles separated by 1 m or less often have very different morphologies, and they can sometimes be distinguished as belonging to different maceral classes. Further, individual kerogen particles often appear to be heterogeneous5-7, as shown in Figure 1, but because of the possible presence of unresolved fluid-filled pores, inorganic mineral particles, bitumen particles and imaging artifacts, it is difficult, if not impossible, to determine how the chemical nature of the kerogen varies within the particle. Molecular models for kerogens of various types and maturities have been developed8-12 based on information gleaned from elemental analysis, spectroscopy and pyrolysis experiments13-17. These molecular models are typically constructed of carbon (C), hydrogen (H), oxygen (O), sulfur (S), and nitrogen (N) atoms that are connected together, in a manner that is consistent with the available information, by covalent bonds, hydrogen bonds and van der Walls forces8, 18-20 with a total of up to few thousand atoms (molecular masses on the order of 10,000 daltons). For comparison, typical commercial thermoplastics (melt-processable polymers) have number average molecular masses of 10,000 to 100,000 daltons, and they are soluble in some nonreactive solvents, particularly in the lower end of the molecular mass range.


5

Energy & Fuels

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 47

The construction of molecular models for various types of kerogen and levels of maturity is an important step forward. However, it is challenging because of the difficulty of isolating kerogen from the host mineral matrix — without significant alteration of the underlying kerogen structure21 and the because of the “volume” averaged nature of chemical analysis methods. Yet, substantial effort has been devoted to the construction of two-dimensional (2D)8, 11 and threedimensional (3D)9-10, 12, 22 representative molecular models for kerogen (often referred to as monomers)10, 12. Although none of these molecular models can completely represent the chemical complexity and heterogeneity of kerogen, they do capture some important chemical characteristics of kerogen, and molecular simulations based on them may provide valuable insights. To some extent, the intramolecular chemical heterogeneity of large macromolecular models may capture some of the larger scale heterogeneities within and between kerogen particles. Molecular simulations can be conducted with an ensemble of kerogen structural models (monomers) with periodic boundary conditions to enable simulations of the bulk properties of kerogen and the transport of fluids through compact and nanoporous kerogen matrices. However, this does not enable the scale of the chemical heterogeneity to be extended beyond the scale of the individual kerogen molecular model (10 nm) unless a variety of kerogen structural models are used. Very often, when kerogen monomers were packed into a periodic cube and used to construct comparatively larger kerogen macromolecular models, the crosslinking among representative component structures was either not accounted for or not explicitly emphasized10, 22-25

. Although these representative monomers are large and they do include intramolecular

crosslinking (loops), they do not capture the complexity and variability of natural kerogen, and


6

Page 7 of 47

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Energy & Fuels

they are not models of extended three-dimensional polymeric networks like those in molecular gels, rubbers, and kerogen, which are formed via chemical crosslinking. A very large body of research on polymers indicates that crosslinking provides structural stability to polymers and influences various polymer properties including solubility, swelling, susceptibility to biodegradation, mechanical properties and the transport of small molecules within the polymer26-28. Likewise, because of its insolubility in a wide range of organic solvents with different solubility parameters, kerogen is widely believed to be a crosslinked geopolymer29, and the crosslinking can be expected to play a role in phenomena that are relevant to primary migration and the recovery of oil and gas from organic-rich rocks such as adsorption/desorption of hydrocarbons and small molecules such as CO2, H2S and H2O in the kerogen matrix30-32, and the adsorption-induced swelling behavior of kerogen33-35. In this article, we describe and discuss the application of reactive force field molecular dynamics to simulate an ensemble of kerogen molecules in order to understand the thermochemical transformation of kerogen (i.e., breaking and formation of chemical bonds), to obtain molecular insights into possible crosslinking pathways, and to develop preliminary quasi-infinite macromolecular networks for kerogens of various thermal maturities, which might be relevant to the development of a better fundamental understanding of the generation, retention and transport of hydrocarbon fluids in kerogen. II. Methodology Molecular dynamics using the general bond-order-dependent ReaxFF reactive force field36 was used to simulate bond breaking and both intra- and inter-molecular bond formation during the rapid, high temperature heating of kerogen. The ReaxFF force field, which has been widely used to study the pyrolysis and combustion of hydrocarbon systems17, 36-39, enables the simulation of


7

Energy & Fuels

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 47

chemical reactions without the computational burden of quantum mechanical molecular dynamics, and systems with thousands of atoms as well as their chemical reactions can be investigated. The total energy in the ReaxFF reactive force field can be described by 𝐸𝑡𝑜𝑡𝑎𝑙 = 𝐸𝑏 + 𝐸𝑜𝑣𝑒𝑟 + 𝐸𝑢𝑛𝑑𝑒𝑟 + 𝐸𝑙𝑝 + 𝐸𝑣𝑎𝑙 + 𝐸𝑡𝑜𝑟 + 𝐸𝑣𝑑𝑤 + 𝐸𝐶𝑜𝑢𝑙

(1)

where Eb, Eover, Eunder, Elp, Eval, Etor, Evdw, and ECoul are the bond, overcoordination, undercordination, lone-pair, valence, torsion, van der Waals and Coulomb energies respectively40-41. The terms in the ReaxFF force field, eq 1, can be divided into those that depend on bond orders, and those that do not, and the bond orders are calculated from the corresponding interatomic distances. Additional details concerning the different energy contributions and the ReaxFF reactive force field can be found in the scientific literture36, 42. Figure 2 illustrates the ReaxFF MD simulation workflow used in the work reported in this paper to simulate the thermo-chemical transformation of kerogen and identify potential crosslinking pathways. At the start of each simulation, 10 randomly oriented monomers of immature type I kerogen (kerogen 1-I) with atomic composition C251H385O13N7S3, low maturity type II kerogen (kerogen 2-L) with atomic composition C234H263O14N5S2 or low maturity type III kerogen (kerogen 3-L) with atomic composition C233H204O27N4 (as proposed and classified by Ungerer et al.12 and shown in Figure 3), with a total of 6590, 5180, and 4680 atoms, were placed in a 7nm x 7nm x 7nm cubic simulation box, with periodic boundary conditions. The molecular structures of the monomers proposed by Ungerer et al.12 were based primarily on matching elemental analysis and spectroscopic data15. For each kerogen model, a multistage simulation was conducted (Table I summarizes the simulation schedule) in order to investigate the mechanisms of kerogen maturation. During the first stage, the ensemble was equilibrated at a temperature of 300 K and a pressure of 100


8

Page 9 of 47

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Energy & Fuels

atmospheres for 25 ps to obtain initial molecular conformations and packings with kerogen densities similar to the densities measured under laboratory conditions (temperatures of about 300 K and pressures of 1 atmosphere) before the kerogen structures were used in the subsequent simulation stages. For kerogen 1-I, the density obtained after 300 K and 100 atmospheric pressure simulation stage was 1.03 g/cm3, which was in the close agreement with the measured density of kerogen 1-I densities (0.95 g/cm3, 43and 0.968 g/cm3,44) but for kerogens 2-L and 3-L the densities after stage 1 (0.86 g/cm3 and 0.848 g/cm3) were somewhat lower than their measured densities (1.18 g/cm3 and 1.25 g/cm3, 45). Subsequently, the temperature of the ensemble was heated linearly to 1400 K over a period of 200 ps to investigate bond breaking processes on the molecular dynamics time scale. The system was then cooled in stages to simulate bond forming processes and investigate whether ReaxFF molecular dynamics simulations could be used to generate larger macromolecules or a quasi-infinite periodic molecular network (through the periodic boundary conditions). The maximum temperature of 1400 K was selected based on the previous work of Liu et al.17, who used reactive force field molecular dynamics to simulate the pyrolysis of the Green River Shale kerogen. Liu et al.17 observed Green River Shale kerogen pyrolysis at temperatures in the 1200–1600 K range (though the maximum pyrolysis temperature used in their study was 2800 K), thus indicating that the kerogen pyrolysis presumably begins at a temperature of about 1200 K on a molecular dynamics time scale. It is worthwhile to discuss further the relevance of high-temperature ReaxFF MD simulations on very short time scales (typically on the order of a few nanoseconds) to the slow thermochemical transformation of kerogen during catagenesis that occurs on geological time scales (millions to


9

Energy & Fuels

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 47

tens of millions of years). The rate constants of chemical reactions can usually be represented quite well by the modified Arrhenius equation

𝑘(𝑇) = 𝐴𝑇 𝑛 𝑒

−𝐸𝑎 ) 𝑘𝐵 𝑇

(

(2)

or the transition state theory expression

𝑘(𝑇) = (𝑇)(𝑘𝐵 𝑇/ℎ)𝑒

−𝐺‡ (𝑇) ) 𝑘𝐵 𝑇

(

,

(3)

where k(T) is the rate constant at absolute temperature T, A is the frequency factor (characteristic of a vibrational frequency), Ea is the activation energy, kB is the Boltzmann constant, h is the Planck constant and G‡ is the quasithermodynamic Gibbs activation free energy. In eq 2, the exponent, n, is typically quite small (usually within the range -1  n  1), and in eq 3, the transmission coefficient (𝑇) is typically weakly temperature dependent, with values on the order of unity. For simple covalent bond breaking processes, Ea  G‡  Eb, where Eb is the bond energy, and because Eb >> kBT except at extremely high temperatures, it is sufficient to replace eqs 1 and 2 by 𝑘(𝑇) = 𝐴′𝑒

−𝐸′𝑎 ) 𝑘𝐵 𝑇

(

or 𝑘(𝑇) = 𝐴′𝑒

−𝐸′′ 𝑎) 𝑅𝑇

(

,

(4)

over a quite large range of the temperatures for the purposes of semiquantitative analysis, where 𝐴′ is the effective frequency factor, 𝐸𝑎′ is the effective bond-breaking activation energy, 𝐸𝑎′′ is the molar bond energy and R is the universal gas constant. Assuming a constant progressive burial rate of 10-3 meters per year, a constant geothermal gradient of 25 K per kilometer and an oil window that extends from 340–425 K, the kerogen in a sedimentary basin would be heated at a rate of 25 K per My, and it would require about 3 My to pass through the oil window. This is 1023 times the physical time scale of typical reactive force


10

Page 11 of 47

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Energy & Fuels

field simulations and almost 1030 times the duration of each time step (0.1–0.5 fs). Based on an analysis of oil production in the subsurface by Quigley and MacKenzie46, most of the oil would have been produced at a temperature of  415 K over a period of  1 My. If it is assumed that the kinetics of kerogen maturation can be represented by eq 4, the activation energy, 𝐸𝑎𝑚 required to achieve the same extent of reaction in a molecular dynamics simulation at a temperate of 𝑇𝑀𝐷 over a simulation time of 𝑡𝑀𝐷 and in the subsurface at a temperate of 𝑇𝑠𝑠 over a geological time of 𝑡𝑠𝑠 is given by

𝐸𝑎𝑚

=𝑅

𝑘 ln 𝑀𝐷

𝑘𝑠𝑠 1 1 − 𝑇𝑠𝑠 𝑇𝑀𝐷

=𝑅

𝑡 ln 𝑠𝑠

𝑡𝑀𝐷 1 1 − 𝑇𝑠𝑠 𝑇𝑀𝐷

(5)

where 𝑡𝑠𝑠 is the time over which most oil production occurs in the subsurface and 𝑡𝑀𝐷 is the time over which most oil production occurs in the molecular dynamics simulation. Based on the temperature-time schedule used in the reactive force field molecular dynamics simulations, it is reasonable to expect that most of the simulated reaction chemistry occurred at the higher temperatures used in the simulations (1300 K). For 𝑡𝑠𝑠 = 1 My (~ 3.154  1013 s), 𝑡𝑀𝐷 = 10-9 s, 𝑇𝑠𝑠 = 415 K and 𝑇𝑀𝐷 = 1300 K, eq 5 predicts that 𝐸𝑎𝑚  62.8 kcal mol-1. Changing 𝑡𝑠𝑠 /𝑡𝑀𝐷 by a factor of 10 changes 𝐸𝑎𝑚 by about 2.75 kcal mol-1, changing 𝑇𝑠𝑠 by 25 K changes 𝐸𝑎𝑚 by about 5.5 kcal mol-1 and changing 𝑇𝑀𝐷 by 100 K changes 𝐸𝑎𝑚 by about 2.25 kcal mol-1. In practice, sediment accumulation rates and the concomitant burial rates can vary greatly from one depositional environment to another47. If the burial rate and heating rate are much lower, say 10-3 meters per year and 0.25 K per My, then 𝑡𝑠𝑠 /𝑡𝑀𝐷 will be much larger. However, if the heating rate is lower, there is a compensating reduction in the temperature at which most oil is produced46, and the activation energy needed to match the extent of reaction during a 1 ns molecular dynamics simulation at a temperature of 1300 K to geological times and temperatures


11

Energy & Fuels

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 47

does not change significantly. Of course, the time dependence of the burial depth and time temperature history of kerogen in sedimentary basins may be quite complex. However, a constant heating rate is often used in laboratory experiments and numerical models, and it facilitates theoretical analysis. Quigley and MacKenzie46 assumed constant heating rates of 0.1 – 100 K per My in their analysis of the temperatures of oil and gas formation in the subsurface. Although a single activation energy, typically 45 – 60 kcal mole-1, 48-49 (188-251 kJ mol-1 or 3.13 – 4.17 × 10-19 J per bond), has often been used to represent the overall kinetics of kerogen pyrolysis, the plethora of reactions that occur during kerogen maturation and pyrolysis have a broad distribution of activation energies and the kinetics of kerogen pyrolysis is also frequently represented by a continuous or discrete distribution of parallel reactions or by frequency factors and activation energies that vary with the extent of pyrolysis. Although the distribution of activation energies obtained by analyzing the data obtained from pyrolysis experiments is typically quite broad when the kinetics is represented by a set of parallel reactions, it is typically dominated by activation energies that also lie in the 45 – 60 kcal mol-1 range48-51. The activation energy of 𝐸𝑎𝑚 = 62.8 kcal mol-1 needed to achieve the same extent of bond breaking after 1 My at 415 K and after 1 ns at 1300 K, based on the simple single activation energy Arrhenius kinetics model (eq 4) is just outside this range. Neither the distribution of activation energies nor the single activation energy, if a one reaction model is used, should be taken too seriously.As mentoned previuously, if the burial rate is ten times slower or ten times faster, the activation energy needed to achieve this time-temperature superposition changes by only a little more than 2.5 kcal mole-1, and if the characteristic oil window temperature is decreased or increased by 50 K, the activation energy changes by about 10.0 kcal mol-1. This indicates that reactive force field molecular molecular dynamics simulations at high temperatures and on very short time


12

Page 13 of 47

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Energy & Fuels

scales can indeed provide insights concerning processes that occur in sedimentary basins on geological time scales in much the same way that pyrolysis (T  700 K, t = 104-105 s) and flash pyrolysis (T  900 K, t = 0.1-1.0 s) have been used to obtain information relevant to kerogen maturation.52-55 However, the relative rates of reactions with different activation energies is highly temperature dependent if the difference in activation energies is large, and high activation energy reactions that are important at temperatures much greater than those at which natural kerogen maturation occurs, at depths less than 10 km, may be unimportant during natural kerogen maturation. The LAMMPS molecular dynamics code56 was used to perform all of the simulations with isothermal-isobaric (NPT) ensembles and periodic boundary conditions in all directions. Time integration of the Nose-Hoover style non-Hamiltonian equations of motion was used to generate atom positions and velocities sampled from the NPT ensembles56. A small time step of 0.25 fs and a simulation time of 1.025 ns were used to capture the bond breaking and formation processes, and thus the crosslinking reactions of kerogen. The atomic positions were recorded every 500-time steps. In addition, the bonds that connect different monomers were identified every 10,000-time steps or every 2.5 ps. Further, the total number of bonds cleaved and a total number of bonds formed was recorded every 50 ps. A similar simulated (P, T) schedule was used for all kerogen types. III. Results A. Identification of molecular moieties participating in bond breaking and bond formation reactions A comparison of the chemical bonds in the initial kerogen ensemble with the chemical bonds in the intermediate and the final kerogen structures, as well as the low molecular mass products


13

Energy & Fuels

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 47

such as small molecules and free radicals generated during the reactive molecular dynamics simulations of kerogen, revealed the possible reaction pathways; productions of radicals and small molecules; and the formations of crosslinking covalent bonds between kerogen monomers. The chemical bonds (in the ReaxFF force field formulation) were determined at every timestep as follows. The atoms in a kerogen ensemble experience various interatomic forces (as described in eq 1 during every simulation timestep, which may influence the interatomic distance and hence the interatomic distance dependent bond order. The bond order is higher at short interatomic distances (thus allowing the formation of chemical bonds between atoms) and decreases with increasing interatomic distance before becoming zero at relatively larger interatomic distances (thus breaking the chemical bond between atoms). Bond formation and bond breaking persist throughout the simulation57. Figure 4 shows various molecular moieties that underwent bond cleavage reactions during the simulations to form highly reactive sites, which may eventually participate in the formation of crosslinks between kerogen monomers. The bond breaking processes shown in Figure 4 are only some of the bond breaking processes that form reactive sites that may participate in subsequent crosslinking reactions during the high-temperature stages of the simulations of rapid kerogen maturation. In addition to participating in crosslinking reactions, reactive sites formed by bond breaking may result in the generation of small molecules, isolated small molecular clusters or the formation of covalent bonds that do not crosslink kerogen molecules. In kerogen 1-I, bond cleavage occurs primarily in long aliphatic chains containing —C—O— and —C—C— covalent single bonds, which typically have bond breaking energies of ~85 kcal mol-1 58-59 unless one or both of the reaction products is stabilized by electron delocalization or the bond energy is reduced by strong steric repulsion, which does not occur in the 1-I, 2-L or 3-L kerogen models


14

Page 15 of 47

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Energy & Fuels

used in this work, and is very unlikely to play a significant role in natural kerogens. The substantially weaker —C—S— bonds with bond strengths of ~60 kcal mol-1 59 are rapidly broken, but there are much fewer —C—S— than —C—O— and —C—C— bonds present in kerogen 1-I and 2-L, and kerogen 3-L contains no sulfur. Cleavage of bonds between aliphatic and aromatic carbons was also observed (Figure 4a). In the case of kerogen 2-L and kerogen 3-L, additional molecular moieties participated in bond cleaving and reactive site generation. Despite the high dissociation energies of the bonds between hydrogen atoms and aromatic moieties, including furanic groups, the scission of these —C—H— bonds also occurred in the simulation of kerogen 2-L (Figure 4b). The frequency factors, A in eqs 2– 4, for the scission of —C—H— bonds is higher than the frequency factors for the scission of —C—C— bonds and there are more bonds between aromatic carbons and hydrogen than single —C—C— and —C—O— bonds. Nevertheless the ratio between aromatic carbons–hydrogen and single —C—C— and —C—O— bond breaking was surprising. Also, thiophenic sulfur is a reactive site in kerogen 2-L. In the case of kerogen 3-L, the scission of carbon–carbon and carbon–oxygen single bonds, which have similar bond dissociation energies was dominant. In particular, the scission of —C—C— bonds in short aliphatic chains linking aromatic moeities played an important role in the generation of reactive radicals in kerogen 3-L. During the last century, a very large body of information about the mechanisms of organic reactions has been acquired, and concepts based on this information and quantum mechanics have provided an empirical and theoretical foundation for understanding chemical transformations that strongly influences the interpretation of experiments on the chemistry of complex organic materials such as kerogen. The majority of the bonds cleaving and bond forming (e.g., crosslinking) reactions are consistent with previous conclusions reported in the


15

Energy & Fuels

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 47

scientific literature. Behar et al.8 constructed kerogen models based on elemental analysis, 13C NMR spectroscopy, thermogravimetry, functional analysis and pyrolysis experiments, and observed that aliphatic —C—C— bonds in immature type I Green River shale kerogen, nitrogen containing functional groups (e.g., amide nitrogen) in the oil window type II kerogen, and oxygen containing functional groups (e.g., ethers) in matured type III kerogen are the dominant bond cleaving sites. On the other hand, Kerichko and Gagarin60 found that the ethers and aliphatic carbon connect aromatic carbon rings, thus confirming the role of ethers and aliphatic carbons as potential crosslinking bonds. Similarly, Bushnev et al.18 performed NMR analysis and showed that sulfide bridges crosslink aliphatic carbons, thus also confirming the role of sulfide bridges in (crosslinking) bond formation. Numerous other studies61-62 led to the similar conclusions. Further, it is likely that numerous bond cleaving reactions that are difficult to identify and characterize experimentally occur during thermal maturation of kerogen. The analysis of newly formed bonds in the various kerogens indicates that the number of —C—S— bonds formed (post initial bond cleavage at higher temperatures) was significantly larger than the number of —C—O— and —C—C— bonds formed in kerogen 1-I and kerogen 2-L, whereas more —C—C— bonds were formed in kerogen 3-L. This suggests that the weaker —C—S— bond in kerogen 1-I and kerogen 2-L configurations underwent substantial bond breaking reactions resulting in the generation of reactive sites and the subsequent formation of new —C—S— bonds. The rates of thermally driven bond cleavage depend primarily on the bond energy. Consequently, it is not surprising that we were able to identify only single bond cleavage processes, and that weak —C—S— bonds and both —C—C— and —C—O— bonds that result in the formation of one or two radicals stabilized by electron delocalization were preferentially cleaved (Figure 4).


16

Page 17 of 47

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Energy & Fuels

B. Thermochemical transformation of molecular moieties at elevated temperature and pressure Figure 5a shows the variation in the numbers of ruptured and newly formed bonds and the change in the total number of bonds during the simulated maturation of kerogen 1-I (with the (P, T) schedule specified in Table 1) from the beginning of the simulations at intervals of 50 ps. As expected, the heating temperature and heating/quenching time play a significant role in the bond breaking and formation process63. Bond cleaving and bond forming processes occurred during the heating of kerogen from 300 K to 1400 K (from 0 to 0.2 ns) — particularly at the higher temperatures (> 900 K), when the amount of vibrational energy associated with each chemical bond fluctuates until sufficient energy to break the chemical bond is accumulated. While the temperature was reduced from 1400 K to 1200 K (from 0.2 to 0.6 ns) a noticeable increase in bond cleavage and net bond formation occurred (as shown in Figure 5a). This seems surprising, but it might be a consequence of the production of reactive free radicals during the very early stages of the simulation that subsequently undergo rapid bond cleavage processes such as  scission (radical depolymerization of alkyl chains) and radical polymerization resulting in multiple bond cleavage and bond formation processes. Additional cooling below 1000 K (beyond 0.8 ns) reduced the bond cleaving/net bond formation rates, and at 300 K there was no reaction on the molecular dynamics time scale. Similar observations were made for kerogen 2-L and kerogen 3-L (Figures 5b and 5c). Apparently, and as expected, the mobility of free radicals and small molecules, and the reaction rates are higher at elevated temperature and decrease with decreasing temperature. The polyaromatic rings, ketone (—C=O), and the majority of the carboxylic acid (—COOH) functional groups remained stable at the elevated temperature, except that a few carboxylic sites


17

Energy & Fuels

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 47

underwent bond cleaving reactions to generate formic acid (HCOOH). Additional chemical reactions that were observed include the dissociation of —OH from —COOH functional groups and its subsequent reaction with hydrogen atoms resulting in the formation of H2O molecules. This reaction is not fully understood at this stage and warrants further investigation to assess whether the ReaxFF forces field is able to correctly capture carboxylic acid reactions in the complex chemical environment of kerogen and to identify whether there are any additional factors that influence the behavior of carboxylic functional groups. Ether (—C—O—), aliphatic carbon (—C—C—), and weak sulfide (—C—S—) bonds underwent bond cleaving reactions resulting in the formation of labile small molecular fragments, free radicals and highly reactive sites on the parent kerogen structure. The —C—S— bonds (from aliphatic sulfur) are weaker than —C—C— and —C—O— bonds and undergo multiple bond cleaving/bond formation reactions before forming a stable crosslinking —C—S—C— bond. In addition, hydrogen atoms covalently bonded to aliphatic and aromatic carbon dissociates to form reactive sites on the carbon atoms. The newly formed free radicals, small molecules, reactive sites on the parent monomer structure are highly unstable at elevated temperatures, they do not remain unreacted for very long, and they form more stable secondary reaction products64. The free radicals, undergo further bond cleavage (particularly via  scission), preferentially, react with available local reactive sites or migrate to more distant reactive sites, which results in either extension of existing kerogen molecules or the formation of non-reactive stable molecules (e.g., H2O and C2H4). The migration of free radicals and small molecules depends on the local chemical environment (physical arrangement and dynamics of various molecular moieties). The radicals and other highly reactive products that are not able to find sites to react with remain in


18

Page 19 of 47

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Energy & Fuels

the kerogen. Similarly, it has been demonstrated experimentally65-66 that demineralized kerogens contain uncoupled electrons at concentrations of 1017 – 1019 cm-3. The newly developed reaction sites (e.g., carbon radicals) can also satisfy their reactivity through covalent bond formation with hydrogen atoms present in the kerogen or with olefinic moieties formed by  scission. We detected ethylene in all of the simulations, but hydrogen atoms were not detected – perhaps because they have very short lifetimes. The secondary products generated during the simulations of kerogen 1-I, kerogen 2-L and kerogen 3-L were also characterized and compared with experimental data available in the scientific literature. Because the degree of pyrolysis was small, it is expected that the secondary pyrolysis products will show some, but not all, of the characteristics of typical pyrolysis/flash pyrolysis experiments in which kerogen is typically completely pyrolyzed. It is also possible that short lived intermediate fragments generated in the simulations were not observed experimentally because they decomposed or reacted with other fragments before they could be detected. Table 2 summarizes the low molecular mass products of the kerogen thermal maturation simulations. As expected, because of its higher oxygen content, kerogen 3-L generates the most H2O. Further, kerogen 2-L and kerogen 3-L generated C3H6 whereas kerogen 1-I produced C4H8. In addition, for kerogen 1-I, water (H2O), ethylene (C2H4), methanol (CH3OH) and the unstable fragment OCHOH were observed as the major reaction products. H2O, C2H4, methanol (CH3OH), and formaldehyde (CH2O) molecules were observed for kerogen 2-L whereas H2O, C2H4, vinyl alcohol (CH2CHOH), CH2O and CH2CHCHOH were produced from kerogen 3-L. Finally, numerous highly reactive products that are not observed experimentally were generated from all of the kerogens. This is not surprising at the high temperatures used in the simulations


19

Energy & Fuels

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 47

since the rates of bond breaking processes increase more rapidly with increasing temperature than the rates of bond forming processes because of their higher activation energies. Kerogen 1-I produces H2O as well as alkanes during the early stages of maturation, and this is consistent with the results of pyrolysis experiments. Further, the reactive force field molecular dynamics simulations of Morwell coal (immature lignite enriched in aromatic moieties) performed by Salmon E., et al.39 shows the generation of H2O, C2H4, and CH2O molecules, which are consistent with the reaction products generated from kerogen 3-L in the present investigation. C. Characteristics of crosslinked kerogen structure Crosslinks are covalent bonds or atoms and bonds that connect different monomers to form a larger 3D polymeric structure as shown in Figure 6. The low molecular mass products that do not contribute to the remaining, more mature, kerogen were disregarded during the characterization of the largest crosslinked kerogen structure using cluster analysis. Essentially, two types of crosslinks were observed — the main crosslinks that connect different partially degraded kerogen monomer to a form larger kerogen structure and the local crosslinks that connect small molecular fragments to the main carbon backbone. The local crosslinks did not significantly contribute to the largest continuous kerogen structure formed in the simulations. The largest kerogen structures formed from kerogen 1-I, kerogen 2-L, and kerogen 3-L, are shown in Figure 6. They form a quasi-infinite periodic molecular networks via a connection in all three-directions across the periodic boundary conditions used in the molecular dynamics simulations. The main crosslinks vary based on the type of kerogen. Sulfide (—C—S—) and carbon-carbon (—C—C—) bonds were the primary crosslinks formed between kerogen 1-I molecules, and between kerogen 2-L molecules (Figures 6a and 6b), whereas —C—C— bonds were the main crosslinks formed between kerogen 3-L molecules (Figure 6c). Interestingly, no —C—O—


20

Page 21 of 47

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Energy & Fuels

crosslinks were formed by kerogen 1-I or by kerogen 2-L molecules, whereas —C—S— crosslinks were not formed between kerogen 3-L molecules because the original kerogen monomer did not contain any sulfur atoms. The newly formed main crosslinks appeared to be randomly distributed in the final crosslinked kerogen structure. Experimental results from the literature show similar trends in kerogen crosslinking. Bushnev et al.18 performed 13C NMR analysis and observed that the sulfide bridges connecting aliphatic carbons help to form a 3D crosslinked kerogen structure. Richnow et al.62 also confirmed the role of sulfur and oxygen atoms containing functional groups (e.g., ethers) as cross-linkages connecting various kerogen building blocks. Moreover, the majority of relevant studies62, 67-68 show that —C—C—, —C—O—, and —C—S— bonds are the primary crosslinks that connect various functional groups to form 3D crosslinked kerogen structures. The effect of thermal maturity on kerogen crosslinking is consistent with experimental observations, which indicate that low maturity type I kerogen (kerogen 1-L) shows the least crosslinking among all kerogen types69. Further, kerogen 2-L displayed the highest crosslink density and resulted in the largest residual kerogen structure (shown in Figure 6b). The H/C and O/C ratios of the residual crosslinked kerogen 2-L did not change substantially during simulated maturation/pyrolysis (Table 3). In other words, the position on the van Krevelen diagram70 did not change substantially. Consequently, the degree of crosslinking that occurs in the simulations is expected to be smaller than the degree of crosslinking in natural kerogen of the same type and of the same maturity that has undergone crosslinking during the entire catagenesis process. In addition, we also determined the carbon-carbon radial distribution function (a.k.a. pair correlation function) g(r) of the resultant crosslinked kerogen 2-L structure up to a radial distance of 10 Å to verify consistency with carbon-carbon radial distribution functions that have


21

Energy & Fuels

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 47

been obtained from other molecular simulations of kerogens that have been reported in the scientific literature12. From Figure 7, it can be seen that the first peak occurs at a radial distance of 1.41 Å, which corresponds to the lengths of aromatic carbon–carbon bonds (the shortest C—C bond lengths since there are no carbon-carbon double bonds or triple bonds). The subsequent peaks at the radial distance around 1.49 Å and 1.55 Å correspond to bonds with relatively larger bond length including aliphatic carbon–aliphatic carbon bonds as well as aliphatic carbon– aromatic carbon single bonds. The peaks beyond a radial distance of 2 Å correspond to pairs of carbon atoms that are linked through two or more bonds. The locations and intensities of the peaks observed in g(r) are in close agreement with those obtained from other molecular simulation that investigated kerogen 2-L12. The shape of g(r) (Figure 7) for r < 2.0 Å suggest that it might be possible to separate this part of the radial distribution function into three, or possible more, overlapping peaks. However, because the peak shapes are not known, the peak areas cannot be reliably determined. However, the area under the peak at a radial distance of 1.41 Å appears larger than the peak areas at 1.49 Å and 1.55 Å thus indicating that the residual crosslinked kerogen 2-L structure contains a higher percentage of aromatic-aromatic carboncarbon bonds than other carbon-carbon bonds and a higher percentage of aromatic carbon atoms than aliphatic carbon atoms. IV. Discussion Kerogens, which are insoluble in all non-reactive solvents, are widely believed to be crosslinked polymers. Experimental estimation of the degree of crosslinking is usually based on the Flory Rehner theory71, which takes into account the entropy of mixing, the configurational entropy of the crosslinked polymers and the enthalpy of mixing (based on regular solution theory). Discrepancies between the Flory Rehner theory and experiments72 are often attributed to


22

Page 23 of 47

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Energy & Fuels

inadequacy of the assumption that the polymer chains that extend from one crosslink to another undergo affine deformation or the assumption that the free energy of mixing and the free energy of deformation of the polymer network can be treated as additive contributions to the total free energy. However, for kerogens, other assumptions of the Flory Rehner conceptual model may be more important. These include: the assumption that the polymer network, including both the crosslinks and the polymer chains between them, is chemically homogeneous and can be characterized by a single solubility parameter; the assumption that the contour lengths of the polymer chains that connect one crosslink to another are the same; and the assumption that the change in free energy of the polymer network is purely entropic. Crosslink densities estimated from swelling and stress-strain measurements using Flory Rehner theory or modifications of Flory Rehner theory73 are based on similar assumptions, and they result in effective crosslink densities. A variety of processes may contribute to the effective crosslink density74. These include covalent chemical bond formation, the formation of ionic bonds by the interaction of functional groups with multivalent cations, the formation of hydrogen bonds, the formation of nanoscale crystallites, and entanglement. In the case of kerogens, it is unlikely that nanoscale crystallinity contributes to the effective crosslink density, but in high maturity kerogens, the stacking of polyaromatic sheets may act much like the formation of crystallites. It is also possible that the heterogeneity of kerogen contributes to its insolubility. Different parts of the kerogen molecules can be expected to have different solubilities, characterized by their solubility parameters within the framework of regular solution theory, and there may be no solvent with a solubility parameter that matches the solubilities of all of the parts of the kerogen well enough to dissolve it (there may be a net attractive kerogen-kerogen interaction energy that is larger enough to prevent dissolution by all non-reactive solvents).


23

Energy & Fuels

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 24 of 47

In practice, little is known about the nature of crosslinking in kerogens because of the complexity and variability of kerogens and the difficulty of measuring the effective crosslinking density. For synthetic crosslinked polymers, which undergo essentially irreversible chemical crosslinking, spectroscopic methods75 may be used to determine the density of the chemical moieties that are characteristic of chemical crosslinking, but in most cases, it is not possible to distinguish between intermolecular and intramolecular crosslinking. In the case of kerogen, which is formed and transformed by both bond breaking and bond formation processes, the formation of new chemical moieties that are characteristic of crosslinking does not provide a reliable estimate of the extent of crosslinking because bond breaking processes may disconnect them from the crosslinked chemical network. Molecular simulations, like those reported here, may provide important insights into the nature of chemical crosslinking in kerogen that could support a more reliable interpretation of spectroscopic data in terms of chemical crosslink density. Most experimental investigations of the thermal alteration of kerogen (e.g., pyrolysis, T  700 K, t = 104-105 s or flash pyrolysis, T  900 K, t = 0.1-1.0 s) are conducted at temperatures that are much closer to the temperatures at which natural kerogen maturation occurs than the temperatures of up to 1,400 K used in our reactive force field simulations, and the chemistry that occurs in the experiments does differ from that which occurs under subsurface conditions. Although the difference in chemistry between natural kerogen maturation and the chemistry observed in our reactive force field simulations is larger, because of the larger difference in temperature and because we rely on a molecular model that cannot capture the variability and full complexity of natural kerogen, the simulations do provide valuable and detailed insights into the possible kerogen maturation reaction pathways on atomic/molecular length scales. This includes information about the location of bonds that are preferentially cleaved as well as bond


24

Page 25 of 47

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Energy & Fuels

connectivity and conformational information that cannot be determined experimentally. The range of extrapolation between the time scale and temperature of kerogen maturation in laboratory experiments and the time scale and temperature of kerogen maturation in ReaxFF force field molecular dynamics simulations is not very different from the range of extrapolation between the time scale and temperature of kerogen maturation in the subsurface time scale and time scale temperature of kerogen maturation in the laboratory. Based on eq 4, the ratio between reaction rate constants at temperature T1 and temperature T2 is given by 𝑘(𝑇1 ) 𝑘(𝑇2 )

=𝑒

𝐸 1 1 ( 𝑎 ( − )) 𝑅 𝑇2 𝑇1

(6)

and 1/T2 – 1/T1 is a measure of the range of extrapolation. In the case of laboratory pyrolysis experiments at a temperature of 700 K and oil formation in the subsurface at a temperature of  400 K, 1/400 – 1/700 = 1.07  10-3, and for flash pyrolysis at a temperature of 900 K, 1/400 – 1/900 = 1.39  10-3. In the case of laboratory pyrolysis experiments at a temperature of 700 K and molecular dynamics simulations at a temperature of 1300 K, 1/700 – 1/1300 = 0.66  10-3, and for flash pyrolysis at a temperature of 900 K, 1/900 – 1/1300 = 0.34  10-3. Consequently, if eq 4 is (approximately) correct, extrapolation between laboratory conditions and subsurface conditions is larger than extrapolation between molecular dynamics simulations and laboratory experiments. Consequently, it would be reasonable to expect that the differences in pyrolysis chemistry between pyrolysis experiments and reactive force field molecular dynamics simulations should be comparable with or smaller than the differences between the pyrolysis chemistry in laboratory experiments and the chemistry of kerogen maturation in the subsurface, providing the confinement conditions and chemical environments are essentially the same, and the ReaxFF force field is accurate enough.


25

Energy & Fuels

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 26 of 47

Despite the ability of reactive molecular dynamics simulations to provide qualitative information about slow chemical transformations by using temperatures that are high enough to enable them to occur on molecular dynamics time scales, it is worthwhile to describe the limitations of this approach76. First, the replication of any realistic macroscopic laboratory experiment or the natural thermal maturation of kerogen (~1016 s) is beyond the capability of any molecular dynamics simulations due to the range of time scales associated with the simulations (physical time steps of ~10-15 s and total simulation times of ~10-9 s if the ReaxFF force field is used, and the system contains a few thousand atoms). Further, even advanced reactive force fields such ReaxFF36, which allow on-the-fly bond formation and bond breaking, suffer from limitations such as unreasonable charge distributions for complex structures like kerogen, which does not have a precisely known equilibrium state, and inaccurate representtion of dispersion interactions76. Also, relatively high computation times are an impediment in the application of ReaxFF force field molecular dynamics to simulate significantly larger kerogen macromolecular structures. In principle, quantum molecular dynamics could be used, but the computation times would be even longer. A strong correlation between reported activation energies and frequency factors has been observed49 – the higher the activation energy, the higher the frequency factor. For a “typical” activation energy of  52.5 kcal mol-1, the frequency factor is, on average, a few times 1013 s-1 (on the order of 1013.5) for a variety of kerogens and organic rich rocks, and it has been observed that the combination of frequency factors and activation energies over the ranges 1010 < A < 1017 s-1 and 40 < Ea < 65 kcal mol-1 give approximately the same rate constant at a temperature of about 725 K49 (the average temperature at which most hydrocarbon is produced in typical kerogen pyrolysis experiments49). This suggests that the difference between the various


26

Page 27 of 47

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Energy & Fuels

combinations of A and Ea reported in the literature may be more a matter of how the experiments are conducted and analyzed than fundamental differences between the kerogen and rock samples. The high value of the activation energy 𝐸𝑎𝑚 (eq 5) required to achieve a similar extent of kerogen maturation in reactive force field molecular dynamics simulations on a time scale on the order of 10-9 s and on geological time scales in the subsurface relative to the activation energies obtained from most pyrolysis experiments does not imply that the effective activation energy in the molecular dynamics simulations was high. On the contrary, the high value of 𝐸𝑎𝑚 can be attributed to the surprisingly fast alteration of the kerogen in the simulations. Based on the amount of bond breaking that occurred over a period of about 500 ps at an average temperature of about 1300 K, a rate constant of about 3  108 s-1 is estimated at this temperature. If an activation energy of 52.5 kcal mol-1 and a frequency factor of 3  1013 s-1, which are consistent with pyrolysis experiments49 are assumed, then a rate constant of 4.47  104 s-1 is predicted, and if a frequency factor of 3  1013 s-1 is assumed, then an activation energy of only 29.7 kcal mol-1 is required to obtain a rate constant of 3  108 s-1 at 1300 K based on eq 4. Similarly, if an activation energy of 52.5 kcal mol-1 is assumed, then a frequency factor of 2.01  1017 is required to obtain a rate constant of 3  108 s-1 at 1300 K. While there is substantial uncertainty in the effective rate constants and the effective temperature in the molecular dynamics simulations reported here and the use of a single activation energy instead of a distribution of activation energies or conversion dependent activation energies and or frequency factors, the effective rate constants estimated from the simulations do appear to be substantially larger than expected. Other molecular dynamics simulations of the pyrolysis of kerogens using a ReaxFF force field have resulted in surprisingly fast pyrolysis rates that do not appear to be compatible with the rates obtained from pyrolysis experiments17, 77-78, and when activation energies and frequency


27

Energy & Fuels

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 28 of 47

factors were estimated for the simulated pyrolysis of Gloeocapsomorpha prisca microfossils (algal oil shale kerogen) a frequency factor of 5.6  1014 and an activation energy of 37  2 kcal mol-1 was estimated at a density of 1 g cm-3, 79. The apparent discrepancies between the rates of pyrolysis in ReaxFF force field molecular dynamics simulations and the extrapolation of experimental pyrolysis kinetics to simulation temperatures based on eq 4 might be attributed to inadequacy of this simple equation. However, the discrepancies are large and there are other possibilities. Since the model kerogen monomers contain no free radicals, production of free radicals via bond breaking plays a critically important role in the simulations, and the difference between gas phase bond energies estimated from the ReaxFF force field may differ from the bond energies obtained from more accurate quantum mechanical calculations by more than 10 kcal mol-1 in some cases80. Consequently, a few artificially weak bonds could have a large effect on the simulated pyrolysis kinetics. Other effects should be considered, such as the effect of density on molecular mobility and the effects of confinement (in the simulations no low molecular mass substances are removed, but in pyrolysis experiments volatile products are often removed and analyzed and the weight loss measured). If the density of the model kerogen is lower than the experimental density of the kerogen, the cage effect that increases the recombination of radical pairs produced by chain scission can be expected to be weaker in the simulations than in the experiments and the time required for a radical to migrate to a reactive site can be expected to be shorter. Both of these effects this could contribute to the high rate of pyrolysis observed in the simulations. In molecular simulations of kerogen pyrolysis, only organic matter is present while in whole rock pyrolysis experiments minerals that might affect kerogen pyrolysis are present, and even after chemical separation small quantities some minerals, possibly formed during chemical separation,


28

Page 29 of 47

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Energy & Fuels

may remain2. However, minerals are believed to have little or no effect or act as catalysts rather than reaction inhibitors81-83. Consequently, it seems unlikely that the absence of minerals can account for the rapidity of pyrolysis in molecular dynamics simulations. Chemical separation may also alter kerogen, but we are not aware of any evidence for the idea that chemical alteration during chemical separation might reduce the rate of kerogen pyrolysis. Finally, differences between an ensemble of identical kerogen monomers and a more heterogeneous cross-linked network could be important. In particular, radicals and small molecules may be much more mobile in an ensemble of monomers undergoing pyrolysis than a cross-linked polymer network undergoing pyrolysis. Although polymer networks were observed in the simulations reported here, it is likely that they are substantially more tenuous than the evolving kerogen network during maturation in the subsurface and in pyrolysis experiments. V. Conclusion Reactive force field molecular dynamics simulations were used to investigate bond breaking and formation processes during the thermal maturation of kerogens and possible crosslinking reaction pathways. The kerogen crosslinking mechanism appears to be highly complex and influenced by kerogen type and maturity. A large number of possible reaction mechanisms result in the generation of highly reactive free radicals, reactive small molecules such as ethylene, formaldehyde and formic acid and much less reactive small molecules, which can be considered to be unreactive during the simulations and natural kerogen maturation. The primary crosslinks that bind kerogen monomers together to form larger 3D macromolecules or macromolecular networks are —C—C— (alkyl—aryl and alkyl—alkyl), —C—O—, and —C—S— bonds.


29

Energy & Fuels

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 30 of 47

Although ReaxFF reactive force field simulations do have some limitations, they provide much needed fundamental insights into kerogen thermal maturation and crosslinking processes. The quasi-infinite kerogen molecular structure developed in the present study is a preliminary step towards the development of more robust 3D macromolecular kerogen structures that will facilitate a better fundamental understandings of the retention and transport of fluid molecules in kerogen matrices. However, high-temperature reactive force field molecular dynamics simulations are unable to provide quantitative information about the kinetics of natural kerogen maturation because relative reaction rates are highly temperature dependent, particularly when activation energies are significantly different. Acknowledgement This work was supported by the Laboratory Directed Research and Development (LDRD) program at the Idaho National Laboratory (INL), which is operated by the Battelle Energy Alliance for the U.S. Department of Energy under Contract No. DE-AC07-051D14517. We would also like to thank Prof. Darryl Butt and Dr. Patrick Price, Boise State University, for the scanning electron microscopy (SEM) characterization of the Eagle Ford shale sample and Prof. Adri van Duin, Pennsylvania State University, for providing the C/H/N/O/S ReaxFF force field parameters.


30

Page 31 of 47

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Energy & Fuels

References 1. Durand, B., Kerogen: Insoluble organic matter from sedimentary rocks. Editions technip: 1980. 2. Vandenbroucke, M.; Largeau, C., Kerogen origin, evolution and structure. Org. Geochem. 2007, 38 (5), 719-833. 3. Ciobanu, M.C.; Burgaud, G.; Dufresne, A.; Breuker, A.; Rédou, V.; Maamar, S. B.; Gaboyer, F.; Vandenabeele-Trambouze, O.; Lipp, J. S.; Schippers, A., Microorganisms persist at record depths in the subseafloor of the Canterbury Basin. ISME J. 2014, 8 (7), 1370-1380. 4. Magnabosco, C.; Ryan, K.; Lau, M. C.; Kuloyo, O.; Lollar, B. S.; Kieft, T. L.; van Heerden, E.; Onstott, T. C., A metagenomic window into carbon metabolism at 3 km depth in Precambrian continental crust. ISME J. 2016, 10 (3), 730-741. 5. Huang, W.; Peng, P. a.; Yu, Z.; Fu, J., Effects of organic matter heterogeneity on sorption and desorption of organic contaminants by soils and sediments. Appl. Geochem. 2003, 18 (7), 955-972. 6. Karapanagioti, H. K.; Kleineidam, S.; Sabatini, D. A.; Grathwohl, P.; Ligouis, B., Impacts of heterogeneous organic matter on phenanthrene sorption: equilibrium and kinetic studies with aquifer material. Environ. Sci. Tech. 2000, 34 (3), 406-414. 7. Xiao, B.; Yu, Z.; Huang, W.; Song, J.; Peng, P., Black carbon and kerogen in soils and sediments. 2. Their roles in equilibrium sorption of less-polar organic pollutants. Environ. Sci. Tech. 2004, 38 (22), 5842-5852. 8. Behar, F.; Vandenbroucke, M., Chemical modelling of kerogens. Org. Geochem. 1987, 11 (1), 15-24. 9. Bousige, C.; Ghimbeu, C. M.; Vix-Guterl, C.; Pomerantz, A. E.; Suleimenova, A.; Vaughan, G.; Garbarino, G.; Feygenson, M.; Wildgruber, C.; Ulm, F. J., Realistic molecular model of kerogen's nanostructure. Nat. Mater. 2016, 15 (5), 576-582. 10. Orendt, A. M.; Pimienta, I. S.; Badu, S. R.; Solum, M. S.; Pugmire, R. J.; Facelli, J. C.; Locke, D. R.; Chapman, K. W.; Chupas, P. J.; Winans, R. E., Three-dimensional structure of the Siskin Green River Oil Shale Kerogen Model: a comparison between calculated and observed properties. Energy Fuels 2013, 27 (2), 702-710. 11. Siskin, M.; Scouten, C.; Rose, K.; Aczel, T.; Colgrove, S.; Pabst Jr, R., Detailed structural characterization of the organic material in Rundle Ramsay Crossing and Green River oil shales. In Composition, geochemistry and conversion of oil shales, Springer: 1995; pp 143-158. 12. Ungerer, P.; Collell, J.; Yiannourakou, M., Molecular modeling of the volumetric and thermodynamic properties of kerogen: Influence of organic type and maturity. Energy Fuels 2014, 29 (1), 91-105. 13. Ballice, L.; Yüksel, M.; Saglam, M.; Schulz, H.; Hanoglu, C., Application of infrared spectroscopy to the classification of kerogen types and the thermogravimetrically derived pyrolysis kinetics of oil shales. Fuel 1995, 74 (11), 1618-1623. 14. Jehlička, J.; Urban, O.; Pokorný, J., Raman spectroscopy of carbon and solid bitumens in sedimentary and metamorphic rocks. Spectrochi. Acta, Part A: Mol. Biomol. Spectrosc. 2003, 59 (10), 2341-2352.


31

Energy & Fuels

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 32 of 47

15. Kelemen, S.; Afeworki, M.; Gorbaty, M.; Sansone, M.; Kwiatek, P.; Walters, C.; Freund, H.; Siskin, M.; Bence, A.; Curry, D., Direct characterization of kerogen by X-ray and solid-state 13C nuclear magnetic resonance methods. Energy Fuels 2007, 21 (3), 1548-1561. 16. Kister, J.; Guiliano, M.; Largeau, C.; Derenne, S.; Casadevall, E., Characterization of chemical structure, degree of maturation and oil potential of Torbanites (type I kerogens) by quantitative FT-ir spectroscopy. Fuel 1990, 69 (11), 1356-1361. 17. Liu, X.; Zhan, J. H.; Lai, D.; Liu, X.; Zhang, Z.; Xu, G., Initial pyrolysis mechanism of oil shale kerogen with reactive molecular dynamics simulation. Energy Fuels 2015, 29 (5), 29872997. 18. Bushnev, D.; Burdel’naya, N.; Mokeev, M.; Gribanov, A. In Chemical structure and 13C NMR spectra of the kerogen of carbonaceous rock masses, Doklady Earth Sciences, Springer: 2010; pp 210-213. 19. Carlson, G., Computer simulation of the molecular structure of bituminous coal. Energy Fuels 1992, 6 (6), 771-778. 20. Kribii, A.; Lemée, L.; Chaouch, A.; Amblès, A., Structural study of the Moroccan Timahdit (Y-layer) oil shale kerogen using chemical degradations. Fuel 2001, 80 (5), 681-691. 21. Thomas, J. J.; Valenza, J. J.; Craddock, P. R.; Bake, K. D.; Pomerantz, A. E., The neutron scattering length density of kerogen and coal as determined by CH3OH/CD3OH exchange. Fuel 2014, 117, 801-808. 22. Falk, K.; Coasne, B.; Pellenq, R.; Ulm, F. J.; Bocquet, L., Subcontinuum mass transport of condensed hydrocarbons in nanoporous media. Nat. Commun. 2015, 6. 23. Pawar, G. Modeling and simulation of the pore-scale multiphase fluid transport in shale reservoirs: A molecular dynamics simulation approach. Ph.D. Dissertation, The University of Utah, Salt Lake City, UT, 2016. 24. Zhang, L.; LeBoeuf, E. J., A molecular dynamics study of natural organic matter: 1. Lignin, kerogen and soot. Org. Geochem. 2009, 40 (11), 1132-1142. 25. Hu, Y.; Devegowda, D.; Striolo, A.; Phan, A.; Ho, T. A.; Civan, F.; Sigal, R., The dynamics of hydraulic fracture water confined in nano-pores in shale reservoirs. J. Unconvent. Oil Gas Res. 2015, 9, 31-39. 26. Pignatello, J. J., Dynamic interactions of natural organic matter and organic compounds. J. Soils Sediment. 2012, 12 (8), 1241-1256. 27. Schneckenburger, T.; Lattao, C.; Pignatello, J. J.; Schaumann, G. E.; Thiele-Bruhn, S.; Cao, X.; Mao, J., Preparation and characterization of humic acid cross-linked with organic bridging groups. Org. Geochem. 2012, 47, 132-138. 28. Shokuhfar, A.; Arab, B., The effect of cross linking density on the mechanical properties and structure of the epoxy polymers: molecular dynamics simulation. J. Mol. Model. 2013, 19 (9), 3719-3731. 29. Freund, H.; Walters, C. C.; Kelemen, S. R.; Siskin, M.; Gorbaty, M. L.; Curry, D. J.; Bence, A., Predicting oil and gas compositional yields via chemical structure–chemical yield modeling (CS-CYM): Part 1–Concepts and implementation. Org. Geochem. 2007, 38 (2), 288-305.


32

Page 33 of 47

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Energy & Fuels

30. 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. Micropor. Mesopor. Mater. 2014, 197, 194-203. 31. Falk, K.; Pellenq, R.; Ulm, F. J.; Coasne, B., Effect of Chain Length and Pore Accessibility on Alkane Adsorption in Kerogen. Energy Fuels 2015, 29 (12), 7889-7896. 32. Gensterblum, Y.; Busch, A.; Krooss, B. M., Molecular concept and experimental evidence of competitive adsorption of H2O, CO2 and CH4 on organic material. Fuel 2014, 115, 581-588. 33. Hruljova, J.; Savest, N.; Oja, V.; Suuberg, E. M., Kukersite oil shale kerogen solvent swelling in binary mixtures. Fuel 2013, 105, 77-82. 34. Kelemen, S.; Walters, C.; Ertas, D.; Kwiatek, L.; Curry, D., Petroleum expulsion Part 2. Organic matter type and maturity effects on kerogen swelling by solvents and thermodynamic parameters for kerogen from regular solution theory. Energy Fuels 2006, 20 (1), 301-308. 35. Savest, N.; Oja, V.; Kaevand, T.; Lille, Ü., Interaction of Estonian kukersite with organic solvents: A volumetric swelling and molecular simulation study. Fuel 2007, 86 (1), 17-21. 36. Van Duin, A. C.; Dasgupta, S.; Lorant, F.; Goddard, W. A., ReaxFF: a reactive force field for hydrocarbons. J. Phys. Chem. A 2001, 105 (41), 9396-9409. 37. Castro-Marcano, F.; Kamat, A. M.; Russo, M. F.; van Duin, A. C.; Mathews, J. P., Combustion of an Illinois No. 6 coal char simulated using an atomistic char representation and the ReaxFF reactive force field. Combust. Flame 2012, 159 (3), 1272-1285. 38. Liu, L.; Bai, C.; Sun, H.; Goddard III, W. A., Mechanism and kinetics for the initial steps of pyrolysis and combustion of 1, 6-dicyclopropane-2, 4-hexyne from ReaxFF reactive dynamics. J. Phys. Chem. A 2011, 115 (19), 4941-4950. 39. Salmon, E.; van Duin, A. C.; Lorant, F.; Marquaire, P. M.; Goddard, W. A., Early maturation processes in coal. Part 2: Reactive dynamics simulations using the ReaxFF reactive force field on Morwell Brown coal structures. Org. Geochem. 2009, 40 (12), 1195-1209. 40. Chenoweth, K.; Cheung, S.; Van Duin, A. C.; Goddard, W. A.; Kober, E. M., Simulations on the thermal decomposition of a poly (dimethylsiloxane) polymer using the ReaxFF reactive force field. J. Am. Chem. Soc. 2005, 127 (19), 7192-7202. 41. Mueller, J. E.; van Duin, A. C.; Goddard III, W. A., Development and validation of ReaxFF reactive force field for hydrocarbon chemistry catalyzed by nickel. J. Phys. Chem. C 2010, 114 (11), 4939-4949. 42. Chenoweth, K.; Van Duin, A. C.; Goddard, W. A., ReaxFF reactive force field for molecular dynamics simulations of hydrocarbon oxidation. J. Phys. Chem. A 2008, 112 (5), 1040-1053. 43. Cheng, B.; Zhao, J.; Yang, C.; Tian, Y.; Liao, Z., Geochemical Evolution of Occluded Hydrocarbons inside Geomacromolecules: A Review. Energy Fuels 2017. 44. Zhang, Z.; Jamili, A. In Modeling the Kerogen 3D Molecular Structure, SPE/CSUR Unconvent. Res. Conf., Society of Petroleum Engineers: 2015.


33

Energy & Fuels

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 34 of 47

45. Jiménez, A.; Iglesias, M. J.; Laggoun-Défarge, F.; Suarez-Ruiz, I., Study of physical and chemical properties of vitrinites. Inferences on depositional and coalification controls. Chem. Geo. 1998, 150 (3), 197-221. 46. Quigley, T.; Mackenzie, A., The temperatures of oil and gas formation in the sub-surface. Nature 1988, 333 (6173), 549-552. 47. Burdige, D. J., Preservation of organic matter in marine sediments: controls, mechanisms, and an imbalance in sediment organic carbon budgets? Chem. Rev. 2007, 107 (2), 467-485. 48. Burnham, A. K., A simple kinetic model of oil generation, vaporization, coking, and cracking. Energy Fuels 2015, 29 (11), 7156-7167. 49. Ungerer, P., State of the art of research in kinetic modelling of oil formation and expulsion. Org. Geochem. 1990, 16 (1-3), 1-25. 50. Behar, F.; Vandenbroucke, M.; Tang, Y.; Marquis, F.; Espitalie, J., Thermal cracking of kerogen in open and closed systems: determination of kinetic parameters and stoichiometric coefficients for oil and gas generation. Org. Geochem. 1997, 26 (5), 321-339. 51. Pepper, A. S.; Corvi, P. J., Simple kinetic models of petroleum formation. Part I: oil and gas generation from kerogen. Marine Petro. Geo. 1995, 12 (3), 291-319. 52. Burnham, A. K.; Happe, J. A., On the mechanism of kerogen pyrolysis. Fuel 1984, 63 (10), 1353-1356. 53. Damsté, J. S. S.; Eglinton, T. I.; De Leeuw, J. W.; Schenck, P., Organic sulphur in macromolecular sedimentary organic matter: I. Structure and origin of sulphur-containing moieties in kerogen, asphaltenes and coal as revealed by flash pyrolysis. Geochim. Cosmochim. Acta 1989, 53 (4), 873-889. 54. González-Vila, F. J.; Amblès, A.; del Rıo, J.; Grasset, L., Characterisation and differentiation of kerogens by pyrolytic and chemical degradation techniques. J. Anal. Appl. Pyro. 2001, 58, 315328. 55. Witte, E.; Schenk, H.; Müller, P.; Schwochau, K., Structural modifications of kerogen during natural evolution as derived from 13C CP/MAS NMR, IR spectroscopy and Rock-Eval pyrolysis of Toarcian shales. Org. Geochem. 1988, 13 (4), 1039-1044. 56. Plimpton, S., Fast parallel algorithms for short-range molecular dynamics. J. Comput. Phys. 1995, 117 (1), 1-19. 57. Aktulga, H. M.; Fogarty, J. C.; Pandit, S. A.; Grama, A. Y., Parallel reactive molecular dynamics: Numerical methods and algorithmic techniques. Parallel Comput. 2012, 38 (4), 245259. 58. Blanksby, S. J.; Ellison, G. B., Bond dissociation energies of organic molecules. Acc. Chem. Res. 2003, 36 (4), 255-263. 59. Kerr, J., Bond dissociation energies by kinetic methods. Chem. Rev. 1966, 66 (5), 465-500. 60. Krichko, A. A.; Gagarin, S. G., New ideas of coal organic matter chemical structure and mechanism of hydrogenation processes. Fuel 1990, 69 (7), 885-891. 61. Allen, D. T.; Gavalas, G. R., Reactions of methylene and ether bridges. Fuel 1984, 63 (5), 586592.


34

Page 35 of 47

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Energy & Fuels

62. Richnow, H.; Jenisch, A.; Michaelis, W., Structural investigations of sulphur-rich macromolecular oil fractions and a kerogen by sequential chemical degradation. Org. Geochem. 1992, 19 (4-6), 351-370. 63. Faulon, J.; Vandenbroucke, M.; Drappier, J.; Behar, F.; Romero, M., 3D chemical model for geological macromolecules. Org. Geochem. 1990, 16 (4-6), 981-993. 64. Leszczynski, J., Handbook of computational chemistry. Springer Science & Business Media: 2012; Vol. 2. 65. Bakr, M.; Akiyama, M.; Sanada, Y., ESR assessment of kerogen maturation and its relation with petroleum genesis. Org. Geochem. 1990, 15 (6), 595-599. 66. Larsen, J. W.; Ashida, R.; Doetschman, D. C., Kerogen chemistry 6. Involvement of native free radicals in kerogen thermal deoxygenation. Energy Fuels 2005, 19 (5), 2216-2217. 67. Bustin, R.; Ross, J.; Rouzaud, J.-N., Mechanisms of graphite formation from kerogen: experimental evidence. Int. J. Coal Geo. 1995, 28 (1), 1-36. 68. Gupta, N. S.; Briggs, D. E.; Landman, N. H.; Tanabe, K.; Summons, R. E., Molecular structure of organic components in cephalopods: Evidence for oxidative cross linking in fossil marine invertebrates. Org. Geochem. 2008, 39 (10), 1405-1414. 69. Parikh, H.; Larsen, J., Macromolecular structure of Type I and II kerogen: Advances in hydrocarbon characterization. Prep-Am. Chem. Soc. Div. Petro Chem. 2000, 45 (4), 556-558. 70. Van Krevelen, D., Graphical-statistical method for the study of structure and reaction processes of coal. Fuel 1950, 29 (12), 269-284. 71. Flory, P. J.; Rehner Jr, J., Statistical mechanics of cross-linked polymer networks I. Rubberlike elasticity. J. Chem. Phys. 1943, 11 (11), 512-520. 72. Neuburger, N.; Eichinger, B., Critical experimental test of the Flory-Rehner theory of swelling. Macromolecules 1988, 21 (10), 3060-3070. 73. Martin Jr, D. L. Crosslink Density Determinations for Polymeric Materials; DTIC Document: 1970. 74. Tillet, G.; Boutevin, B.; Ameduri, B., Chemical reactions of polymer crosslinking and postcrosslinking at room and medium temperature. Prog. Poly. Sci. 2011, 36 (2), 191-217. 75. Tonelli, A. E., NMR spectroscopy and polymer microstructure: the conformational connection. Titles Supplied by John Wiley & Sons Australia: 1989. 76. Han, Y.; Jiang, D.; Zhang, J.; Li, W.; Gan, Z.; Gu, J., Development, applications and challenges of ReaxFF reactive force field in molecular simulations. Front. Chem. Sci. Eng. 2016, 10 (1), 1638. 77. Salmon, E.; van Duin, A. C.; Lorant, F.; Marquaire, P.-M.; Goddard, W. A., Thermal decomposition process in algaenan of Botryococcus braunii race L. Part 2: Molecular dynamics simulations using the ReaxFF reactive force field. Org. Geochem. 2009, 40 (3), 416-427. 78. Wang, H.; Feng, Y.; Zhang, X.; Lin, W.; Zhao, Y., Study of coal hydropyrolysis and desulfurization by ReaxFF molecular dynamics simulation. Fuel 2015, 145, 241-248.


35

Energy & Fuels

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 36 of 47

79. Zou, C.; Raman, S.; van Duin, A. C., Large-scale reactive molecular dynamics simulation and kinetic modeling of high-temperature pyrolysis of the gloeocapsomorphaprisca microfossils. J. Phys. Chem. B 2014, 118 (23), 6302-6315. 80. Bauschlicher Jr, C. W.; Qi, T.; Reed, E. J.; Lenfant, A.; Lawson, J. W.; Desai, T. G., Comparison of ReaxFF, DFTB, and DFT for phenolic pyrolysis. 2. Elementary reaction paths. J. Phys. Chem. A 2013, 117 (44), 11126-11135. 81. Dembicki, H., The effects of the mineral matrix on the determination of kinetic parameters using modified Rock Eval pyrolysis. Org. Geochem. 1992, 18 (4), 531-539. 82. Mango, F. D., The light hydrocarbons in petroleum: a critical review. Org. Geochem. 1997, 26 (7-8), 417-440. 83. Yang, S.; Horsfield, B., Some predicted effects of minerals on the generation of petroleum in nature. Energy Fuels 2016, 30 (8), 6677-6687.


36

Page 37 of 47

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Energy & Fuels

Figures

Figure 1. A scanning electron microscope (SEM) image of Eagle Ford shale. The disparity in the pore sizes and their distributions within kerogen particles can be clearly observed. Moreover, kerogen particle 1 (at the bottom right) appears to be denser, and it does not show any clearly visible pores when compared with kerogen particle 2, which shows abundant pores of various sizes. Such differences in kerogen structure and kerogen morphology in kerogen particles that are separated by only about1 m or less provide a clear indication of the structural complexity of kerogen and its chemical heterogeneity.


37

Energy & Fuels

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 38 of 47

Figure 2. The workflow used to perform the reactive force field molecular dynamics simulations. The final altered kerogen was analyzed to identify various bond cleavage and bond formation reactions, to determine the crosslinks formed between atoms that were initially located in different monomers, and to identify the secondary products generated during the thermochemical transformation of kerogen.


38

Page 39 of 47

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Energy & Fuels

Figure 3. The kerogen monomers proposed by Ungerer et al.12 to represent: a) kerogen 1-I with elemental composition C251H385O13N7S3, b) kerogen 2-L with elemental composition C234H263O14N5S2, and c) kerogen 3-L with elemental composition C233H204O27N4. The gray, white, red, blue, and yellow colors represent carbon (C), hydrogen (H), oxygen (O), nitrogen (N), and sulfur (S) atoms respectively. Additional details concerning the kerogen monomers can be found elsewhere12.


39

Energy & Fuels

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 40 of 47

4a. Kerogen 1-I bond cleaving sites (shown by dashed lines) i. Aliphatic carbon—aliphatic carbon bond

iii. Aliphatic carbon—oxygen bond (ether)

ii. Aliphatic carbon—aromatic carbon bond

iv. Aliphatic carbon—sulfur bond (thioether)

4b. Kerogen 2-L bond cleaving sites (shown by dashed lines) i. Aromatic carbon—hydrogen bond

iii. Aliphatic carbon—aromatic carbon bond

v. carbon  to oxygen in furanic group— hydrogen bond

ii. Aliphatic carbon—aliphatic carbon bond

iv. Aliphatic carbon—oxygen bond

vi. Aliphatic carbon—sulfur bond

vii. Aliphatic carbon—oxygen bond (ether)

4c. Kerogen 3-L bond cleaving sites (shown by dashed lines) i. Aliphatic carbon—aliphatic carbon bond (in short aliphatic chains)


40

Page 41 of 47

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Energy & Fuels

ii. Aliphatic carbon (in short aliphatic carbon chain)—aliphatic carbon bond

iii. Aliphatic carbon connecting aromatic iv. Aliphatic carbon—aliphatic carbon and aliphatic rings bond

v. Aromatic carbon—oxygen bond

vi. Aliphatic carbon—oxygen bond

vii. Aliphatic carbon—oxygen bond

Figure 4. Identification of bond cleavage sites (shown by dashed lines) during reactive molecular dynamics simulations. A variety of possible reaction pathways based on the cleavage of aliphatic carbon-aliphatic carbon, thioether, and ethers bonds that generate highly reactive radicals that subsequently participate in bond forming reactions are shown. Atoms are color coded for visual distinction. (a) kerogen 1-I bond cleaving sites, (b) kerogen 2L bond cleaving sites, (c) kerogen 3-L bond cleaving sites. The moieties highlighted in gray become free radicals stabilized by electron delocalization after bond cleavage.


41

Energy & Fuels

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 42 of 47

Figure 5. A variation in the number of bonds cleaved and bonds formed in a) kerogen 1-I, b) kerogen 2L, and c) kerogen 3-L during the (P, T) schedule given in Table 1. The net bonds formed is the difference between the number of bonds cleaved and number of bonds formed. Further, each colored compartment corresponds to a simulation stage from simulation schedule described in Table 1.


42

Page 43 of 47

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Energy & Fuels

Figure 6. Perspective view of the final crosslinked kerogen structures. These structures are connected to themselves through the periodic boundaries used in the simulations, and they can be considered as the unit cells of the quasi-infinite cubic network. The gray, blue, red and yellow colors represent carbon (C), nitrogen (N), oxygen (O), and sulfur (S) atoms respectively. The locations of only newly formed crosslinks among kerogen monomers are shown as balls for the clarity. a) kerogen 1-I with the —C—S— crosslinking bonds, b) kerogen 2-L with —C—S— and —C—C— crosslinking bonds, and c) kerogen 3L with —C—C—crosslinking bonds. Both alkyl-alkyl and alkyl-aryl carbon-carbon crosslinks were formed.


43

Energy & Fuels

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 44 of 47

Figure 7. The carbon-carbon pairwise distribution function g(r) in the residual crosslinked kerogen 2-L structure. The peak location of g(r) are 1.41 Å (aromatic carbons–aromatic carbon bond with shorter bond length), 1.49 Å and 1.55 Å (aliphatic carbon–aliphatic carbon and aliphatic-aromatic single carbon-carbon bonds with comparatively larger bond length). The g(r) peak shapes are not well defined and cannot be separated but the peak area at 1.41 Å appears larger than the peak areas at 1.49 Å and 1.55 Å thus indicating that the residual cross-linked kerogen structure contains more aromatic-aromatic carbon-carbon bonds than that other carbon-carbon bonds.


44

Page 45 of 47

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Energy & Fuels

Table 1. Pressure – temperature (P, T) schedule for each kerogen ensemble. All simulations were performed using isobaric-isothermal (NPT) ensembles. A time step of 0.25 fs with a total simulation time of 1.025 ns was used to perform the reactive molecular dynamics simulations. The pressure (P) and the temperature (T) are specified in units of standard atmospheres (~0.1 MPa) and K respectively.

Stage

Start (P, T)

End (P, T)

Simulation time for each stage (ns)

1

1, 300

100, 300

0.025

2

100, 300

100, 1400

0.2

3

100, 1400

100, 1300

0.2

4

100, 1300

100, 1200

0.2

5

100, 1200

100, 1000

0.2

6

100, 1000

100, 800

0.1

7

100, 800

100, 500

0.05

8

100, 500

1, 300

0.05

(P, T) simulation schedule

Total simulation 1.025 time (ns)


45

Energy & Fuels

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 46 of 47

Table 2. Secondary reaction products generated during the thermal evolution of kerogen (at the end of the simulation schedule described in Table 1).

Molecule/

Number of molecules/fragments

Fragment

kerogen 1-I

kerogen 2-L

kerogen 3-L

C2H4

3

5

1

C3H6

0

2

1

C4H8

1

0

0

C4H6

0

0

1

CH3OH

2

1

0

CH2CHOH

0

0

2

H2O

5

6

24

H3O

0

0

1

CH2O

0

1

3

CH3CHCH3

1

0

0

CH2OCH2

0

0

1

CH3CHOH

0

0

1

CH2OH

0

0

1

OCHOH

1

3

1

CH2CHCHOH 0

0

1

OHCHOH

0

0

3

OCHCH2OH

0

0

1

OH

0

0

1


46

Page 47 of 47

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Energy & Fuels

Table 3. Comparison of the elemental composition of the starting kerogen 2-L ensemble and the final residual crosslinked kerogen 2-L structure. The elemental composition of the starting kerogen structure is taken from Ungerer et al.12, whereas the elemental composition of the final crosslinked kerogen structure was determined from the residual kerogen 2-L.

Elemental composition

Starting kerogen Structure

Residual Crosslinked kerogen structure

H/C

1.1200

1.10670

O/C

0.0600

0.05457

N/C

0.0220

0.02154

S/C

0.0090

0.00860


47

Reactive Molecular Dynamics Simulation of Kerogen Thermal

Recommend Documents