Thermal decomposition of Huadian oil shale: Part 2. The role of bonds in the kerogen skeleton chain during pyrolysis evaluated by ReaxFF MD simulations

Abstract

This study successfully applied reactive force field (ReaxFF) molecular dynamic simulations to study the two-stage transformation of kerogen: initially transforming into pyrolytic bitumen and then into final products. It was found that the carbon–oxygen bond in the kerogen skeleton chain breaks first, while further transformation involves significant carboncarbon bond breakage. Specifically, carbon–carbon bond breakage contributes almost 66.7% to pyrolytic bitumen formation, and nearly 80% to the final product formation. In addition, higher heating rates favor higher bond rupture speed and lead to a decline in heteroatoms within shale oil. In summary, this work provides more atomic insights of oil shale kerogen decomposition.

Headnote

Abstract. This study successfully applied reactive force field (ReaxFF) molecular dynamic simulations to study the two-stage transformation of kerogen: initially transforming into pyrolytic bitumen and then into final products. It was found that the carbon–oxygen bond in the kerogen skeleton chain breaks first, while further transformation involves significant carboncarbon bond breakage. Specifically, carbon–carbon bond breakage contributes almost 66.7% to pyrolytic bitumen formation, and nearly 80% to the final product formation. In addition, higher heating rates favor higher bond rupture speed and lead to a decline in heteroatoms within shale oil. In summary, this work provides more atomic insights of oil shale kerogen decomposition.

Keywords: kerogen, decomposition, bond rupture, molecular dynamics, reactive force field.

(ProQuest: ... denotes formulae omitted.)

1. Introduction

Oil shale is attracting increasing attention due to its vast global reserves and potential as a substitute for crude oil through conversion into shale oil. I n industries, oil shale is usually retorted above ground to produce shale oil, the yield and quality of which are closely related to the kerogen content within oil shale. Ke rogen is an organic matter with a skeleton chain mainly consisting of carbon–oxygen (C–O) and carbon–carbon (C–C) bonds that are insoluble in organic solvents. Based on the hydrogen/carbon (H/C) ratio in descending order, kerogen is classified into three types: I, II and III, with type I kerogen producing the most shale oil. Both the kerogen within the oil shale and the retorting system influence the quality and quantity of the obtained shale oil.

Given this, various parameters of the retorting system, such as the heating rate [1], retorting temperature [2], particle size [3], residence time [4], and pyrolytic atmosphere [5], have been researched to shed light on the kerogen decomposition mechanism and thus improve the quantity and quality of shale oil. The obtained results are beneficial for industrial applications, but provide limited insights into the pyrolytic reaction mechanisms of kerogen. Th erefore, detailed kerogen decomposition mechanisms continue to be a subject of ongoing research.

It is generally accepted that the decomposition of kerogen involves two stages: first, kerogen transforms into pyrolytic bitumen at around 350 °C, and if further heated at a higher temperature, pyrolytic bitumen decomposes into volatiles and coke. The first stage is highly important for further pyrolytic re actions and product formation [6]. Therefore, the chemical structure of the reactants in the first stage, Huadian oil shale kerogen, and the resulting product, pyrolytic bitumen, were studied experimentally [7]. I t was found that compared with kerogen, pyrolytic bitumen showed obvious differences in aliphatic carbons, aromatic carbons, and oxygen (O), nitrogen (N), and sulfur (S) heteroatoms, indicating that there existed complex chemical reactions during kerogen pyrolysis.

To further investigate functional group transformations, in situ Fourier transform infrared (FTIR) analyses were conducted on kerogen pyrolysis, starting from room temperature and heating to 600 °C [8]. It was noted that carbonyl (C=O) groups on the kerogen surface first increase and then decline. To investigate the reaction pathways of functional groups containing O, this work also selected some representative reactants containing hydroxyl (–OH), ether (C–O), and carbonyl (C=O) groups and applied quantum mechanics (QM) calculations to them. The simulation results indicated that ether groups in kerogen aliphatic chain, in the form of R1 –CH2 –O–CH2 –R2 , generated carbon monoxide (CO) through intermediates containing C=O groups rather than H2 O through intermediates containing hydroxyl groups. It can be summarized that the single pathway of representative reactants can be obtained by QM. However, the empirical selection of representative reactants is a trade-off aiming to reduce costs, and thus the statistical characteristic of reactions, such as the proportion of C–O bond contribution to product formation, cannot be obtained.

To simulate chemical reactions without artificial intervention, the reactive force field (ReaxFF) was developed and applied to molecular dynamics (MD) simulations by van Duin et al. [9]. Notably, they calculated the bond dissociation energy of different hydrocarbons using ReaxFF MD [9] and found that the results from ReaxFF MD are much closer to those from the density functional theory (DFT) than those from other MD methods. Subsequently, ReaxFF MD simulation was applied to the oxidation process of hydrocarbons without pre-defined pathways [10], and the results were also consistent with QM calculations. Therefore, ReaxFF MD simulation is a reliable approach to simulate chemical reactions.

So far, ReaxFF MD has been applied to investigate various thermochemical reactions of fuels, such as desulphurization of coal [11], combustion of coal char [12], pyrolysis of Liulin coal [ 13], high-density polyethene [14], and Green River oil shale [15, 16].

In terms of oil shale, Liu et al. [15, 16] investigated the complete decomposition of Longkou and Green River oil shale kerogens into final products, revealing that the ether bridge bond was the first to break, followed by the rupture of the paraffin chain. However, the extent to which the ether bridge (C–O bond) and paraffin chain breakage (C–C bond) contribute to kerogen transformation remains uncertain. Furthermore, bond rupture and product evolution might be influenced by the heating rate, but the latter was only considered for validating the simulation in this study [16], and the influence of heating rates on bond breakage still requires further investigation.

Gi ven this, the current study used the large-scale atomic/molecular massively parallel simulator (LAMMPS) to conduct a series of MD simulations. The main objectives were to determine the impact of the breakdown of skeleton chain bonds on kerogen transformation and the effect of heating rates on bond breakage.

2. Methods

2.1. The construction of the kerogen macromolecular model and simulation cell

As the first step to conduct simulation, the construction of the kerogen molecular model was based on a series of characteristic experiments, including ultimate analysis and nuclear magnetic resonance (NMR) analysis. The constructed kerogen model's theoretical NMR spectra were obtained by DFT and closely resembled the experimental NMR spectra [17]. The 2D plot of the kerogen molecule is shown in Figure 1. Evidently, in the aliphatic chain, C–O bonds were much rarer than C–C bonds, while both played significant roles in kerogen transformation.

Subsequently, the kerogen molecule was packed into a simulation cell. Specifically, five kerogen molecules were packed into a cell, the size of which was based on the chosen density. The density of the cell was set to 1.0 g/cm3 , which fell within the density range (0.93–1.05 g/cm3 ) used in references [18–20]. The obtained cell was further optimized before exporting data to LAMMPS for simulation. The detailed process can be found in reference [17].

2.2. Simulation details in LAMMPS

The ReaxFF used in this study was developed by Liu et al. [21], and the overall simulation process was quite similar to the one used in reference [17]. Specifically, the general process of simulations consisted of four stages. In addition, two ultimate temperatures were set for the third and fourth stages to simulate the breakdown of skeleton chain bonds contributing to the two stage decomposition. The lower temperature, 1900 K, was set to simulate the transformation of kerogen into bitumen, while the higher temperature, 2500 K, was set to simulate the transformation of kerogen into final products [17]. The 1900 and 2500 K in the simulation corresponded to 364 and 506 °C in the experiments, respectively, as established by Equation (1) provided by Xu et al. [17, 22]:

. . . ( 1)

where Texp is the temperature of the thermogravimetric (TG) experiment, Tsim is the temperature of the simulation, Ea is the activation energy of oil shale (kJ/mol), and R isthe gas constant (kJ·mol–1·K–1). In addition, t sim and t exp are the duration times of the simulation works and the TG experiment. Specifically, t sim / t exp equals 10–12 (10–3 s / 109 s), and Ea is 240 kJ/mol, considering that nearly half of the kerogen is decomposed into final products.

Furthermore, this work applied different heating rates ranging from 6 to 1200 K/ps to the third heating stage. Specifically, three heating rates (30, 60, and 1200 K/ps) were applied to the simulation of kerogen transformation into pyrolytic bitumen, and six heating rates (6, 8, 12, 24, 48, and 1200 K/ps) were applied to the simulation of kerogen decomposition into final products. Notably, the total duration of the third and fourth stages was set to 450 ps at 6 K/ps, because the third heating stage had already taken 400 ps, while the total duration at other heating rates was 350 ps. Based on the heating rates and ultimate temperature applied to the third heating stage, the simulations were labeled as Temperature@HeatingRate. For example, the label 1900K@1200 indicated that the ultimate temperature and heating rate applied to the third stage were 1900 K and 1200 K/ps, respectively. In addition, another simulation labeled 3100K@12 was conducted to validate the char mass loss behavior during the final product formation.

2.3. Data post-processing

The obtained data were first processed using Python packages released in reference [23] and subsequently by a self-developed program. The intermediate products were classified into five categories according to their carbon number. In addition, the bond number in the simulation system was monitored using a self-developed program based on OVITO software [24], which applied a cutoff distance of 1.5–1.6 Å and thus could not distinguish between single and double bonds. However, C–C and C–O bonds are quite common in oil shale kerogen and are more likely to break than their corresponding double bonds [30]. Therefore, the CO and CC bond variations obtained by the self developed program and plotted in the figures could represent C–O and C–C bond variations.

2.4. Validation

Although ReaxFF molecular dynamic simulation has been widely used to investigate combustion and pyrolysis reactions of numerous materials, this study still conducted validation on the simulation process. First, the bond breakage order of the 1900K@1200 simulation was summarized and compared with the energy barriers of the corresponding bond rupture as calculated by the DFT. Second, the C–O bond transformation was further investigated and compared with findings from FTIR experiments. Third, during the formation of final products, the char mass profiles of the simulation system (2500K@12 and 3100K@12) were compared with those from TG experiments. These validation steps confirmed that the simulations were reliable.

3. Results and discussion

3.1. Validation of simulations

3.1.1. ReaxFF MD and DFT: aligning bond rupture sequences and energy barriers

Du ring bitumen formation, C–O and C–C bonds in aliphatic chains break. Bonds containing O and N were paid greater attention and labeled A–G, as shown in Figure 1. The breakage order during bitumen formation (1900K@1200) is summarized in Figure 2a. A–G bonds are found in every kerogen, and the rupture of these bonds can occur more than once because the simulation cell consists of five kerogens. The corresponding time steps of every individual bond rupture are plotted in Figure 2a as scatter points, and the continuous line shows the mean time step. From the mean time step, the bonds on the right side of Figure 2a break later than those on the left side, except for the H' bond. Specifically, ether bonds (B and C, C–O–C), ester bonds connecting methine (A, D and E, R–CH(–R1)–C(=O)–O–R2), hydroxyl bonds (H', R–O), and amino bonds (F, R–N) break in order. These findings are consistent with those in references [5, 6], which state that bonds in connection with ester, ether, and hydroxyl groups break first.

In addition, bond rupture order is closely related to energy barriers. Therefore, Figure 2b summarizes the energy barriers for the rupture of ether bonds in aliphatic chains (B and C), hydroxyl bond (H'), amino bond (F), and amide bond (G) as calculated by the DFT. It can be observed that the energy barrier values for B and C, H', F, and G bonds increase in this order, which also coincides with the bond breakage order found by ReaxFF MD simulations.

3.1.2. ReaxFF MD and FTIR: consistency of C–O bond transformation

The broken bonds in ether and ester are consistently found to be C–O bonds. Consequently, the resulting radical containing C–[O] is reactive and tends to transform into another radical, which can be detected by FTIR experiments. Therefore, attention was given to the resulting fragments from the rupture of C–O bonds during the simulation process. As Figure 2c shows, the R–C–[O] radical generated by the breakage of C–O bonds within ether, ester, and cyclic ether gradually transforms into a carbonyl group and then remains stable until further pyrolysis reactions occur. This transformation of C–O bonds into C=O bonds was also observed in the formation of formaldehyde during coal pyrolysis [15], the formation of furan during the pyrolysis of ɑ-cyclodextrin (a small molecule surrogate of cellulose) [25], and the formation of intermediates during the pyrolysis of L-rhamnose monohydrate (a monomer of polysaccharides in seaweed) [26]. As a result, the carbonyl content within the system increased, which is consistent with FTIR findings showing that carbonyl content gradually increased during the initial stages of pyrolysis [8].

3.1.3. ReaxFF MD and TG: partial consistency of char mass profiles

During the simulation, the char mass fraction was a key focus. Among the six simulations of kerogen transformation into final products, the one with a heating rate of 12 K/ps, 2500K@12, was compared with a TG experiment, which had a heating rate of 10 K/min and a final temperature of 1273.15 K [27]. The final temperature in the TG experiment was much higher than 779.15 K, corresponding to the ultimate temperature of 2500 K in the ReaxFF MD simulation. Therefore, a higher ultimate temperature of 3100 K was also applied in the ReaxFF MD simulation, labeled as 3100K@12. In addition, the used macromolecule of kerogen in these simulations was based on the kerogen from the Dachengzi diggings located in Huadian, Jilin Province, which is also the source of the kerogen used in the TG experiment.

It was found that the char mass fractions of both 2500K@12 and 3100K@12 showed only slight differences compared to the TG experiment until reaching 30.52%. Furthermore, when the char mass fraction in 3100K@12 reached 30.52%, the system temperature (Tsim) was 2632.8 K; in contrast, Texp was 772.7 K in the TG experiment, when the char mass fraction reached 30.52%. Assuming an activation energy of 255 kJ/mol, which averaged the activation energy of kerogen transformation into final products as reported in the references [17, 28], Tsim of 2632.8 K corresponded to Texp of 781.1 K, which was also very close to the observed 772.7 K.

As the simulation proceeded further, the discrepancy between the mass fraction in the TG experiment and that in the simulation became increasingly apparent. This discrepancy has also been noted in reference [29] and can be attributed to the accumulation of small-molecule products. These intermediate products with low molecular weight might favor radical propagation, leading to more intense secondary reactions, including further decomposition. In the TG experiments, these small-molecule intermediate products were released from char as volatiles, preventing them from contributing to further decomposition. As a result, the discrepancy between the TG experiments and the simulation increased.

In summary, the results from the conducted simulations show partial consistency with findings from other experimental studies and simulations. Thus, the decomposition of kerogen into final products can be divided into two stages: the initial stage, consistent with TG experiments, which can provide additional molecular insights into kerogen decomposition, and the later stage, which explores the potential mechanisms triggered by the accumulation of small-molecule radicals.

3.2. The role of C–C and C–O bonds in the kerogen skeleton chain during pyrolytic bitumen formation

An ultimate temperature of 1900 K was applied to the simulation of kerogen transformation into pyrolytic bitumen over 350 ps. During the simulation, the char mass fraction was calculated and plotted in Figure 3a. As shown, the char mass fraction decreased from 100 to 82.08% (1900K@30), 72.08% (1900K@60), and 78.51% (1900K@1200). In addition, the char fraction at all heating rates decreased slightly by the end of the simulation, indicating the transformation of kerogen into bitumen. Furthermore, the slight discrepancy in char mass fraction between 1900K@60 and 1900K@1200 proved that a heating rate of 1200 K/ps is feasible for this simulation.

Figure 3b shows that the number of C–C bonds decreased by 42 (1900K@30) and 49 (1900K@60 and 1900K@1200) with the char mass fraction declining. At the same time, products of shorter chains (C<20) were formed through random scission reactions or chain-end scission reactions, where C–C bonds were the dominant bonds broken due to their lower dissociation energy compared to C=C or C≡C bonds [30].

However, as shown in Figure 3a, the increase in the mass fraction of shorter-chain products only partially compensated for the decrease in the char mass fraction. Specifically, the mass fraction of shorter-chain products increased from 0 to 9.81% (1900K@30), 12.58% (1900K@60), and 14.14% (1900K@1200), which accounted for just 50% of the char mass loss. Therefore, it can be inferred that the char mass loss is also linked with the breakage of C–O bonds.

As Figure 3b shows, the number of C–O bonds declined by 23 (1900K@30), 32 (1900K@60), and 22 (1900K@1200) as the char mass decreased. In terms of the total number of bond changes, the C–C bond decreases twofold compared to the C–O bond. Figure 3c illustrates that the combined breakage of C–O and C–C bonds was closely matched but slightly bigger than the total number of hydrocarbons produced. This suggests that the breakage of C–O and C–C bonds primarily resulted in the formation of pyrolytic bitumen, with a minor contribution to the formation of organic gas and ring rupture. More specifically, C–C bond breakage contributed nearly 67% to the char mass loss during pyrolytic bitumen formation, while C–O bond breakage accounted for 33%. In addition, the more gradual char mass loss observed in the later stages of the simulation indicated that kerogen had nearly completed its transformation into bitumen.

Kerogen decomposition is a radical reaction, and radical reactions usually start with homolytic fission, which is affected by thermal radiation and strongly related to the temperature applied in the system. The resulting radicals participate in subsequent radical propagation reactions, which evolve over time. In general, these reactions are influenced not only by temperature but also by time. Therefore, the bond variation speed was calculated and plotted in Figure 3d with respect to the heating speed. Specifically, Figure 3d demonstrates the bond variation speed of the system at 1900 K. It can be seen that the rupture speed of C–O and C–S bonds shows an escalating trend as the heating rate increases, whereas that of C–C and C–N bonds does not show a clear trend. Previous findings indicate that the energy barrier for the C–N bond is higher than that for the C–O bond. If the applied temperature provides sufficient energy to overcome this energy barrier, higher heating rates could favor higher bond rupture speed.

3.3. The role of C–C and C–O bond breakage during the formation of final pyrolytic products

An ultimate temperature of 2500 K was used to simulate the transformation of kerogen into final pyrolytic products. The resulting char fraction was calculated and plotted in Figure 4a. The char mass declined with the rupture of C–C and C–O bonds, with the total bond rupture shown in Figure 4b. Nearly half of the C–O bonds were in the form of carboxyl located at the end of the kerogen chain, while the other half were located in the middle of the kerogen skeleton chain, as shown in Figure 1. Evidently, the rupture of the C–O bond in carboxyl can result in H2 O formation. As shown in Figure 5a and b, the number of H2 O molecules was almost consistent with a half of the C–O bond decrease. However, the maximum accumulation of H2 O was 45 at 2500K@48, meaning that H2 O formation related to C–O bond rupture only accounted for no more than a 3% decrease in the char mass loss. Moreover, Figure 4b shows that the total number of broken C–C bonds was almost 300% higher than that of broken C–O bonds. The combined number of broken C–O and C–C bonds (as shown in Figure 4b) was also 20% bigger than the total increase in the number of hydrocarbons (Figure 4d), regardless of whether the char mass fraction was at 30 or 5%. This indicates that C–O and C–C bond breakage plays a dominant role in the formation of final pyrolytic hydrocarbon products, with C–C bond breakage contributing nearly 80%.

Interestingly, C–C bond breakage showed a declining trend in the later stages of the simulation at 2500K@1200, 2500K@48, 2500K@24, and 2500K@12, as shown in Figure 4b. This trend was also consistent with the char mass fraction re-increases of 2500K@1200 and 2500K@48, indicating the start of coking reactions in these systems. Based on these findings, the whole process was divided into two stages: the first stage involved a decline in char from 100 to 30%, and the second stage involved further reactions, including decomposition and coking reactions. In the first stage of 2500K@12, the char mass profile was nearly consistent with the TG experimental results; therefore, the overall chemical reactions in this stage were more likely representative of real experimental reactions. Even when the char mass profile showed a difference compared to the TG experiment, the chemical reactions found in ReaxFF MD provide valuable insights into the decomposition mechanisms.

3.4. Enhanced bond rupture with increased heating rates

Figure 4a and b shows that char mass, C–C bond breakage, and C–O bond breakage vary with different heating rates. To illustrate how heating rates influence chemical reactions, Figure 4c summarizes the bond rupture speed variation with heating rates in the first decomposition stage. During this stage, the rupture speed of C–C bonds increased with heating rates, whether measured when the pyrolytic system first reached 2500 K or when the char mass fraction equaled 30 wt%. This trend was also observed for C–O and C–S bond rupture speeds. Apparently, higher heating rates favored increased bond rupture speeds, provided the applied temperature supplied sufficient energy for bond rupture. With higher C–C bond rupture speeds, there were more concentrated radicals of medium molecular weight (carbon chain 5≤C<10). These radicals mainly originated from random C–C bond rupture rather than chain-end scission, and could further promote radical reactions. Consequently, Figure 4c shows that as char mass decreased further to 30 wt%, the C–N bond broke more quickly.

The radicals obtained in the first stage favored further decomposition in the second stage, where the char mass loss increased. Figure 4d shows that as the char mass decreased from 30 to 5%, the number of total resulting radicals, especially the number of short chain radicals whose total carbon number was less than 10, would occur at higher heating rates. These radicals underwent further radical reactions. As the reactions progressed towards completion, the char mass loss fraction for 2500K@1200 and 2500K@48 increased again, as evident in Figure 4a. Figure 4e plots the changes in the number of char and medium radicals from 5 wt% to the end of the simulation. It can be seen that, with the char number increasing, the number of radicals of medium molecular weight (C<20) decreased for 2500K@1200 and 2500K@48. Thus, these obtained radicals of medium molecular weight (C<20) were more likely to rearrange with each other and result in coke. To avoid coking reactions, the concentrated radicals of medium molecular weight (C<20) resulting from pyrolysis at higher heating rates should be removed promptly.

As previous results indicate, higher heating rates favored bond rupture and resulted in an increased formation of short carbon chain radicals during the transformation of kerogen. These short carbon chain radicals played a significant role in the propagation of radical reactions, which further influenced the distribution of heteroatoms (O, S, and N). Figure 5a–d summarizes the number of these products at the initial time of coking and at the end of the simulation. The initial time of coking was defined as the moment when C–C bonds started rebounding. Specifically, the corresponding indexes at this time are plotted in the blue area and the corresponding tick label contains the suffix "-c". For simulations at 2500K@6 and 2500K@8, no rebound of C–C bonds was observed, and thus the initial time of coking was supposed to be the same as the end time of the simulation.

Inorganic heteroatom O was mainly found in CO2 and H2 O, with the corresponding numbers plotted in Figure 5a. Furthermore, the total number of O atoms in inorganic gas was calculated and shown in Figure 5a, where the value of "reference O" equals the number of CO2 molecules. As Figure 5a shows, the total number of O atoms in inorganic gas increased with escalating heating rates, which indicates that higher heating rates favor the accumulation of gaseous inorganic O.

In practice, the organic products of small molecules (C<3) are also gaseous and can be easily separated from the major organic contents, while the remaining liquid organic products are classified as pyrolytic oil. To investigate the effect of heating rate on the O content of pyrolytic oil, the O content of gaseous products, including both inorganic and organic gases, was plotted in Figure 5b. It can be observed that the O content of gaseous products shows an increasing trend with rising heating rates. In other words, higher heating rates reduce the O content in pyrolytic oil.

Similarly, heteroatoms S and N in inorganic gas were found in H2 S and NH3 , respectively, as shown in Figure 5c and d. H2 S exhibited an increasing trend with escalating heating rates, peaking at 12 K/ps. In contrast, NH3 did not show an obvious trend. This behavior can be ascribed to the high reactivity of the resulting NH3 , which would undergo further reactions.

3.5. Formation pathways of small molecules during the transformation of kerogen into final products

As said above, some small molecular weight species were also released during the transformation process. Therefore, the formation pathways of these small molecules were explored in detail.

3.5.1. H2 O, H2 S and NH3 formation pathways H2 O formation is strongly linked to the hydroxyl radical resulting from the dissociation of hydroxyl groups. As Figure 1 shows, the kerogen macro molecule contains only one H' bond. The single hydroxyl radical resulting from H' bond rupture forms just one H2 O molecule by capturing hydrogen. This pathway for H2 O formation is also shown in references [15, 16]. How ever, Figure 5b illustrates that the final accumulative number of H2 O is far more than one per kerogen molecule. This indicates that there should be other pathways for H2 O formation. An examination of all the serial tests, shown in Figure 6, revealed alternative formation pathways for H2 O formation. Apart from G bond breakage, the hydroxyl radical can also result from carboxyl rupture. Then, the resulting hydroxyl radical seizes the hydrogen atom from the paraffin chain or H radical to form H2 O.

The formation pathways for NH3 and H2 S are also presented in Figure 6. The amino groups within oil shale kerogen exist in two forms: as an amide group and as an amino group at the end of the aliphatic chain. The corresponding bonds connecting these two groups are presented in Figure 1 and labeled as G bond and F bond, respectively. According to reference [31], each amino group in kerogen can decompose into an amino radical. This amino radical may further catch hydrogen from the aliphatic chain, forming NH3 . In this study, the amino group connecting the kerogen chain through the F bond consistently followed this pathway, resulting in the formation of one NH3 per kerogen. However, the amide group, including the Hb bond, seldom broke to form an amino radical and contribute to NH3 formation. Furthermore, the sulfhydryl at the end of the paraffin chain was found to break frequently.

The resulting sulfhydryl radical could convert into H2 S either by combining with an H radical or by seizing a hydrogen atom from other intermediates, as also reported in reference [11]. The resulting NH3 was very reactive and could react with other radicals. Therefore, it could be seen that the number of NH3 molecules declined as the system proceeded with further reactions. Reference [32] indicates that NH3 can react with species containing carboxyl to transform into species containing ammonium, which may further produce species containing amide and finally become nitrile species. These nitrile species are highly likely to remain in shale oil, leading to an increase in N content. However, the previous analysis indicates that higher heating rates can favor the hydroxyl dissociation from carboxyl and thus prevent carboxyl from reacting with NH3 to form nitrile species. This implies that heating rates may reduce both N and O content in shale oil.

3.5.2. Formation pathways for other small species

The previous analysis indicates that the formation of NH3 , H2 O, and H2 S is strongly related to their corresponding functional groups and results from radical reactions. In summary, these radical reactions lead to the rupture of aliphatic chains in polymers during both pyrolysis and degradation [14, 33]. Notably, aliphatic chains are common in kerogen, which means that the decomposition of kerogen involves radical reactions.

Generally, radical reactions include random scission reactions and chain end scission reactions. The former usually result in a significant decrease in molecular weight of the original products, whereas the latter release volatile products without a substantial decrease in the molecular weight of the original products. Figure 6 summarizes the formation of volatile products of small molecular weight through chain-end scission reactions during kerogen pyrolysis. Specifically, the C–C bond connecting carbonyl or close to carbonyl breaks to form CO2 , ethenone (C2 H2 O), and radicals containing carbonyl. The release of these molecules can reduce the carbonyl content within the oil shale char. This pathway could explain the FTIR experimental findings, which show that the carbonyl percentage within kerogen declines as the temperature exceeds 350 °C [8].

In addition, as Figure 6 shows, the chain-end scission reactions usually generate one or two radicals. These resulting radicals can further adsorb the H atom from the aliphatic chain, resulting in reactive radicals, which then contribute to further radical reactions. This phenomenon, known as radical propagation, is inevitably related to C–C bond rupture. Most importantly, Figure 4c shows that C–C bond rupture is more favorable at higher heating rates, leading to an increase in volatile products. Consequently, the radical propagation initialized by these volatile products is also enhanced by higher heating rates, suggesting that higher heating rates accelerate the decomposition of kerogen. If the accumulation of volatile products is sufficient, these radicals may combine directly with each other and result in bigger molecule products, such as char, which accounts for why the char mass fraction rebounds earlier at higher heating rates, as shown in Figure 4a.

4. Conclusions

To explore the role of bonds within the skeleton chain on the thermal decomposition of kerogen, this study used ReaxFF MD to simulate the pyrolytic reactions of Huadian oil shale kerogen. The results are as follows.

1. Extensive efforts indicate that the simulations in this paper are reliable. Firstly, the bond rupture order is consistent with the corresponding energy barrier calculated by the DFT. Secondly, it was found that the R–C–[O] radical, resulting from C–O breakage in ether, cyclic ether, and ester, gradually transforms into R–C=O, which results in an increase in C=O, consistent with the experimental FTIR findings. Thirdly, the char mass loss observed in simulations at specific heating rates partially aligns with the TG experiment results.

2. The skeleton of kerogen mainly involves C–O and C–C bonds, whose breakage plays a significant role in bo th pyrolytic bitumen formation and final pyrolytic product formation. C–C bond breakage contributes 66.7% to pyrolytic bitumen formation, while C–O bond breakage accounts for the remainder. In addition, some C–O bonds, such as those in ester close to indole and connecting methine, ether bonds (C–O–C), ester bonds connecting methine (R–CH(–R1)–C(=O)–O–R2), and ester bonds connecting methylene (R–CH2 –C(=O)–O–R1), break successively. Furthermore, the breakage of C–C bonds within paraffin chains plays a more significant role in final product formation, contributing nearly 80%. C–O bond breakage is associated with the formation of small molecule species containing O, such as H2 O.

3. In addition, increasing heating rates favor random C–C bond rupture speed, resulting in the accumulation of medium molecular weight radicals. These radicals are beneficial for radical propagation, potentially involving chain-end scission reactions that release small molecule species containing heteroatoms, thus reducing the content of heteroatoms within pyrolytic oil. This study quantitatively determines the contribution of bond rupture in the kerogen skeleton to product formation and sheds light on the effect of heating rates on bond rupture and corresponding product formation at the atomic scale.

Acknowledgments

This study was supported by th e N ational Natural Science Foundation of China (grants No. 52 200155, No. 51 876122, and No. 51766004). In addition, the reviewers and editors provided many valuable comments to improve this study. The publication costs of this article were partially covered by the Estonian Academy of Sciences.

Sidebar

Received 20 April 2024, accepted 18 October 2024, available online 31 October 2024

Footnote

* Corresponding authors, [email protected], [email protected]

References

References

1. Wang, S., Liu, J., Jiang, X., Han, X., Tong, J. Effect of heating rate on products yield and characteristics of non-condensable gases and shale oil obtained by retorting Dachengzi oil shale. Oil Shale, 2013, 30(1), 27–47. https://doi. org/10.3176/oil.2013.1.04

2. Wang, S., Jiang, X., Han, X., Tong, J. Effect of retorting temperature on product yield and characteristics of non-condensable gases and shale oil obtained by retorting Huadian oil shales. Fuel Process. Technol., 2014, 121, 9–15. https:// doi.org/10.1016/j.fuproc.2014.01.005

3. Ahmad, N., Williams, P. T. Influence of particle grain size on the yield and composition of products from the pyrolysis of oil shales. J. Anal. Appl. Pyrolysis, 1998, 46(1), 31–49. https://doi.org/10.1016/S0165-2370(98)00069-2

4. Salhi, N., Bennouna, C., Bitar, H., Sergent, M., Luu, R. P. T. An experimental design to optimize pyrolysis conditions of Timahdit (Morocco) oil shale. Fuel, 1996, 75(5), 633–640. https://doi.org/10.1016/0016-2361(95)00258-8

5. Hershkowitz, F., Olmstead, W. N., Rhodes, R. P., Rose, K. D. Molecular mechanism of oil shale pyrolysis in nitrogen and hydrogen atmospheres. In: Geochemistry and Chemistry of Oil Shales (Miknis, F. P., McKay, J. F., eds). American Chemical Society, Washington, 1983, 301–316. https://doi. org/10.1021/bk-1983-0230.ch015

6. Wen, C. S., Kobylinski, T. P. Low-temperature oil shale conversion. Fuel, 1983, 62(11), 1269–1273. https://doi.org/10.1016/S0016-2361(83)80008-8

7. Li, Q., Han, X., Liu, Q., Jiang, X. Thermal decomposition of Huadian oil shale: Part 1. Critical organic intermediates. Fuel, 2014, 121, 109–116. http://dx.doi. org/10.1016/j.fuel.2013.12.046

8. Chen, B., Han, X., Jiang, X. In situ FTIR analysis of the evolution of functional groups of oil shale during pyrolysis. Energy Fuels, 2016, 30(7), 5611–5616. https://doi.org/10.1021/acs.energyfuels.6b00885

9. Duin, A. C. T., van, Dasgupta, S., Lorant, F., Goddard, W. A. ReaxFF: a reactive force field for hydrocarbons. J. Phys. Chem. A, 2001, 105(41), 93969409. https://doi.org/10.1021/jp004368u

10. Chenoweth, K., Duin, A. C. T., van, Goddard, W. A. ReaxFF reactive force field for molecular dynamics simulations of hydrocarbon oxidation. J. Phys. Chem. A, 2008, 112(5), 1040–1053. https://doi.org/10.1021/jp709896w

11. 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. https://doi.org/10.1016/j.fuel.2014.12.074

12. Castro-Marcano, F., Kamat, A. M., Russo, M. F., Jr., Duin, A. C. T., van, 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. https://doi.org/10.1016/j.combustflame.2011.10.022

13. Zheng, M., Li, X., Liu, J., Wang, Z., Gong, X., Guo, L., Song, W. Pyrolysis of Liulin coal simulated by GPU-based ReaxFF MD with cheminformatics analysis. Energy Fuels, 2013, 28(1), 522–534. https://doi.org/10.1021/ef402140n

14. Liu, X., Li, X., Liu, J., Wang, Z., Kong, B., Gong, X., Yang, X., Lin, W., Guo, L. Study of high density polyethylene (HDPE) pyrolysis with reactive molecular dynamics. Polym. Degrad. Stab., 2014, 104, 62–70. https://doi.org/10.1016/j. polymdegradstab.2014.03.022

15. 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), 2987–2997. https://doi.org/10.1021/acs. energyfuels.5b00084

16. Qian, Y., Zhan, J.-H., Lai, D., Li, M., Liu, X., Xu, G. Primary understanding of non-isothermal pyrolysis behavior for oil shale kerogen using reactive molecular dynamics simulation. Int. J. Hydrog. Energy, 2016, 41(28), 12093–12100. https://doi.org/10.1016/j.ijhydene.2016.05.106

17. Mu, M., Han, X., Wang, S., Jiang, X. Interactions of oil shale and hydrogen rich wastes during co-pyrolysis: co-pyrolysis of oil shale and waste tire. Energy Fuels, 2023, 37(6), 4222–4232. https://doi.org/10.1021/acs.energyfuels.2c02954

18. Pan, S. Molecular Simulation on Oil Shale Kerogen Molecular Structure during Pyrolysis Reaction. Master's thesis. Northeast Electric Power University, China, 2020.

19. Wang, Q., Cheng, F., Pan, S. Chemical bond concentration and energy density of oil shale kerogen. J. Fuel Chem. Technol., 2017, 45(10), 1209–1218.

20. Ru, X., Cheng, Z., Song, L., Wang, H., Li, J. Experimental and computational studies on the average molecular structure of Chinese Huadian oil shale kerogen. J. Mol. Struct., 2012, 1030, 10–18. https://doi.org/10.1016/j.molstruc.2012.07.027

21. Liu, L., Liu, Y., Zybin, S. V., Sun, H., Goddard, W. A. ReaxFF-lg: correction of the ReaxFF reactive force field for London dispersion, with applications to the equations of state for energetic materials. J. Phys. Chem. A, 2011, 115(40), 11016–11022. https://doi.org/10.1021/jp201599t

22. Xu, H., Yu, H., Fan, J., Xia, J., Liu, H., Wu, H. Formation mechanism and structural characteristic of pore-networks in shale kerogen during in-situ conversion process. Energy, 2022, 242, 122992. https://doi.org/10.1016/j. energy.2021.122992

23. Döntgen, M., Przybylski-Freund, M.-D., Kröger, L. C., Kopp, W. A., Ismail, A. E., Leonhard, K. Automated discovery of reaction pathways, rate constants, and transition states using reactive molecular dynamics simulations. J. Chem. Theory Comput., 2015, 11(6), 2517–2524. https://doi.org/10.1021/acs.jctc.5b00201

24. Stukowski, A. Visualization and analysis of atomistic simulation data with OVITO – the Open Visualization Tool. Model. Simul. Mater. Sci. Eng., 2010, 18(1), 015012. https://doi.org/10.1088/0965-0393/18/1/015012

25. Mettler, M. S., Mushrif, S. H., Paulsen, A. D., Javadekar, A. D., Vlachos, D. G., Dauenhauer, P. J. Revealing pyrolysis chemistry for biofuels production: conversion of cellulose to furans and small oxygenates. Energy Environ. Sci., 2012, 5(1), 5414–5424. https://doi.org/10.1039/c1ee02743c

26. Yuan, C., Jiang, D., Wang, S., Barati, B., Gong, X., Cao, B., Zhang, R., Zhang, C., Odey, E. A. Study on catalytic pyrolysis mechanism of seaweed polysaccharide monomer. Combust. Flame, 2020, 218, 1–11. https://doi. org/10.1016/j.combustflame.2020.04.022

27. You, Y., Wang, X., Han, X., Jiang, X. Kerogen pyrolysis model based on its chemical structure for predicting product evolution. Fuel, 2019, 246, 149–159. https://doi.org/10.1016/j.fuel.2019.02.075

28. Bai, F., Guo, W., Lü, X., Liu, Y., Guo, M., Li, Q., Sun, Y. Kinetic study on the pyrolysis behavior of Huadian oil shale via non-isothermal thermogravimetric data. Fuel, 2015, 146, 111–118. https://doi.org/10.1016/j.fuel.2014.12.073

29. Feng, W., Zheng, M., Jin, L., Bai, J., Kong, L., Li, H., Bai, Z., Li, W. Co pyrolysis behaviors of coal and polyethylene by combining in-situ Py-TOF MS and reactive molecular dynamics. Fuel, 2023, 331, 125802. https://doi. org/10.1016/j.fuel.2022.125802

30. Hou, S., Xu, Y. Organic Chemistry. Higher Education Press, Beijing, 2015.

31. Han, X., Chen, B., Li, Q., Tong, J., Jiang, X. Organic nitrogen conversion during the thermal decomposition of Huadian oil shale of China. Oil Shale, 2017, 34(2), 97–109. https://doi.org/10.3176/oil.2017.2.01

32. Evans, E. J., Batts, B. D., Cant, N. W., Smith, J. W. The origin of nitriles in shale oil. Org. Geochem., 1985, 8(5), 367374. https://doi.org/10.1016/01466380(85)90015-4

33. Singh, B., Sharma, N. Mechanistic implications of plastic degradation. Polym. Degrad. Stab., 2008, 93(3), 561–584. https://doi.org/10.1016/j. polymdegradstab.2007.11.008

34. Peng, K., Wu, Q., Tipeng, W., Zongqiang, F., Liwei, J., Zhongfu, T. Generation mechanism of NOx and N2 O precursors (NH3 and HCN) from aspartic acid pyrolysis: a DFT study. Int. J. Agric. Biol. Eng., 2016, 9(5), 166–176. https:/doi. org/10.3965/j.ijabe.20160905.2559