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 (Scontaining) 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, 13 C 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 13 C nuclear magnetic resonance ( 13 C NMR) and 29 Si 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 c o n s t r u c t e d w i t h t h e e m p i r i c a l f o r m u l a s o f C 100 H 160 N 2.2 S 0.68 O 9. 22 and C 100 H 155 N 3.10 S 0.62 O 2.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 threedimensional (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 twodimensional 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
2.1. 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.
2.2. 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 equation where m is the total weight of the crucible (g), m a and m b 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 (%). 2.3. Experimental Methods. 2.3.1. Solid-State 13 C NMR. Solid-state 13 C NMR spectra were recorded on a Bruker 13 C CP MAS TOSS-4 mm-AVANCE NEO system, operating at 100.63 MHz at room temperature. The 13 C 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 ZrO 2 rotors. Significantly, the most relevant molecular motions in kerogen are expected to be the rotation of CH 3 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 13 C NMR spectra were curve-resolved. The definition of 13 C NMR structural parameters and the chemical shift ranges are the same as those reported previously 28,29 and are listed in Table 4. Furthermore, based on the results given in Table 4, 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 5).
2.3.2. 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 H 2 O 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.
2.3.4. 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 halfmaximum values.
2.4. Computational Methodology. 2.4.1. 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.
2.4.2. 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
3.1. Structural Characteristics of Kerogen. 3.1.1. 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 1, 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 1. 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 3.1.2. XRD Analysis. Figure 1 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.
3.1.3. FTIR Analysis. Table 2 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 2). In addition, the in-plane vibration region of −(CH 2 ) n − (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, FeS 2 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 HNO 3 or FeSO 4 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 HNO 3 or sulfur-bearing reagent was not used for the preparation of kerogen, and the pyrite component also was not involved in the kerogen model.
3.1.4. XPS Analysis. Table 3 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 ACS Omega http://pubs.acs.org/journal/acsodf Article 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 (C 1 −C 4 ) 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. 3.1.5 Structural Parameters of Carbons. Figure 2 shows the solid-state 13 C NMR spectra for kerogen. Table 4 summarizes the identification results of carbons from solidstate 13 C 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 (f al H ) carbon, accounting for carbons. In comparison, the relative content of aromatic methyl (f al A ) 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 (f al D = 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 (f a C ) and carbonyl ( f a O ) 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.
Owing to the complexity and heterogeneity of the kerogen structure, it is still difficult to accurately understand its molecular structure. Therefore, structural parameters calcu-  ) and the bridged aromatic carbon rate (f a B ) 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 (X BP ) 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 X BP of the Huadian kerogen sample was 13.05%; this value is approximately close to that of naphthalene (X BP = 20%), while it is far less than that of anthracene (or phenanthrene, X BP = 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.
3.2. 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.
(a) 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 C 233 H 357 N 3 O 25 S 3 (relative molecular mass, 3691). The level of the selected molecular size was comparable to those reported previously. 53 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 6 lists the resulting distribution and number of carbon-

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 3 shows the corresponding model structures. Owing to the influence of the number of grids, 3D molecular structures of kerogen exhibited various crosslinking 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.
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 4. 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 7 lists the detailed potential energy components.
The total energy of the kerogen molecular model (E total ) mainly comprised bonding energy (E valence + E crossterm ) and nonbonding energy (E non-bond ). E valence made the highest contributions to the E total of the kerogen model (399.120 kcal/mol). The E valence 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 (E angle ). The E non-bond of the kerogen model was −464.507 kcal/mol. The electrostatic energy E Coulomb 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 E Coulomb of the kerogen molecular model accounts for a large proportion. Compared with E Coulomb , the van der Waals (E vdW ) 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 E vdW 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/cm 3 . 40,57 Therefore, 18 multimolecular structure models of kerogen with density values from 0.40 to 1.30 g/cm 3 are constructed. Figure 5 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/ cm 3 , 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/cm 3 , indicating that the constructed Huadian kerogen model in our study can reasonably reflect the true physical density of the kerogen sample.
3.4. 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 6). The Mayer bond orders of all chemical bonds in the kerogen structure were calculated and analyzed by Hirshfeld atomic charges. Table 8 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

ACS Omega
http://pubs.acs.org/journal/acsodf Article 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 7 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 bondorder 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

ACS Omega
http://pubs.acs.org/journal/acsodf Article environment, which were prone to multistep reactions during thermal decomposition. 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 8 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 8). 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.
3.5. 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   Figure 9). The HOMO−LUMO energy gap value (E gap = E HOMO − E LUMO ) was 0.0884 eV. Molecules with a large E gap are generally stable or unreactive, while molecules with a small E gap are reactive. According to the HOMO− LUMO energy gap value of the kerogen model, the E HOMO and E LUMO 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 E HOMO and E LUMO 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.
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  10). 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 10) 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 10) 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 8), where the bond-order values of aromatic(C)−aromatic(C) are greater than those of alkane(C)−aromatic(C) in the kerogen model.

CONCLUSIONS
In this study, based on the experimental results obtained from solid-state 13 C 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 13 C NMR spectra revealed that the X BP 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 E angle . The multimolecular structure models of kerogen with different densities were established to evaluate the accuracy of

ACS Omega
http://pubs.acs.org/journal/acsodf Article 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 E gap 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.