Experimental and Molecular Simulation Studies of Huadian Oil Shale Kerogen.

Microscopic details on the intrinsic chemical reactivity of Huadian oil shale kerogen associated with electron properties of kerogen were investigated by the combination of experimental analyses and molecular simulations. Multimolecular structure models of kerogen with different densities were constructed for examining the accuracy of the proposed kerogen model. Results revealed that the simulated density of the kerogen model is in good agreement with the experimental value. Evaluation of the kerogen model revealed that the energy optimization process is mainly related to the change in the bond angle caused by atom displacement. According to the results from the Hirshfeld analysis of atomic charges, the S atoms in thiophene and S=O structures exhibit positive charges. By contrast, the concentration of electrons on the S atom led to the electronegativity of the sulfhydryl group. To investigate the distribution characteristics of electrons in kerogen, the molecular electrostatic potential (MEP) of a complete kerogen molecule was calculated. Notably, this is the first report of an MEP diagram of the kerogen model that can provide valuable information on the determination of electrophilic or nucleophilic reaction sites for kerogen.

Introduction

Oil shale is a sedimentary rock that can be used as a nonrenewable solid fossil fuel.[1,2] Oil shales are mainly composed of inorganic minerals and organic matter.[3] As the main component of organic matter compactly mixed in minerals (such as montmorillonite, calcite, etc.),[4−6] kerogen primarily comprises five elements of carbon, oxygen, nitrogen, hydrogen, and sulfur.[7,8]

In terms of elemental composition, the hydrogen-to-carbon (H/C) and oxygen-to-carbon (O/C) ratios are the most important physicochemical indexes. Kerogens can be categorized into three types: Type I (lacustrine); Type II (marine) or IIS (S-containing); and Type III (humic). In Type I, the H/C atomic ratio is frequently greater than 1.5, and the O/C atomic ratio is between 0.03 and 0.1. Aliphatic compounds are dominant in Type I kerogen, such as the Green River oil shale in America.[9] In Type II, H/C and O/C atomic ratios are typically ∼1.3 and 0.15, respectively. Compared to Type I kerogen, Type II kerogen comprises highly aliphatic compounds, containing more aromatic hydrocarbon compounds, and thiophene-type compounds, as well as a large amount of cyclic hydrocarbons, such as the Ghareb Formation oil shale in Jordan.[10] Jordanian marinite of the Upper Cretaceous Ghareb Formation contains type-IIS kerogen and is rich in a complex distribution of organic sulfur compounds, among which the alkylthiophenes can be used as molecular biomarkers to indicate paleoenvironmental changes.[11] In this case, type II kerogens are further classified as belonging to Type IIS (S-containing) kerogens.[12] Type II kerogen has intermediate characteristics between Type I and Type III and is mainly derived from organic matter deposited in a marine setting, such as the New Albany Shale within the Illinois Basin.[13] Type III kerogen is derived from terrestrial plant materials and has a low H/C ratio. In general, Type III kerogens are more gas-prone and coal-like and therefore they do not generate commercial quantities of liquid hydrocarbons upon pyrolysis in most cases.

From the perspective of the chemical structure model, the molecular models of oil shale or kerogen that have been developed over the past 30 years are productively investigated via the development and comprehensive utilization of various experimental technologies. Research on the molecular structure of kerogen started late compared with that on the chemical structure model of coal. Behar and Vandenbroucke[14] have characterized the functional groups and thermal decomposition properties of three classical types of kerogens with different thermal maturities (i.e., Green River kerogen for Type I, Toarcian kerogen for Type II, and Douala kerogen for Type III) using elemental analysis, electron microscopy, 13C nuclear magnetic resonance (NMR) spectroscopy, and thermogravimetry analysis. In their study, two-dimensional (2D) molecular structure models of kerogens were constructed on the basis of the statistical analytic method. Behar and Vandenbroucke have comparatively proposed a relatively complete modeling process of kerogen chemical models and provided quantitative and qualitative analysis results for structural units in terms of aromatic clusters, alkyl chains, and heteroatom-containing compounds, such as porphyrins. As presented in their kerogen models, the molecular components are mainly connected by covalent bonds, such as C–C and C–heteroatoms (O, N, and S). Meanwhile, as weak interactions, hydrogen bonds, van der Waals forces, and intramolecular π bonds also play a role in connecting the molecular structure. Siskin et al.[15] have employed solid-state 13C nuclear magnetic resonance (13C NMR) and 29Si NMR spectroscopy for the chemical identification of organic matter in Rundle Ramsay Crossing oil shale and Green River oil shale. From analytical data of the carbon skeleton structure of kerogen, the structural characteristics of molecular fragments and shale oil products were quantitatively obtained. Representative structural models of Rundle Ramsay Crossing and Green River oil shales were constructed with the empirical formulas of C100H160N2.2S0.68O9.22 and C100H155N3.10S0.62O2.17, respectively. The chemical structure parameters of molecular models, such as element composition, aromatic carbon ratio, and distribution of hydrocarbon and heteroatom functionalities, are in good agreement with the experimental characterization results.

In fact, the organic chemical structure of kerogen is a three-dimensional (3D), highly cross-linked macromolecule comprising different functional groups. The configuration of the kerogen molecular structure in geometric space can reflect a number of important characteristics of oil shale or kerogen, such as energy structures, inter-/intramolecular interactions, bond orders, and electron density, which are extremely difficult to be detected by experiments. Recently, by utilizing material modeling technology and advanced analytical detection methods, studies on chemical reactivities of kerogen have focused on the model construction and property prediction of the thermal evolution of molecular kerogen. Based on chemical structure analysis, supplemented with fast pyrolysis results, Tong et al.[16] have proposed six constitutional isomers of two-dimensional molecular models of the Huadian oil shale kerogen. The simulated results indicated that van der Waals interactions have an important effect on the physical density of the kerogen molecular structure. Faulon et al.[17] have defined a chemical representation of the kerogen molecular model based on physicochemical analyses. Their work presents a modeling method of kerogen structure in three successive steps using classical chemical symbols. Salmon et al.[18] have developed a new method for evaluating maturation pathways for an important kerogen-forming geopolymer, involving reactive modeling based on quantum mechanics to reproduce maturation. Their study reports that maturation at early stages involves mainly reorganization of the kerogen and production of high-molecular-weight products and low amounts of oxygenated gaseous products. On the basis of the lignite structural model, Salmon et al.[19] have generated the macromodel of Morwell Brown coal using the ReaxFF reactive force field. Their work on using reactive dynamics simulations to describe the molecular processes of coal successfully presents the main reactions of thermal decomposition of the functional model. Liu et al.[20] have investigated the initial pyrolysis mechanism of oil shale kerogen using ReaxFF molecular dynamics (MD). As presented in their NVT simulations, the initial thermal decomposition of the Green River kerogen model was in good agreement with the result of gas chromatograph-mass spectrometer (GC-MS) analysis. Freund et al.[21] proposed a coupled chemical structure-chemical yield model, which was based on detailed solid-state chemical characterizations. Their model can be used to predict oil and gas compositional yields of kerogens during the thermal process.

With respect to intrinsic chemical reactivity (e.g., cleavage of chemical bonds and tendency of pyrolysis), marginal information about the electronic behavior of kerogen is available. Several relevant studies have reported on the electronic behavior by employing thermogravimetric analysis (TGA)[22−24] or reactive force field simulations[25,26] to characterize the thermal decomposition products of a kerogen sample or a kerogen molecular model, but modeling studies that describe the electron distribution in kerogen are not sufficient, which plays a key role in understanding the chemical mechanism of hydrocarbon generation. The essence of all chemical reactions in kerogen pyrolysis can be attributed to the breaking and bonding processes of chemical bonds that comprise electrons. Therefore, electron-related properties involved in the thermal conversion process of kerogen, including electron density, charge distribution, and molecular orbitals, are key characteristics for revealing and describing the chemical details of the microscopic mechanism of kerogen reactions. Simulation calculations using molecular dynamics and quantum chemistry (QC) are effective for estimating the electronic parameters of kerogen, whereby the intrinsic relationship between the chemical structure of kerogen and its microscopic chemical reactivity can be constructed to gain insights into the pyrolysis mechanism.

In this study, to obtain the missing information of the microscopic details on the intrinsic chemical reactivity associated with the electron density of the kerogen molecule, the Huadian oil shale kerogen was analyzed on the basis of experimental data; then, the molecular model of Huadian kerogen was constructed by molecular dynamics and quantum chemistry methods. Simultaneously, a modeling process for reconstructing the macromolecular structure of kerogen was proposed, which would maximally reduce the effect of kerogen isomers on energy minimization for kerogen molecular models. The detailed structure and morphology of kerogen models were created and calculated by MD simulations. In addition, the density evaluation of kerogen models was performed to obtain structural properties of the kerogen molecule. The Mayer bond order and valence for each bond in the kerogen model structure were calculated, which permit deep understanding of the chemical stability of kerogen. Finally, the molecular orbitals and electrostatic potential of the kerogen model were computed by a quantum-chemistry-based method that is available for large-scale simulations. In this study, results of experimental analyses of oil shale kerogen and studies of the structural and electronic properties of the kerogen molecular model are provided, thereby providing insights into the cleavage of chemical bonds in the kerogen molecular model and the relationship between reactivity and chemical features of kerogen.

Experiments and Methods

Samples

The oil shale sample was obtained from the Huadian mine located in Jilin, China. Before conducting demineralization, raw oil shales were ground into powder (<200 μm) and dried at 100 °C to remove water. Next, kerogen was isolated from minerals. As a standard method to isolate kerogen, HCl/HF/HCl treatment was employed.[27] To ensure that the isolated kerogen was free from HCl, HF, or other salt compounds, the kerogen sample was washed with distilled water until the pH of the washed filtrate was neutral.

Determination of Ignition Loss

To evaluate the demineralization effect, the kerogen sample (1 g) was placed into a crucible, followed by burning in a muffle furnace for 1 h at 800 °C. The ignition loss of the prepared kerogen sample was calculated by the following equationwhere m is the total weight of the crucible (g), ma and mb are the total weights of the crucible plus kerogen (g) and the total weight of crucible plus ash (g), respectively, and ω is the mass loss rate at ignition (%).

Experimental Methods

Solid-State 13C NMR

Solid-state 13C NMR spectra were recorded on a Bruker 13C CP MAS TOSS-4 mm-AVANCE NEO system, operating at 100.63 MHz at room temperature. The 13C NMR experiment was performed at a spinning speed of 15 kHz. The kerogen sample was packed into a double-resonance probe head using 4 mm ZrO2 rotors. Significantly, the most relevant molecular motions in kerogen are expected to be the rotation of CH3 groups attached to aliphatic and aromatic structures, segmental motions including the full and partial rotation of long-chain aliphatic components, and the torsional oscillation within aromatic rings.[28] The kerogen sample was characterized by a crosspolarization magic-angle spinning experiment performed at a contact time of 2 ms and a recycle delay time of 6 s. Over an acquisition time of 20 ms, a total accumulation of over 20 000 transients was made and 2000 data points were collected for the kerogen sample. To quantify the relative proportion of different carbon types in this kerogen sample, the 13C NMR spectra were curve-resolved. The definition of 13C NMR structural parameters and the chemical shift ranges are the same as those reported previously[28,29] and are listed in Table . Furthermore, based on the results given in Table , several structural parameters[30] were calculated and then used for the identification and quantification of the carbon skeletal structure in the Huadian kerogen sample (see Table ).

Table 4

Assignment of Carbon from Solid-state 13C NMR Spectra

Table 5

Structural Parameters Calculated from Solid-State 13C NMR Results

structural parameter	definition	value
aromatic carbon rate (fa)	fa = farH + farB + farC + farO	12.15
aliphatic carbon rate (fal)	fal = falM + falA + falH + falD + falO	84.66
protonated aromatic carbon rate (faH)	faH = farH	4.81
substitution of the aromatic ring (D)	D = (farB + farC + farO)/fa	60.41
alkyl-substituted carbon rate (faS)	faS = farC	3.67
bridged aromatic carbon rate (faB)	faB = farB	1.40
aromatic cluster size (XBP)	XBP = faB/faCP	13.05
alkyl-chain branching (Bi)	Bi = falD/fal	4.62
phenolic hydroxyl group or ether-oxygen-bound rate (faP)	faP = falO + farO	4.06
nonprotonated aromatic carbon rate (faN)	faN = faS + faB + faP	9.13
methylene percentage of aliphatic (Ai)	Ai = falH/fal	88.84
average carbon number of the methylene chain (Cn)	Cn = falH/farC	20.49

Fourier Transform Infrared (FTIR) Spectroscopy

FTIR spectra were recorded on a Nicolet 6700 spectrometer (Thermo Nicolet Corporation). First, the kerogen sample (2 mg) was mixed with KBr (200 mg). Second, the mixtures were ground into powder (particle size < 0.074 mm) and pressed into discs at 12 MPa for 10 min. Finally, the prepared sample was placed into a drying oven (45 °C) for 1 h to avoid the influence of the H2O absorption peak on the infrared spectra of kerogen. FTIR spectra for the sample were obtained at 4 cm–1 resolution and collected in the 4000–400 cm–1 wavenumber range. The number of scans was 32.

X-ray Diffraction (XRD)

XRD diffraction patterns of oil shale and kerogen were recorded on a D8 VENTURE X-ray diffractometer (Bruker Corporation) with the following conditions: X-rays, Cu Kα operating voltage = 40 kV; operating current = 30 mA; scanning speed = 2°/min; scanning range = 2–90°; acquisition speed = 0.3 s/step; interval step = 0.02°.

X-ray Photoelectron Spectroscopy (XPS)

XPS analysis of the Huadian oil shale kerogen sample was performed with ESCALAB 250XI X-ray photoelectron spectroscopy (Thermo Fisher Scientific Inc.), equipped with a microfocusing monochromator and a charge compensation system. The monochromatic Al Kα radiation was used, which was operated at 150 W, and the spot size was 500 μm. The correction for binding energy was made to account for kerogen sample charging based on the carbon (1s) peak at 284.8 eV. Pass energies of the survey and narrow spectra were fixed at 100 and 30 eV, respectively. The oxygen species and corresponding contents can be calculated on the basis of the XPS carbon (1s) signal of adjacent carbon atoms. Four peaks at 284.8, 286.3, 287.5, and 289.0 eV were used to curve-resolve the XPS carbon (1s) signal. Nitrogen forms and sulfur speciation in the kerogen sample were defined and quantified referring to the parameters reported by Kelemen et al.[29] These signals were simultaneously curve-resolved using components with fixed binding energy positions and the full width at half-maximum values.

Computational Methodology

Molecular Dynamics

The two-dimensional (2D) molecular structure model is the starting point for examining the microscopic reactivities of kerogen. By geometry optimization (i.e., energy minimization), the 2D molecular model of kerogen can be converted into a three-dimensional (3D) structure. The optimized 3D structure corresponds to a minimum in the potential energy surface. Therefore, geometry optimization constrains the coordinates of atoms until it satisfies certain specified criteria; meanwhile, the total energy of the molecular model is calculated. Therefore, when the stable molecular conformation with the minimum energy is output and determined, the energy minimization stage is deemed to be accomplished.

In this study, MD and QC simulation calculations were performed using Materials Studio 2017 software, and the overall quality for simulations was fine. To improve the modeling efficiency, the established kerogen model and its isomers were subjected to a geometry optimization task to search for a minimum energy structure with the following calculation details: algorithm for geometry optimization, Smart; convergence tolerance: energy < 1 × 10–4 kcal/mol; force < 0.005 kcal/mol/Å; displacement < 5 × 10–5 Å; maximum number of cycles: 100 000; calculation method for atomic charges: force field assigned. The first ab initio force field COMPASSII was selected as the high-quality force field to consolidate parameters of atoms in kerogen models.[31−33]

The atoms in the kerogen model exhibit a high degree of freedom of displacement; as a result, a minimum point in the molecular configuration within a certain displacement range may be present. However, the potential energy for the molecular configuration may be at a new minimum after a relatively large displacement. Therefore, to avoid the molecular structure falling into the local energy minimum point (potential well) during geometry optimization, anneal dynamics was adopted to search for the lowest-energy structure of the kerogen model, which allows for the conformation to overcome energy barriers to shift to other low-energy areas via the periodic control of the simulation temperature. Dynamics details: annealing cycles = 20; initial temperature = 300 K; mid-cycle temperature = 500 K; heating ramps per cycle = 50; dynamics steps per ramp = 100; time step = 1 fs; ensemble, NVT.

Quantum Chemistry Simulations

In this study, QC simulations of the kerogen molecular model were performed using the DMol[3] program based on the density functional theory (DFT). This simulation method can be employed to evaluate electron-dependent properties with respect to structural changes, which is critical to examine chemical systems and large-scale molecular structures.[34] Hence, the electronic properties of the kerogen model (Mayer bond order, molecular orbitals, and electrostatic potential) were computed using the DMol[3] module with the following calculation details: functional: GGA and B3LYP;[35,36] SCF tolerance: 1 × 10–6; basis set: DNP; max. SCF cycles: 50; multipolar expansion: hexadecapole. All electrons are included in the calculation.

Results and Discussion

Structural Characteristics of Kerogen

Ultimate Analysis

After minerals were removed, the oil shale kerogen sample was first evaluated by the determination of ignition loss. According to the calculation results of eq , the ash content and mass loss rate at ignition (ω) for the kerogen sample were 0.92 wt % and 99.08%, respectively, which can satisfy quality requirements for sample preparation.[37] Then, oil shale kerogen was characterized by ultimate analysis, as summarized in Table . The carbon content of kerogen was 74.88%. The atomic ratios H/C and O/C were calculated for classifying the type of kerogen. On the basis of the van Krevelen diagram,[38] the Huadian oil shale kerogen belonged to Type I. The H/C and O/C ratios of Huadian kerogen were close to those of Green River oil shale kerogen (H/C, 1.59 and O/C, 0.12).[39] Compared to Types II and III, the Huadian kerogen exhibited a higher potential for petroleum generation based on the method defined by Vandenbroucke and Largeau.[40]

Table 1

Ultimate Analysis of Huadian Kerogen (Dry Ash-Free Basis, wt %)

sample	C	H	Oa	N	Sb	H/Cc	O/Cc
kerogen	74.88	9.67	10.71	1.32	2.75	1.55	0.11

By difference.

Total sulfur.

Atomic ratio.

XRD Analysis

Figure shows the XRD patterns of oil shale and kerogen. The diffraction intensities reflect that the oil shale sample comprises six major types of minerals. Clay and carbonate minerals were the main minerals in terms of the relative quantity, mainly containing montmorillonite, kaolinite, illite, and calcite. In addition to the diffraction peaks of montmorillonite, each characteristic peak exhibited a fine peak shape, indicative of the presence of crystalline minerals in the sample. In the XRD spectra of kerogen, peak characteristics of other types of minerals were absent, except for a small amount of pyrite. Owing to the removal of clay and carbonate minerals, peak characteristics of pyrite were more clear and prominent than those observed for the oil shale sample. Besides, a clear broad peak was observed at 2θ = 20°, mainly corresponding to hydrocarbons, indicating that the kerogen chemical structure comprises amorphous organic matter. By comparing the XRD patterns of oil shale and kerogen samples, the minerals with a crystalline structure exerted a significant masking effect on the diffraction peak of the organic matter in the raw oil shale.

Figure 1

XRD patterns of oil shale and kerogen (M, montmorillonite; K, kaolinite; I, illite; C, calcite; Q, quartz; P, pyrite).

FTIR Analysis

Table summarizes the detailed functional groups in the kerogen sample. The position and intensity of characteristic peaks depended on the vibration form of chemical bonds and the surrounding chemical environment. According to the distributions of the characteristic groups at 4000–400 cm–1, the organic functional groups present in the kerogen sample were categorized into nine groups (Table ). In addition, the in-plane vibration region of −(CH2)– (n > 4) observed at ∼720 cm–1 reveals that the length of the methylene carbon chain is no less than four in terms of carbon number in the organic structure of kerogen. Notably, FeS2 was identified by FTIR analysis, as well as detected in the XRD spectra of kerogen. As the major residual mineral that is frequently present as the constituent of kerogen, pyrite is concentrated after chemical elimination (i.e., acid treatment). Typically, oxidative reagents such as[41] HNO3 or FeSO4 are used to remove pyrite. This treatment can lead to the oxidation of the organic matter functional groups in kerogen.[42] Meanwhile, by the incorporation of sulfur-bearing reagents, the sulfur content of kerogen increases undesirably. In view of this, in this study, a HNO3 or sulfur-bearing reagent was not used for the preparation of kerogen, and the pyrite component also was not involved in the kerogen model.

Table 2

FTIR Analysis of the Functional Groups in Kerogen

functional group	wavenumber (cm–1)	functional group	wavenumber (cm–1)
FeS2	425	–C=C–	1622
–(CH2)n–	720	–C=O–	1706
–C–O–, −C–OH	1286	–CH2	2846
–CH3	1363	–CH2	2925
–CH3, −CH2	1459	–OH	3410

XPS Analysis

Table summarizes the identification results for the organic forms of oxygen, nitrogen, and sulfur in kerogen. XPS is capable of providing detailed information about the amount and form of aromatic carbons present in complex carbonaceous materials. Significantly, XPS is a surface method and does not provide a true measure of the bulk moiety distribution. The oxygen species and corresponding contents can be calculated on the basis of the XPS C 1s peak.[43] The relative content of each functional group is defined as the mole fractions of organic C, N, and S, which corresponds to the ratio of each peak area to the total area in the XPS spectrum.[44] According to the curve fitting results of the XPS C 1s spectrum, aromatic and aliphatic carbons were the main forms of organic carbon in kerogen with a relative content of 88.44 mol %. Ether and hydroxyl groups (C–OH, O–C–O) constituted the main components, with the highest mole percent of 8.06%. Carboxyl (O=C–O) and carbonyl (O=C) groups exhibited relatively lower contents of 2.42% and 1.08 mol %, respectively, reflecting that Huadian kerogen is immature organic matter formed in the first stage of deposition. For N 1s peaks, nitrogen was mostly present as pyrrolic and amine functional groups, followed by pyridinic nitrogen. The pyrrolic and pyridinic N structures were thought to form and become predominant during the petroleum formation stage, generating soluble nitrogenous compounds in oil products.[45,46] Moreover, a low content of nitrogen oxides also was detected (10.60%). With the increasing grade of maturity, the nitrogen species of kerogen changed. As a result, pyridinic nitrogen compounds become the main form of nitrogen and are frequently found in immature sedimentary organic matter. In this study, immature Huadian type I kerogen comprised 40.94 mol % pyridinic nitrogen, which corresponded to the highest amount of N species in the kerogen sample. This result is in agreement with previously reported studies.[47] Aromatic sulfur, aliphatic sulfur, thiophene, sulfoxide, and pyrite were found in sulfur. Wiltfong et al.[48] have characterized sulfur speciation in six different kerogens by XANES spectroscopy and suggested that oxidized sulfur forms are more prevalent in the Type I kerogens, which is in agreement with the Huadian kerogen sample used in this study. Thiophene as well as aromatic and aliphatic sulfur corresponded to organic sulfur fractions, accounting for 7.89 and 21.94 mol %, respectively. Thiophenes consist of sulfur in an aromatic carbon ring and are relatively thermally stable.[49] According to our earlier study,[50] the linear carbon skeletons in aliphatic sulfur compounds tended to undergo aromatization, affording short linear alkyl (C1–C4) thiophene compounds during the early stage of sediment burial and secondary thermal reactions; hence, aliphatic sulfur considerably contributes to the formation of thiophene units in kerogen. Apart from organic sulfur, pyrite was also detected in the kerogen sample. The quantitative separation of pyrite-sulfur from organic sulfur is a difficult problem.[51] Pyrite was mainly present in the kerogen sample as it was extremely difficult to be completely removed by acid treatment or physical methods; thus, a certain amount of pyrite is present in the organic structure of kerogen and exists as inorganic sulfur. For the sulfoxide groups in kerogen, according to research reported by Pomerantz et al.,[52] during cracking of kerogen, weak bonds between carbon and reduced forms of sulfur crack readily to form reactive sulfur radicals, and then, these sulfur radicals can combine with nearby oxygen sources such as water or organic oxygen to form sulfoxide groups.

Table 3

Organic Forms of Oxygen, Nitrogen, and Sulfur in Kerogen

3.1.5 Structural Parameters of Carbons

Figure shows the solid-state 13C NMR spectra for kerogen. Table summarizes the identification results of carbons from solid-state 13C NMR spectra. Structural parameters of carbons in the kerogen sample were classified on the basis of chemical shifts, and the corresponding relative content of each carbon functional group was calculated by area normalization. The highest amount of methylene (falH) carbon, accounting for 75.21%, was observed for the kerogen sample, indicating that aliphatic hydrocarbon compounds are the dominant type of carbons. In comparison, the relative content of aromatic methyl (falA) carbon was 2.19%, reflecting that only a small amount of aromatic carbons is present in the organic structure of kerogen. The proportional relationship between aliphatic and aromatic carbons is a typical feature of type I kerogen. The Huadian kerogen sample contained a certain amount of methine, quaternary carbon (falD = 3.91%), indicating that tertiary or quaternary carbon atoms are present in the organic carbon skeleton connecting to the aliphatic carbon chains and subsequently leading to the formation of a cross-linked structure of the kerogen molecule. The carboxyl (faC) and carbonyl (faO) carbons accounted for 2.42 and 0.76%, respectively, indicating that C=O double bonds are mainly present in kerogen as carboxyl functional groups in comparison to carbonyl functional groups. Meanwhile, the result indicated that the carbonyl carbon accounts for the lowest proportion, which is in good agreement with the results obtained from the XPS C 1s peak.

Figure 2

Solid-state 13C NMR spectra for kerogen.

Owing to the complexity and heterogeneity of the kerogen structure, it is still difficult to accurately understand its molecular structure. Therefore, structural parameters calculated by mathematical statistics are often used to comprehensively describe the average structural characteristics of the kerogen macromolecule. In this study, statistical constitution analysis was employed to calculate the statistical structural parameters of Huadian kerogen. Table lists the corresponding results. According to the structural parameters calculated from solid-state 13C NMR results, organic carbon structure parameters in kerogen were categorized into three groups: aromatic carbons, aliphatic carbons, and oxygen-containing carbons. For the same kerogen sample, the aromatic carbon ratio (fa) was inversely proportional to the aliphatic carbon rate (fal), which can reflect the content distributions of aromatic and aliphatic compounds; meanwhile, fa and fal values were affected by the maturity of kerogen. The alkyl-substituted carbon rate (faS) and the bridged aromatic carbon rate (faB) accounted for 3.67 and 1.40%, respectively, indicating that the aromatic structure units are mainly connected to the alkyl functional groups, i.e., alkyl-substituted groups compared to aryl-substituted groups. Aromatic cluster size (XBP) reflects the ratio of aromatic bridge carbons to the total aromatic carbons, which is used to determine the degree of aromatic condensation in kerogen. The XBP of the Huadian kerogen sample was 13.05%; this value is approximately close to that of naphthalene (XBP = 20%), while it is far less than that of anthracene (or phenanthrene, XBP = 30%). This result suggests the presence of 1 or 2 aromatic rings in the kerogen sample, namely, benzene and naphthalene ring structures. Hence, a ternary aromatic-ring structure is unlikely to be present in Huadian kerogen or is present in an extremely low amount.

Determination of Molecular Structural Units

According to the experimental results with respect to structural characteristic parameters, the molecular structure of kerogen mainly comprised aromatic compounds, heteroatom-containing organic compounds, and hydrocarbons with different saturations. Aromatic and cyclic structure clusters in the kerogen molecular model were connected by alkyl/olefinic chains, with different lengths and branching extents, finally affording a complete kerogen molecule. The specific construction process was as follows.

First, the five elements of C, H, O, N, and S, as well as their corresponding relative contents, were determined on the basis of the results of ultimate analysis. However, the kerogen molecules do not exhibit a fixed or repetitive monomer structure. Therefore, it is crucial to determine the molecular mass of the kerogen model. Although a model with a low molecular weight can be easily calculated for quantum chemistry and molecular dynamics simulations, notably, the molecular weight should not be extremely low; the molecular weight is supposed to be sufficiently high to gain a statistical model that can reflect the reasonable complexity of the real kerogen structure. In view of this, the initial chemical formula of the Huadian kerogen model was determined to be C233H357N3O25S3 (relative molecular mass, 3691). The level of the selected molecular size was comparable to those reported previously.[53,54]

For the structural distribution of carbon, the structure and content of aromatic compounds in the kerogen model were determined by fa, faH, faS, and XBP, which was calculated by 13C NMR experimental data. The aromatic structural units were randomly arranged and connected by alkane or olefin chains. The carbon chains of various lengths were designed on the basis of Bi, fal, Ai, and Cn.

For the distribution of functional groups, according to FTIR analysis results, C–O, −OH, C=C, C=O, and −CH3 also were present and differently added to the kerogen model for the construction of aliphatic chains on the basis of the position of characteristic functional groups. In addition, each functional group can be involved in the construction of aliphatic and aromatic structures to ensure that the model values are consistent with the experimental values.

In the kerogen model, sulfur and nitrogen species were incorporated as heterocycles, such as pyrrolic, pyridinic, thiophene, and sulfoxide compounds, according to the chemical parameters obtained by XPS analysis. Meanwhile, the relative contents of sulfur and nitrogen compounds were limited by ultimate analysis, which were supposed to be consistent with the elemental compositions of the kerogen sample.

In fact, the kerogen molecular structure model exhibited a three-dimensional configuration in space. Owing to the presence of a large number of isomers, there are a large number of molecular configurations that satisfy the same structural parameters of the kerogen sample. Accordingly, for a kerogen model, the molecular configuration with the minimum energy obtained by geometry optimization is not necessarily the structure with a reasonable geometry because the optimized structure corresponding to the minimum in the potential energy surface may exist in other kerogen isomers. In view of this, a design concept of the grid method was adopted herein. The grid method involves changing the network structure of the molecular structure of the kerogen model, i.e., to build a certain number of closed grid spaces in the molecular structure, and then, the molecular configuration with the lowest energy was selected. This method can effectively avoid a considerable number of isomers, and the details of the method used here are described in our earlier study.[55]Table lists the resulting distribution and number of carbon-containing structures in our kerogen model, as well as the comparison between model and experimental values.

Table 6

Main Structural Parameters of the Kerogen Molecular Model

carbon types	kerogen model	parameters	model values	experimental values
aliphatic methyl	5	fa	14.02	12.15
methylene	132	fal	81.97	84.66
oxy-aliphatic carbon	16	Cn	18.13	20.49
methine, quaternary	19	D	66.66	60.41
bridged aromatic carbon	8	H/C	1.54	1.55
alkyl-substituted aromatic carbon	13	O/C	0.11	0.11
oxygen-substituted aromatic carbon	7	C	75.71	74.88
carboxyl carbon	2	N	1.14	1.32
carbonyl carbon	2	S	2.60	2.75

Structural Properties of Kerogen

First, models of isomers were subjected to geometry optimization; then, the lowest-energy structures were output by dynamic annealing simulations. Figure shows the corresponding model structures. Owing to the influence of the number of grids, 3D molecular structures of kerogen exhibited various cross-linking degrees during energy minimization. After dynamic optimization and calculation, the molecular models of kerogen became dense and compact, exhibiting anisotropic spatial stereo structures. The aromatic hydrocarbon structure was buried in aliphatic compounds and randomly disordered in a spatial arrangement. The aromatic units were connected by methylene chains with different lengths to form a network skeleton, which was conducive to the storage of small organic molecules or free compounds in the molecular structure of kerogen. Notably, chainlike structures in the models, especially aliphatic chains, were highly twisted and curled; meanwhile, end-carbon chains connected with the kerogen molecular structure tended to undergo cyclization during energy minimization. In addition, the highly cross-linked molecular structure in space is one of the main reasons for the high insolubility of kerogen in general solvents. At the same time, as aromatic or ring structure units are surrounded by aliphatic carbon chains, the aromatic-ring structure cannot easily dissolve during solvent extraction; hence, extraction products are mainly hydrocarbons with a low carbon number or free small-molecule compounds.

Figure 3

Three-dimensional models of geometric isomers for the Huadian kerogen after geometry optimization. The number of grids in the isomer models ranges from 2 to 13. The 12 kerogen isomer models (Grid-2 to Grid-13) have the same molecular formula of C233H357N3O25S3. The atom colors are as follows: C, gray; O, red; N, blue; S, yellow; H, white.

To investigate the influence of grid number and isomers on the energy of the kerogen molecular model, energy variation curves of each isomer model with different grid numbers are plotted in Figure . With the increase in the number of grids, the corresponding energy of the kerogen molecular model fluctuated within a certain range. The energy of kerogen models and the number of grids did not exhibit a linear relationship. Overall, the energy change exhibited a decreasing trend first and then an increasing trend. A unique inflection point was present in the molecular structure model of kerogen under a certain number of grids, i.e., the energy minimum point, which corresponded to an isomer structure of kerogen. For the kerogen molecular system, the lower the structural energy of the isomer, the higher the molecular stability. When the grid number of the kerogen molecular model was 9, the corresponding isomer exhibited the lowest energy. Consequently, the isomer with nine grids (Grid-9) was selected as the final molecular model for the Huadian kerogen sample. Meanwhile, the final kerogen model was subjected to dynamic calculations. Table lists the detailed potential energy components.

Figure 4

Total energy for the kerogen models with different grids.

Table 7

Potential Energy Components of the Kerogen Model

contributions to the total energy (kcal/mol)	energy components
valence interactions (Evalence)	399.120
bond stretching (Ebond)	225.199
valence angle bending (Eangle)	454.311
dihedral angle torsion (Etorsion)	–286.660
inversion (Eoop)	6.269
valence cross terms (Ecrossterm)	–123.652
stretch–stretch (Es-s)	0.797
stretch–bend–stretch (Es-b-s)	–15.996
stretch–torsion–stretch (Es-t-s)	–17.876
separated–stretch–stretch (Esp-s-s)	0.927
torsion–stretch (Et-s)	–21.657
bend–bend (Eb-b)	0.829
torsion–bend–bend (Et-b-b)	–24.941
bend–torsion–bend (Eb-t-b)	–45.734
nonbonding energy (Enon-bond)	–464.507
van der Waals (EvdW)	–131.005
electrostatic (ECoulomb)	–333.502
total energy (Etotal)	359.966

The total energy of the kerogen molecular model (Etotal) mainly comprised bonding energy (Evalence + Ecrossterm) and nonbonding energy (Enon-bond). Evalence made the highest contributions to the Etotal of the kerogen model (399.120 kcal/mol). The Evalence value reflects that the organic structural units in the kerogen molecule are connected by covalent bonds, which played a key role in stabilizing the molecular structure. Energy composition reveals that the energy optimization process for the kerogen molecular structure is mainly caused by the change in the bond angle caused by atom displacement, leading to the variation in the valence angle bending (Eangle). The Enon-bond of the kerogen model was −464.507 kcal/mol. The electrostatic energy ECoulomb played a dominant role in the nonbonding energy of the kerogen molecular system, which was −333.502 kcal/mol. This result is related to the fact that when the molecular system exhibits multiple charges or multipole moments, the electrostatic interactions within the molecule will become stronger.[56] Therefore, the ECoulomb of the kerogen molecular model accounts for a large proportion. Compared with ECoulomb, the van der Waals (EvdW) in the model was relatively small, which was −131.005 kcal/mol (negative energy indicates mutual attraction). When π bonds are present in molecular structures, such as delocalized π bonds in benzene, the intermolecular interactions increase. Kerogen is a type of organic matter with low maturity. The content of aromatic components was less than that of aliphatic components, and the aromatic-ring structures were wrapped by different lengths of methylene chains. Therefore, the EvdW of the kerogen molecule is small due to the spatial hindrance effect.

Multimolecular structure models of kerogen with different densities were established to evaluate the accuracy of the constructed kerogen model. Owing to the influence of maturity and evolution origin, the density of actual kerogen typically ranges from 0.79 to 1.51 g/cm3.[40,57] Therefore, 18 multimolecular structure models of kerogen with density values from 0.40 to 1.30 g/cm3 are constructed. Figure shows the variation tendency of the total energy of the kerogen model at different densities. With the increase in density, the energy of the multimolecular structure model exhibited a decreasing trend first and then an increasing trend, with an inflection point in the density-energy diagram. At a density of 0.95 g/cm3, the condensed kerogen model exhibited the lowest energy of 372.404 kcal/mol. For comparison, the experimental density value for the Huadian kerogen sample used herein was measured by the specific gravity bottle method. As a result, the simulated density of the kerogen model is in good agreement with the experimental density, 0.97 g/cm3, indicating that the constructed Huadian kerogen model in our study can reasonably reflect the true physical density of the kerogen sample.

Figure 5

Total energy of the kerogen model at different densities (a). Final 3D model of kerogen with a density of 0.95 g/cm3 (b). Atom colors are the same as those shown in Figure .

Mayer Bond-Order Analysis

For a deeper understanding of the chemical reactivity of the kerogen sample, the cleavage behavior of bonds was investigated on the basis of the Huadian kerogen model (Figure ). The Mayer bond orders of all chemical bonds in the kerogen structure were calculated and analyzed by Hirshfeld atomic charges. Table summarizes the calculation results for the Mayer bond orders. The Mayer bond-order analysis estimates valences that are close to the corresponding classical values; therefore, the Mayer bond orders can be used to describe binding characteristics in the same or similar molecular system.[58] The distributions of Mayer bond orders in the kerogen model varied with the types of chemical bonds, ranging from 0.886 to 1.827. For the carbon skeleton structures, the bond orders of C–C single bonds mainly focused on 1.000. The bond-order values of aromatic(C)–aromatic(C) and alkane(C)–alkane(C) were ∼1.300 and 0.900, respectively, indicating that the aromatic structures are more stable than aliphatic structures due to the presence of π bonds in a conjugate system. Furthermore, for hydrocarbon structures, the bond-order value of alkene(C)=alkene(C) was 1.688–1.870, which was significantly greater than that of alkane(C)–alkene(C) (0.985 to 1.024); this bond-order value reflects a remarkable influence in that the C=C double bond can strengthen the stability of carbon chains in the kerogen sample. Furthermore, the general average value of bond orders in alkene(C)–alkene(C) (1.100) was slightly greater than that of alkane(C)–alkene(C) (1.005). This is another evidence for the fact that the C=C bond plays a key role in the stability of the molecular structure. Meanwhile, the bond-order value of the C–C bond increased due to the number of surrounding C=C bonds, i.e., olefinic compounds. In conjugated systems with alternating C–C and C=C bonds, although the bond order between atoms is not exactly the same, the difference in the bond-order value between single and double bonds is reduced due to electron delocalization, leading to average bond-order values. The longer the carbon chain in conjugated systems, the lower the difference in the bond-order value between single and double bonds. In addition, the energy of the conjugated system decreased due to electron delocalization. Hence, the conjugated system (ring-shaped conjugated structures and linear-chain conjugated structures) is more stable than the nonconjugated system in the kerogen molecular model. For the sulfur-containing chemical bonds, S=O in the sulfoxide structure exhibited the highest bond-order value (1.827). The bond-order values for S–C bonds (1.188 and 1.184) in thiophene were similar, which was related to the effect of conjugated π bonds. In contrast, the bond orders of sulfoxide(S)–alkane(C) and sulfydryl(S)–alkane(C) were significantly lower, indicating that the sulfur-containing groups tend to break at the junction with alkane carbon atoms and released as −SH and −SO groups. For nitrogen-containing chemical bonds, the bond orders of N–C in pyridine are greater than those in pyrrole, suggesting that six-ring structures are more stable than five-ring structures in terms of nitrogen heterocyclic compounds.

Figure 6

Two-dimensional model of Huadian kerogen (Grid-9). The atoms are to be sequentially numbered in Arabic numerals.

Table 8

Calculation Results for Mayer Bond Orders

chemical type	bond type	bond order	chemical type	bond type	bond order
sulfoxide(S)–sulfoxide(O)	S3 = O1	1.827	alkane(C)–alkane(C)	C2–C193	0.973
sulfoxide(S)–alkane(C)	S3–C4	0.975		C4–C5	0.981
	S3–C2	0.976		C5–C6	0.990
sulfydryl(S)–alkane(C)	S187–C186	1.053		C7–C8	0.995
thiophene(S)–thiophene(C)	S219–C213	1.188		C8–C9	0.985
	S219–C220	1.184		C10–C11	0.972
pyridine(N)–pyridine(C)	N63–C64	1.413		C11–C13	0.980
	N63–C62	1.390		C13–C14	1.004
pyrrole(N)–pyrrole(C)	N155–C151	1.106		C14–C15	0.986
	N155–C156	1.161		C15–C16	1.002
amide(N)–alkane(C)	N224–C223	1.272		C16–C17	1.003
ether(O)–ether(C)	O88–C89	0.982		C17–C18	0.989
	O91–C89	0.996		C18–C19	0.998
	O91–C98	0.886		C19–C20	0.953
	O33–C32	0.915		C29–C30	0.975
	O33–C34	0.911		C34–C37	0.962
	O36–C34	0.936		C37–C38	0.987
ether(O)–aromatic(C)	O40–C24	1.032		C38–C39	1.006
	O113–C119	1.031		C42–C43	0.994
	O212–C218	0.988		C50–C51	0.991
oxhydryl(O)–C	O90–C89	1.059		C52–C53	0.984
	O48–C47	0.956		C52–C76	0.969
	O161–C116	1.040		C53–C54	0.993
oxhydryl(O)–aromatic(C)	O68–C67	1.056		C55–C56	0.978
	O158–C154	0.978		C55–C58	0.963
carboxyl(O)–carboxyl(C)	O74–C73	1.172		C58–C59	0.997
	O75 = C73	1.870		C72–C73	0.967
	O136–C135	1.159		C76–C77	0.990
	O139 = C135	1.839		C93–C12	1.008
carbonyl(O)–carbonyl(C)	O225 = C223	1.607		C94–C95	0.992
	O104 = C103	1.731		C97–C98	1.009
oxhydryl(O)–cyclane(C)	O131–C121	0.972		C106–C107	0.988
alkene(C)=alkene(C)	C84 = C85	1.701		C109–C110	0.994
	C125 = C126	1.688		C132–C133	0.977
	C143 = C144	1.717		C133–C134	0.960
	C145 = C146	1.697		C133–C138	0.921
	C167 = C168	1.836		C147–C170	0.976
	C201 = C202	1.804		C171–C172	0.993
	C210 = C211	1.787	pyridine(C)–alkane(C)	C64–C99	0.959
	C227 = C228	1.809		C66–C69	0.961
	C243 = C244	1.813		C61–C62	0.985
alkane(C)–aromatic(C)	C10–C254	0.964	alkene(C)–alkene(C)	C83–C84	1.145
	C26–C29	0.969		C126–C127	1.142
	C153–C159	0.976		C144–C145	1.109
	C215–C222	0.990		C242–C243	1.002
	C216–C234	0.981	alkane(C)–alkene(C)	C85–C86	1.008
pyrrole(C)–pyrrole(C)	C156–C157	1.488		C124–C125	1.001
pyrrole(C)–alkane(C)	C156–C178	0.975		C142–C143	1.024
pyridine(C)–pyridine(C)	C62–C67	1.279		C146–C147	0.992
	C64–C65	1.361		C200–C201	1.006
	C65–C66	1.333		C202–C203	1.021
	C66–C67	1.252		C226–C227	0.985
thiophene(C)–thiophene(C)	C214–C221	1.183	aromatic(C)–aromatic(C)	C22–C23	1.260
	C220–C221	1.494		C24–C25	1.308
cyclane(C)–cyclane(C)	C114–C120	0.969		C25–C26	1.361
	C115–C123	0.986		C26–C27	1.278
	C121–C122	0.952		C23–C27	1.257
	C122–C123	0.968		C149–C150	1.238
cyclene(C)–cyclene(C)	C259–C260	1.802		C152–C157	1.211
	C259–C264	1.002		C153–C154	1.297
	C260–C261	0.998		C213–C214	1.191
	C261–C262	0.978		C216–C217	1.353
	C263–C264	0.985		C253–C254	1.399
	C256–C257	1.059		C254–C255	1.421

Figure shows the distribution regions of the Mayer bond order for different types of chemical bonds in the kerogen model. The cyclene(C)–cyclene(C) and carboxyl(O)–carboxyl(C) bonds exhibited a larger distribution range of bond orders compared with other types of chemical bonds. The C–C bonds in aliphatic chains and cyclohydrocarbon compounds, including alkane(C)–alkane(C), alkane(C)–alkene(C), and cyclane(C)–cyclane(C), were mainly located on the left side of the distribution diagram and with bond-order values of less than 1.1. In contrast, the bond-order values of aromatic structures such as pyrrole, pyridine, and thiophene were greater than 1.1, representing that the aromatic structures are more stable during the thermal process compared to hydrocarbon structures. In general, the chemical bonds with a narrow bond-order distribution range are mainly located in the low-bond order region (at ∼1). The chemical bonds with a large bond-order distribution range exhibited stronger chemical stability due to the influence of the chemical environment, which were prone to multistep reactions during thermal decomposition.

Figure 7

Distribution regions of the Mayer bond order for different types of chemical bonds.

The electronic structure of the kerogen molecular model is crucial as it can provide chemical details for the prediction of chemical bond cleavage in different reaction environments as well as potential chemical reaction sites. Hence, for a quantitative description of charge distributions in our kerogen model, the Hirshfeld atomic charges of all atoms were calculated. Figure shows the corresponding distribution trends. The Hirshfeld charge analysis[59] is a series of algorithms to calculate the interaction energy between parts of the same molecular system and atomic charge densities, which is defined relative to the deformation density. Most of the carbon atoms exhibited negative electricity, and the Hirshfeld atomic charges of those concentrations were distributed at around −0.4 (Figure ). In contrast, hydrogen atoms only contain positive charges, indicating that hydrogen atoms are involved as hydrogen free radicals as well as electron acceptors in thermal decomposition reactions. For the heteroatoms (e.g., O, N, and S) in the kerogen molecule, all N and O atoms exhibited negative charges; meanwhile, O atoms exhibited stronger nucleophilicity compared to N and S atoms, indicating that the oxygen-containing functional groups as electron donors are prone to nucleophilic reactions. Notably, the S atoms in thiophene and S=O structures exhibited positive charges. On the contrary, for the −SH group, the concentration of electrons on the S atom led to the electronegativity of the sulfhydryl group, suggesting that the sulfhydryl group is more prone to electrophilic reactions, including attack by hydrogen radicals.

Figure 8

Distribution of Hirshfeld atomic charges for all atoms in the kerogen model.

Electronic Properties of Kerogen

Molecular orbitals refer to a single-electron wave function in a molecule, which is typically formed by the linear superposition of different atomic orbitals. The shape and composition of molecular orbitals play a key role in the qualitative analysis of the electronic properties, bonding properties, intermolecular interactions, and chemical activity of molecules. Analysis of the frontier molecular orbitals (i.e., the highest occupied orbital, HOMO, and the lowest unoccupied orbital, LUMO) is especially crucial to evaluate the chemical reactivity of the kerogen molecule. In view of this, the orbitals of the Huadian kerogen model were calculated using the Dmol[3] module in Materials Studio 2017 software, and then, the HOMO and LUMO were created on the isosurface of the kerogen molecular model.

The HOMO (EHOMO) and LUMO (ELUMO) energy values for the kerogen molecule were −4.4022 and −4.3138 eV, respectively (Figure ). The HOMO–LUMO energy gap value (Egap = EHOMO – ELUMO) was 0.0884 eV. Molecules with a large Egap are generally stable or unreactive, while molecules with a small Egap are reactive. According to the HOMO–LUMO energy gap value of the kerogen model, the EHOMO and ELUMO of kerogen are close, i.e., the overlap range is large, indicating that the molecular structure of Huadian kerogen is prone to intramolecular reactions. The HOMO was located in the chain alkane structure unit, while the LUMO was mainly formed by the benzene ring with branched hydrocarbon groups. Furthermore, the conjugation degree of the benzene ring structure in HOMO was significantly greater than that of carbon chains in LUMO. The EHOMO and ELUMO of kerogen can characterize its electron-donating and electron-accepting abilities, respectively. As a result, the reaction sites of kerogen are more concentrated in the chain alkane and benzene ring structure. For the HOMO, the bonding orbitals of carbon atoms were mainly sp hybrid orbitals. Affected by the conjugate effect, σ orbitals were formed in the benzene ring structure.

Figure 9

HOMO and LUMO isosurfaces of the kerogen model. The atom colors are the same as those shown in Figure .

To visually observe the distribution characteristics of electrons in kerogen, the MEP of a complete kerogen molecule was calculated. Electrostatic potential is defined as a unit of positive charge to detect the size of electrostatic interactions around the molecule, which can directly reflect the characteristics of the charge distribution of the molecule. MEP reflects the charge distribution of the entire molecule (including the nucleus and electrons), while the electron density reflects the local distribution of electrons in the molecule. The MEP diagram is beneficial for a qualitative prediction of specific chemical properties of the kerogen molecule, e.g., the electrophilic or nucleophilic reaction site. The 3D MEP diagram of the Huadian kerogen model was calculated and mapped to the isosurface of the total electron density (Figure ). The blue and red regions in the MEP diagram represent negative and positive electrostatic potentials, respectively. MEP was positive in the inner region of the kerogen molecule because of the higher influence of the positive nuclear charge, while the MEP of the external surface of the kerogen molecule exhibited positive and negative electrostatic potentials. When the kerogen molecules were close to each other during thermal decomposition, the spatial distribution of the electrostatic potential plays a key role in controlling the reaction mechanisms of active functional groups, i.e., electrophiles always tend to attack the most negative sites (blue regions in Figure ) on the kerogen molecule surface, while the action mechanism of nucleophiles was opposite to electrophiles. As shown in the MEP diagram, the red regions were mainly concentrated on the surface of oxygen atoms in the kerogen molecule, indicating that the oxygen atoms exhibit a low ability to attract electrons. This result is consistent with those reported in our previous study in that oxygen atoms exhibit a higher electron density compared with that of other atoms in kerogen.[60] Notably, the MEP of the benzene ring structure (black circle areas in Figure ) revealed a clear and uniform distribution due to the effect of conjugated π bonds, indicating that the benzene ring structure is relatively stable and that the functional groups attached to the benzene structure are more likely to be reaction sites. This similar feature is also reflected in the result of Mayer bond-order analysis (Table ), where the bond-order values of aromatic(C)–aromatic(C) are greater than those of alkane(C)–aromatic(C) in the kerogen model.

Figure 10

(a–c) MEP diagrams of the kerogen model. The atom colors are the same as those shown in Figure .

Conclusions

In this study, based on the experimental results obtained from solid-state 13C NMR, FTIR, XPS, and XRD, the structural characteristics of Huadian kerogen were analyzed, and then, the molecular model of Huadian kerogen was constructed using MD and MC methods. A modeling process for reconstructing the macromolecular structure of kerogen was adopted, which would maximally reduce the effect of kerogen isomers on energy minimization for kerogen molecular models. The detailed structure and morphology of kerogen models were created and calculated by a series of MD simulations. To investigate the distribution characteristics of electrons in kerogen, the MEP of a complete kerogen molecule was calculated. In addition, the Mayer bond-order analysis and molecular orbital were conducted to estimate the electronic properties of kerogen.

Structural parameters calculated by 13C NMR spectra revealed that the XBP of the Huadian kerogen sample is 13.05% and that one or two aromatic rings are present in the kerogen sample; hence, the ternary aromatic-ring structure is unlikely to be present in Huadian kerogen. Then, a design concept of the grid method was adopted for the construction of 3D isomer models of kerogen molecules. The energy composition of the kerogen model reveals that energy optimization is mainly caused by the change in the bond angle caused by atom displacement, leading to the variation in Eangle. The multimolecular structure models of kerogen with different densities were established to evaluate the accuracy of the constructed kerogen model; hence, the simulated density of the kerogen model shows a good agreement with the experimental value. The Mayer bond-order analysis indicates that the aromatic structures are more stable than aliphatic structures due to the presence of π bonds in a conjugate system. By the analysis of Hirshfeld atomic charges, the S atoms in thiophene and S=O structures exhibited positive charges. On the contrary, for the −SH group, the concentration of electrons on the S atom led to the electronegativity of the sulfhydryl group, suggesting that the sulfhydryl group is more prone to electrophilic reactions. According to the Egap of the kerogen model, the molecular structure of Huadian kerogen is prone to intramolecular reactions. To investigate the distribution characteristics of electrons in kerogen, the MEP of the kerogen molecule was calculated. MEP was positive in the inner region of the kerogen molecule, while the MEP of the external surface of the kerogen molecule exhibited positive and negative electrostatic potentials. Electrophiles tended to attack the most negative sites on the kerogen molecule surface, while the action mechanism of nucleophiles was opposite to electrophiles.

This paper provided results of experimental and simulation analyses of kerogen, providing insights into the cleavage behavior of chemical bonds in the kerogen model. Moreover, based on the theoretical study, the relationship between reactivity and chemical features of kerogen also was established. Notably, this is the first report presenting an MEP diagram of the kerogen model, which can provide valuable information on the determination of electrophilic or nucleophilic reaction sites for kerogen. Hence, this study will be beneficial to understand the reactivity of the kerogen structure and to subsequently enrich the study of electronic properties of the kerogen molecule.