Early maturation processes in coal. Part 2

Early maturation processes in coal.Part 2:Reactive dynamics simulations using the ReaxFF reactive force field on Morwell Brown coal structuresElodie Salmon a ,Adri C.T.van Duin b ,François Lorant a ,Paul-Marie Marquaire c ,William A.Goddard III d,*aInstitut Français du Pétrole,BP 311,92506Rueil-Malmaison Cedex,FrancebDepartment of Mechanical and Nuclear Engineering,Penn State University,University Park,PA 16802,USA cDépartement de Chimie Physique des Réactions (DCPR),Nancy-Université,CNRS 1rue Grandville,BP 20451,F-54001Nancy,France dMaterials and Process Simulation Center (MC 139-74),California Institute of Technology,Pasadena,CA 91125,USAa r t i c l e i n f o Article history:Received 20April 2009Received in revised form 31August 2009Accepted 2September 2009Available online 6September 2009a b s t r a c tThis paper reports reactive dynamics (RD)simulations of a macro-model of Morwell Brown coal using the ReaxFF reactive force field.We find that these reactive MD simulations successfully reproduce thermal decomposition processes of defunctionalization,depolymerization and rearrangement of the residual structure observed in various experimental studies.For example,our simulations indicate that the decar-boxylation and dehydroxylation of the lignin side chain of the Morwell model involves the formation of double bonds conjugated with the aromatic rings.The process of defunctionalization of the methoxy functions involving the formation of phenolic structures in the residue has been confirmed.We also observe that gaseous hydrocarbons are generated by cleavage of C–C bonds of the lignin side chain.The success in using ReaxFF RD to describe the molecular processes underlying the kinetics in pyrolysis of this model of coal plus the success of a similar previous study on the algaenan of Botryococcus braunii race L biopolymer model of kerogens suggests that such computation can be useful in providing molec-ular based kinetic models for other pyrolysis processes underlying the organic transformations in sedi-mentary materials.Ó2009Elsevier Ltd.All rights reserved.1.IntroductionThere is great interest in determining a kinetic model to de-scribe oil and gas generation from the thermal decomposition of fossil organic matter and of shale and coal.Such information can be invaluable in assessing migration pathways of these precursors evolving to produce oil and gas that would help to locate the opti-mum locations for extraction.This has led to a number of experi-mental studies (Huck and Karweil,1955;Karweil,1955;Vassoevich et al.,1969;Lopatin,1971;Waples,1980)to under-stand and estimate the chemical processes and the kinetics of oil and gas generation.However,it is difficult to reproduce the low temperature,high pressure conditions in practical laboratory experiments.This has made it imperative to derive computational models to predict oil generation.One of the earliest such model by Pitt (1962)consists of a set of global reactions assumed to occur simultaneously and assumed to be independent of each other.Ki-netic equations of all prototypical reactions were then definedexperimentally for each source rock.This allowed the bulk rate of oil and gas generation for a source rock to be accounted -ter,experimental thermal decomposition studies of model com-pounds and structural models of kerogen and coal were used to obtain a molecular level description of the formation of petroleum (Solomon et al.,1988;Savage and Klein,1989;Freund,1992;Payne and Ortoleva,2002;Kelemen et al.,2004;Freund et al.,2007).In addition,Solomon et al.(1988)used a structural model of coal to reproduce the devolatilization of coal observed in pyrolysis exper-iments.Other approaches have simulated the defunctionalization processes generating gas compounds and the bond breaking/cross linking processes involved in tar formation.These molecular mod-els explain defunctionalization of coal by combining the functional group composition of the coal with kinetic coefficients of each gas generated.Alternatively,the formation of the tar has been de-scribed using a statistical distribution of kinetically constrained,one step ter,statistical models of the transformation time for thermal decomposition and pyrolysis and natural evolu-tion of sedimentary organic matter were proposed for asphalthene (Savage and Klein,1989),lignite (Payne and Ortoleva,2002)and kerogen (Freund,1992;Kelemen et al.,2004;Freund et al.,2007).These previous approaches to developing chemical models have necessarily used a limited number of chemical steps involved in such pyrolysis processes,which necessarily limits the accuracy of0146-6380/$-see front matter Ó2009Elsevier Ltd.All rights reserved.doi:10.1016/geochem.2009.09.001*Corresponding author.Address:Materials and Process Simulation Center (MC 139-74),California Institute of Technology,1200East California Blvd.,Pasadena,CA 91125,USA.Tel.:+16263952731;fax:+16265850918.E-mail addresses:wag@ ,shirley@ (W.A.God-dard).Organic Geochemistry 40(2009)1195–1209Contents lists available at ScienceDirectOrganic Geochemistryjournal homepage:www.elsevie r.c o m /l o c a te /o r g g e o c h emthe predictions from these models.For this reason,we embarked on a collaborative experimental and theory-computational program to provide a more accurate molecular based model of the processes.Thefirst phase of this program reported the results of both experiment(Salmon et al.,2009a)and ReaxFF reactive dynamics computational modeling(Salmon et al.,2009b)on ther-mal decomposition process in algaenan of Botryococcus braunii race L.This current paper is the theory part of a second series of similar experimental and theory-computational studies,but now on lig-nite,a model system for coal.As in Salmon et al.(2009b)we use the ReaxFF reactive force field(van Duin et al.,2001)in reactive dynamics(RD)simulations for realistic models of lignite.[For other reactive forcefield ap-proaches to organic geochemical polymers see Stoliarov et al. (2004)and Yin et al.(2006).]The ReaxFF(van Duin et al.,2001) reactive forcefield has been demonstrated to retain nearly the accuracy of the QM calculations from which it is derived while en-abling RD at computational costs nearly as low as for ordinary non-reactive forcefields.The accuracy and efficacy of ReaxFF RD has been demonstrated for a wide variety of materials including tran-sition metal oxides and alloys,homogenous and heterogeneous molecular catalysts,and semiconductor oxides interactions(e.g., van Duin and Sinninghe Damsté,2003;Strachan et al.,2003; Chenoweth et al.,2005,2008;Chen et al.,2005;van Duin et al., 2005;Han et al.,2005;Buehler et al.,2006;Goddard et al., 2006a,b;Ludwig et al.,2006;Leininger et al.,2008).ReaxFF RD has been applied to simulations of105–106million atoms(using massively parallel computers)and with a single CPU is quite prac-tical for simulations with up to4000atoms per periodic cell(Nom-ura et al.,2007;Nakano et al.,2008).The results presented here used just a single CPU.Our previous ReaxFF RD study of the ther-mal decomposition of an aliphatic biopolymer of algaenan of B. braunii race L led to excellent comparison with experiment(Sal-mon et al.,2009a)as does the current application to lignite.The present paper uses ReaxFF RD to describe the initial steps of the thermal decomposition of the kerogen aromatic geopolymer. These simulations aim to predict the molecular decomposition of the Morwell brown coal model defined in parallel experimental study by Salmon et al.(2009c).Thefirst section describes the sim-ulation system,a molecular structure containing500atoms made up from representatives of the structural units of the coal model (e.g.,guaiacyl,syringyl and p-hydroxyl units).This allows the ther-mal decomposition processes of these specific structures and func-tional groups to be described simultaneously.Based on these results we analyze the reproducibility of the thermal decomposi-tion mechanism as a function of the initial spatial positions of the molecules.Then we present the ReaxFF RD results for a macro-model containing2692atoms of Morwell coal.Wefind that the results agree qualitatively with experimental data from the lit-erature,validating the methodology.2.Modeling method2.1.Selection of initial molecular modelsThe Morwell coal is lignite enriched in aromatic chemical struc-tures.This structure,which is similar to kerogen,is selected for the present study.A molecular model of the immature lignite samples was proposed in Salmon et al.(2009c).This model,displayed in Fig.1,was con-strained by(1)experimental characterization(elemental analysis,1196 E.Salmon et al./Organic Geochemistry40(2009)1195–1209Infrared spectra,modern1-and2-D Nuclear Magnetic Resonance)of the brown coal sample from Morwell Open Cut(Gippsland Forma-tion,Australia)and(2)the theoretical transformation of an angio-spermous lignin structure found in the literature(Nimz,1974).To set up the simulation procedure for this type of complex organic matter,wefirst selected small molecular models of 24–33atoms to represent specific substructures of the lignite macro-model(Fig.1).The lignite structure is composed of various aromatic units containing a p-hydroxyl aromatic ring substituted by a C1–C6aliphatic side chain and oxygenated functional groups (methoxy,carboxyl or carbonyl functional groups).Hence,the structural models selected were representative of the different aro-matic units observed in the macro-models of lignite(Fig.2B–F).2.2.Molecular dynamics modeling using the reactive forcefield, ReaxFFA brief description of the ReaxFF reactive forcefield and of the simulation procedure(Fig.3)is presented below;a more detailed description is available in Salmon et al.(2009b).E.Salmon et al./Organic Geochemistry40(2009)1195–12091197For the present study,the simulations aimed at reproducing the experimental results of off line pyrolysis of the Morwell lignite per-formed by Salmon et al.(2009c).To initiate the simulations a peri-odic cubic box was built containing the molecular models described in the previous section.Each box contained up to 15identical functional models,initially at a low density of 0.1kg/dm 3,resulting in system sizes of around 500atoms.We also built one macro-model of lignite containing 2692atoms at low density.Subsequently,a low temperature reactive dynamics (for functional models at 2K and the macro-model at 50K)was performed on these initial low density configurations.Thereafter,the system was equilibrated at 300K using NVT-dynamics (2.5ps =picosec-onds)and a Berendsen thermostat (Berendsen et al.,1984).This system was then compressed to the experimental solid density at room temperature of d =1.4g dm À3by means of re-scaling the cell parameters during a 1.25ps,MD-simulation at T =300K.Then the system was equilibrated at a pressure of 1bar and a temperature of 750K using NPT-dynamics (25ps)with the Berendsen thermostat and barostat.This causes an expansion of the box size to a density of about 0.7–0.8kg/dm 3for the functional models and of 0.68kg/dm 3for the macro-model (Table 1).Finally,we performed heat up simulations from 300–2500K at rates of 88,44and 22K/ps for the functional models (Fig.4)to determine the onset tempera-ture of thermal decomposition at the picosecond time scale of our simulations.For the macro-model A the onset temperature of ther-mal decomposition was determined by heat up simulation per-formed from 700–3100K for 40ps corresponding to a rate of 60K/ps.We observed that the thermal decomposition of functional models began between 1700and 2500K and that of the macro-model at around 1200K.To assess the kinetics of the thermal decomposition reactions,we also performed constant temperature molecular dynamics simulations on each functional models for 50ps at T =2000,2100and 2200K.The macro-model was isother-mally heated as well,for 40ps at 1000,1200,1400,1600,1800and 2000K.The time scale of the simulation (50picoseconds)is many mag-nitudes lower than that employed in the experiments (9hours in Salmon et al.,2009c ).Thus,for the numerical simulations we arti-ficially increased the temperature to be in the range 500–2200K,allowing the chemical reactions to occur within a few picoseconds,but hopefully with similar mechanism to laboratory experiments which take hours at much lower temperature (i.e.,473–573K).We previously demonstrated that despite these time and temper-ature differences between simulation and experiment,we obtained good qualitative agreement between the initial reaction products on aliphatic structure (Salmon et al.,2009b ).3.Simulations on the thermal decomposition of the lignite geopolymer3.1.Thermal decomposition of functional models of ligniteFig.5gives an overview of the kinetics of degradation of each functional model for NVT-dynamic simulations at 2000,2100and 2200K for 50ps.As expected the chemical systems tend to be more reactive at higher temperature.We do reversible reactions (noted R in Fig.5)involving either addition reactions or reverse dismutations involving hydrogen exchanges.For example:ðB ÞC 12H 16O 2þC 8H 8O 2¢C 12H 17O 12þC 8H 7O 2ðC ÞC 10H 13O 5þC 10H 11O 5¢2C 10H 12O 5ðD ÞC 11H 16O 3þCH 3¢C 11H 15O 3þCH 4ðE Þ2C 9H 10O 2¢C 18H 20O 4¢C 9H 11O 2þC 9H 9O 2ðF Þ2C 10H 12O 2¢C 20H 24O 4The aliphatic cycle (F),methoxy (D)and acid (C)functional models decompose faster than the hydroxyl (B)and the ketone (E)models.The acid (C)model leads to earlier initiation reactions compared to the aliphatic cycle (F)than the methoxy (D)models.Table 1Simulation conditions for constant temperature NVT RD simulations and the final results.Models Numbers Volume (Å3)Pressure (Mpa)Initial density (kg/dm 3)Final density (kg/dm 3)Nature Structures Number atoms Number cases Macro-model A Lignite 269215536010 1.430.68Functional modelB Hydroxyl 3315655010 1.430.80C Acid2718900910 1.430.70D Methoxy 3017703710 1.430.79E Ketone2124698510 1.430.86FAliphatic cycle24217159101.430.801198 E.Salmon et al./Organic Geochemistry 40(2009)1195–1209In Tables 2and 3,the products of decomposition are partitioned into four classes of compounds.The aromatic compounds contain aromatic rings and larger struc-tures that have no more carbons than their respective initial function models.Those compounds represent residual struc-tures that are resistant to thermal stress.The alkyl compounds are structural moieties with less than four carbons and derived from the decomposition of the aliphatic side chains.This fraction is composed of oxygenated compounds like formaldehyde,ethanol,propanol and small hydrocarbons like butane,propene,ethene,methyl radical and methane.Other compounds represent inorganic compounds like water,carbon dioxide or H 2.Condensed compounds in Table 3,contain more than the ten car-bons of the initial acid (C 10H 12O 5)and aliphatic cycle functional models (C 10H 12O 2).The first main steps deduced from the reaction products are dis-played in Figs.6–8.Those reactions begin with an initiation reac-tion followed by beta-scission decomposition and/or methathesis (H transfer)with the reactant.For each initiation the bond dissoci-ation energy (BDE)proposed by Luo (2003)is given indicating the strength of the bond that is cleaved and so giving an idea of the occurrence of each reaction.For all the functional models we ob-serve decomposition of the aliphatic side chain and of other func-tional groups,but no aromatic ring is degradated,as expected.The specific decomposition of the aliphatic cycle functional model,dis-played as reaction F in Fig.8,shows that depolymerization and condensation reactions may evolve simultaneously.The generation of unstable biradicals led to the strong rearrangements in the car-bon skeleton structure like those displayed in Fig.9.The occur-rence of high molecular weight compounds resulted from addition or termination reactions.However those reactions are less likely to occur when the aliphatic cycle is substituted by functional groups or connected to other aromatic units or side chains,as is common in lignite macro-structures.For all functional models,we observed radicals generated via a depolymerization of the ali-phatic side chain.In terms of defunctionalization,we observed dehydroxylation of both alkyl and aryl hydroxyl,i.e.,reaction B3and C (Figs.6and 7),releasing both hydroxyl radicals and water.Sometimes,b -hydroxyl functions were decomposed into such oxy-genated compounds as formaldehyde,ethanol and propanol.Decarboxylations releasing carbon dioxide were obtained by the decomposition of the acid functional model C (Fig.7).It is notice-able that the ketone group (reaction E in Fig.7)exhibits a high resistance to thermal stress.This stability is increased by the ben-zyl position of the ketone group in the functional model.Reaction D reports the demethoxylation process generating methyl and phenolic radicals.Thereafter,hydrogenation reactions were at the origin of methane and of the formation of catechol moieties.3.2.Reproducibility of functional model resultsFor each functional model,as defined in Fig.2,we constructed periodic boxes containing 15–24duplicates,depending on the510152025102030405060O H OH O H OOH C H 32100K 8/182000K 3/182200K 10/185101520250102030405060O HOC H 3C H 32100K 7/242000K 0/242200K 4/240510152025102030405060OH O HC H 3O HC H 32100K 0/152000K 1/152200K 7/15510152025O HOH H 3COCH 3C H 32100K 7/172000K 2/172200K 15/175101520250102030405060102030405060Time (ps)Time (ps)OH OH 2100K 15/212000K 8/212200K 18/21CDBEFRRRRRRRRRRRRN u m b e r o f m o l e c u l e sN u m b e r o f m o l e c u l e sN u m b e r o f m o l e c u l e sN u m b e r o f m o l e c u l e sN u m b e r o f m o l e c u l e sFig.5.Thermal of decomposition of functional models for NVT RD simulations at 2000,2100and 2200K.E.Salmon et al./Organic Geochemistry 40(2009)1195–12091199model size(Table1).For each of these boxes thermal decomposi-tion simulations were performed starting fromfive independent initial configurations,thus allowing us to evaluate the reproduc-ibility of the observed chemical events.Results of these reproduc-ibility experiments are displayed in Fig.10for the decomposition and the total production of the hydroxyl functional model at three temperatures(2000,2100and2200K).Similar trends were observed for both the decomposition and the products generation at2000and2100K.At these tempera-tures,respectively up to4and7initial models are decomposed generating up to10and11new molecules.At2200K,the reaction rate is higher than at lower temperature.About5to10hydroxyl functional units are decomposed into11–21new molecules.A de-tailed list of compounds generated after50ps at2000,2100and 2200K,also displayed in Fig.8,indicates that the overall reaction always creates the same products whatever the temperature.At 2000and2100K,some compounds are formed by recombination or addition reactions,either between two initial reactantsTable2Table3List of compounds obtained from thermal decomposition of functional models bearing carboxyl group and an aliphatic cycle at2000,2100and2200K during50ps of NVT RD1200 E.Salmon et al./Organic Geochemistry40(2009)1195–1209(C24H36O6),or between an initial reactant and an intermediate product,or between two intermediate products(C24H35O5, C18H20O4,C14H22O2,C13H21O4).These MD results establish molecu-lar balances for the early decomposition of the hydroxyl structure (Table4).We found thatfive bonds are especially weak in this structure,generating butyl,methyl,hydroxyl,methylhydroxyl rad-icals and propene molecules.This means that the dispersion of products observed previously in Fig.10,is controlled by the prob-ability of occurrence of the reaction as a function of the initial spa-tial conformation of the system.Detailed analysis of products formed after50ps shows that similar early reactions occur,irre-spective of the initial spatial conformation.From these reproducibility experiments,we conclude that the RD simulations performed here provide qualitative descriptions of the thermal decomposition of the organic structures studied. Quantitative results would require statistical analysis of reactions, which requires larger systems or many independent simulations for small systems.Results provided hereafter were all obtained using one spatial conformation and multiple temperatures,so they are restricted to qualitative analysis.3.3.Thermal decomposition of the macro-models of ligniteWe simulated the decomposition of the lignite macro-model A with a heat up ramp of60K/ps.Fig.11displays the formation of the main products generated during this simulation.Thefirst products formed are carbon dioxide and concomitant higher molecular weight compounds having more than three carbon atoms.Formation of CO2starts at700K after only$.10ps and higher molecular weight compounds are generated starting from 1500K at12.5ps.Carbon monoxide is observed at t=18ps corre-sponding to a temperature of$1800K.At about2200K,low molecular weight hydrocarbons,C1–C2oxygenated compounds and water appear.At the end of the simulation,4molecules of water,29molecules of C1–C3hydrocarbons and29molecules of C1–C3oxygenated compounds were generated.Thefirst hydrocar-bon molecules(methane and ethylene)appear after14ps ($1500K)followed by acetylene(C2H2)at33ps($3000K)due to the high temperature.Di-hydrogen(H2)is generated at tem-peratures higher than2500K.Molecular C1–C4oxygenated com-pounds formed at3209K/40ps contain methanol(CH4O)and formaldehyde(CH2O).Analysis of the molecular compounds formed at different time steps,listed in Table5,indicates that high molecular weight spe-cies at lower temperature(61947K)are mainly composed of C30+aromatic structures.Beyond2000K the proportion of C40Àcompounds increases.Since high temperatures favor decomposi-tion rather than addition and recombination,no additional C100+ high molecular weigh structures are observed beyond3000K.Un-der our MD simulation conditions,oxygenated compounds inE.Salmon et al./Organic Geochemistry40(2009)1195–12091201C1–C6are formed before C1–C4hydrocarbons.More CO2and methyl radicals are generated compared to CO and methane.We observe the formation of hydroxyl radicals at1947K -ter on(24ps),water appears in the system.In addition to ramped heating runs,we performed constant temperature RD simulations between500and2000K for79ps.It seemed that no reaction happened during the simulation at 500K/79ps.But closer analysis of the chemical system during the simulation shows some structural rearrangement does occur at this temperature.Indeed,two molecules of CO2are released in the system and then reincorporated into the main structure by an addition process.A few kerogen decomposition reactions gener-ate CO at1600,1800and2000K.Fig.12displays the overall forma-tion of CO2as a function of time and temperature.We observe that the defunctionalization process is dominant compared to the reac-tions of CO2.The compound classes observed at50ps in the simulations at 500,1000,1200,1400,1600,1800and2000K are displayed in pounds are lumped according to mass fractions fre-quently used in experimental petroleum studies.We observed that C100+compounds are dominant in the system.At1800K,these compounds represent62wt%of the mass of the initial structure but at2000K only29wt%of C100+compounds remain in the system.The main products formed during RD simulation are C14–C100 compounds that may be ascribed to heavy oxygenated compounds1202 E.Salmon et al./Organic Geochemistry40(2009)1195–1209E.Salmon et al./Organic Geochemistry40(2009)1195–12091203experimentally observed in the C 14+fraction extracted from the pyrolysis residue.At 1600K about as many C 6–C 14as C 1–C 6com-pounds are formed,$3–5wt%.At 2000K,C 1–C 6oxygenated com-pounds make up 3wt%of the mass of the structure,but only 0.5wt%of light hydrocarbons are generated.Formation of inorganic compounds such as CO 2,CO and water,increases with temperature to about 9wt%at 2000K.Predicted chemical structures at the end of the RD simulations at 1600,1800and 2000K show an increase in aromaticity.We also observed an increase in double bonds on the side chains conjugated with aromatic rings.Most of those dou-ble bonds are in the trans conformation.Additional ring structures are also observed in the final system configurations.Fig.14gives an overview of these ring structures.3.4.Discussion:comparison of experimental and reactive dynamics resultsOur RD simulations on small systems (<500atoms)reproduce the processes related to the formation of low molecular weight products.As observed in experiments,the RD finds that CO 2is gen-erated by defunctionalization of carboxylic groups.Initially,we thought that ketone functional groups might be responsible for CO or CO 2release;however the RD do not show any thermal decomposition of ketone functional groups.Instead we find that methyl radicals are generated by cleavage of C–C bonds on the ali-phatic side chain or by cleavage of C–O bonds of methoxy groups.Defunctionalization of the methoxy group by release of methyl radical has been proposed by Hatcher and Clifford (1997)based on observations of natural samples and pyrolysis experiments.Subsequently,the methyl radicals may convert into heavier mole-cules by addition to other molecules or radicals or by metathesis and methyl radicals can be transformed into methane.This mech-anism was strongly suggested though artificial maturation studies of Behar and Hatcher (1995)and Salmon et al.(2009c).We find that the RD confirms this mechanism for functional model D bear-ing methoxy groups.Our simulations also showed ethylene to be formed by cleavage of C–C bonds on the aliphatic side chain of the aromatic units.For example,we observe that the naphtheno-aromatic units unsubstituted present in functional models F,are able to generate a large amount of ethylene.However in lignite macro-structures it is likely that those naphtheno-aromatic units are connected to other aromatic units or side chains,which would lower the rate of ethylene generation.RD simulations of the thermal decomposition of functional model B bearing a hydroxyl group by MD show that C–C bonds in a ,b or c position of the hydroxyl group are weaker than C–O bonds involving the hydroxyl groups.Cleavage of the C–C bonds results in the formation of oxygenated compounds such as formal-dehyde,ethanol or propanol.During pyrolysis experiments (Sal-mon et al.,2009c ),these compounds would appear in the NSO fraction,which contain high molecular weight compounds.How-ever the chemical structure of the NSO fraction was not identified by Salmon et al.(2009c)for two reasons:first,the amount NSO fraction recovered by extraction was too low;and second this frac-tion was dried before quantification in order to remove any trace of solvent,which would evaporate low molecular weight compounds associated to this fraction.As such,we cannot directly link the sim-ulations to the experiments for these small oxygenated compounds.Table 4Major reactions found in RD of thermal decomposition of functional models hydroxyl (B)based on the balance of structures recorded during dynamic simulations at 2000,2100and 2200K,50ps.2000K/50ps 2100K/50ps 2200K/50ps Reactant ProductsReactant ProductsReactant Products C 12H 18O 3?C 11H 15O 2+CH 3O C 12H 18O 3?C 12H 18O 3?C 12H 17O 3+HC 11H 15O 2+CH 2O +HC 11H 15O 2+CH 2O +H C 11H 14O 2+CH 2OC 10H 12O 2+CH 3+CH 3O C 10H 12O 2+CH 3+CH 3O C 9H 10O 2+C 3H 6+H 2O C 9H 10O 2+C 3H 6+H 2O C 9H 11O 2+C 3H 6+HO C 8H 10O 2+C 4H 8C 8H 9O 3+C 4H 9C 8H 9O 3+C 4H 9C 8H 9O 3+C 3H 6+CH 3C 8H 9O 3+C 3H 6+CH 3C 8H 8O 2+C 4H 9+OHC 8H 8O 2+C 4H 9+OHC 8H 8O 2+C 4H 9+OHC 8H 6O 2+C 3H 6+CH 4+H 2O0204060801001201401601802000510152025303540Time (ps)510152025303540Time (ps)N u m b e r o f m o l e u c l e s .500100015002000250030003500TotalproductionCO 2H 2O (C 1-C 2) oxygenated compounds H 2(C 1-C 2) hydrocarbons CO Heavy compoundsTemperature (K)T (K)05101520253035404550N u m b e r o f m o l e c u l e s .500100015002000250030003500Temperature (K).H 2CH 4C 2H 4TotalproductionFig.11.Products generated during the thermal decomposition of the macro-model A by NVT simulation under heat up temperature of 60K/ps.1204 E.Salmon et al./Organic Geochemistry 40(2009)1195–1209。