Hydrodesulfurization of benzothiophene on Ni2P surface

The study of benzothiophene hydrodesulfurization reaction path contribute to clarifying the mechanism of hydrodesulfurization (HDS) of heavy oil. In this work, experiments and simulations were combined to study the reaction pathway of benzothiophene hydrodesulfurization catalyzed by Ni2P. In experimental part, Ni2P catalyst was prepared and characterized. Then, the catalytic property of the catalyst for benzothiophene hydrodesulfurization was evaluated. The substance types and contents in the liquid phase products were detected to verify the accuracy of the simulation results. Dmol3 module of the Materials Studio (MS) simulation software was used to simulate the adsorption and hydrodesulfurization of benzothiophene on the surface of Ni2P catalyst and explore the most probable reaction path. The results showed that the most stable adsorption configuration of benzothiophene on the surface of Ni2P was Ni-hcp. In addition, indirect desulfurization of benzothiophene was more advantageous than direct desulfurization. The most possible path for indirect desulfurization was Benzothiophene (BT) – Dihydrobenzothiophene (DHBT) – C8H9S2 – 2-phenylethyl mercaptan (PET) – Ethylbenzene (EB), while that of direct desulfurization was Benzothiophene (BT) – C8H7S2 – Styrene thiol (CMT) – Styrene (ST) – Ethylbenzene (EB).


Introduction
Utilization of the abundant heavy crude oil resources has been accepted as an effective approach to tackle the shortage of light oil supply in the coming decades (Montes et al., 2019). To solve the problems of high density, high viscosity, poor fluidity, and difficult exploitation of heavy oil, and to achieve the purpose of environmental protection, HDS is the most important and extensively used hydrotreating process (Peng et al., 2019). The conventional HDS is the catalytic removal of sulfur from sulfur-containing compounds in petroleum-based feed stocks, but compounds like thiophene, benzothiophene, dibenzothiophene and their derivatives were highly resistant to desulfurization due to their aromatic character and the weak donor ability of the sulfur atom (Willey et al., 1981). Among them, benzothiophene (BT) is a particularly interesting substrate since it is one of the most abundant compounds in heavy oils. Based on hydrogenation and desulfurization of BT, various mechanisms were proposed. Bartsch and Tanielian (1974) and Kilanowski and Gates (1980) observed that ethylbenzene (EB) was the only product of HDS of BT and suggested styrene (ST) to be the possible intermediate. De Beer et al. (1976) and Geneste et al. (1980) suggested that hydrogenation of BT should occur prior to desulfurization and the key intermediate should be dihydrobenzothiophene (DHBT). Daly (1978) and Van Parijs and Froment (1986) proposed that HDS of BT could have two reaction mechanisms of hydrogenation and hydrogenolysis.
To understand these processes, studies have been focused on the interaction scope and strength of sulfur-containing compounds with the active sides of catalyst surfaces leading to desulfurization. At present, most of the catalysts in industrial use were Co-Mo or Ni-Mo catalyst systems supported by alumina and silica. Some new nitride carbide and phosphide catalysts have excellent performance and good prospect, but they have not been applied to production practice. In transition metal phosphide catalysts, when the smaller non-metallic atoms P enter the metal lattice, they are iso-interstitially formed to form an octahedral interstitial compound. The non-metal atoms of this compound participate in the formation of covalent bonds, and the metal atoms in the compound have metal bonds, which make the compound not only have the electromagnetic properties of metal atoms, but also have good mechanical strength. Compared with transition metal nitrides and carbides with the same structure, transition metal phosphides show better hydrodesulfurization activity and better sulfur poisoning resistance. Compared with commercial catalysts, it has better hydrogenation selectivity and low oxygen consumption. Ni 2 P is an isotropic solid that exposes metal atoms on all surface planes, resulting in a higher density of active sites (Bussell, 2017;Jiang et al., 2019aJiang et al., , 2019bTopalian et al., 2019;Yu et al., 2019;Zhao et al., 2019;Zhu et al., 2019). Oyama et al.(2002) compared the HDS and hydrodenitrogenation (HDN) activities of different phosphides in model raw materials, and found that the order of activity of these phosphides is Fe 2 P <CoP <MoP <WP <Ni 2 P.
According to the investigation, it was found that benzothiophene is a representative sulfur-containing compound in petroleum and is of great significance for HDS research. As a new catalyst, Ni 2 P has higher activity than other traditional catalysts. Therefore, in this paper, Ni 2 P was selected as catalysts, and benzothiophene was used as a model compound for HDS experiments. MS software was used to simulate the HDS reaction of benzothiophene on Ni 2 P (001) surface. The various reaction mechanisms of the hydrodesulfurization process was theoretically analyzed, and the most possible reaction pathway was found.
In the experimental part, the catalyst was prepared and characterized. After the catalyst was successfully prepared, the catalyst was used for hydrodesulfurization of benzothiophene model oil, the gas and liquid products after the reaction were separated, and the composition and content of the product were detected. In the simulation part, the adsorption configuration was first determined, and then the reaction path was simulated based on the determined adsorption configuration, and the activation energy of different reaction paths was calculated through the transition state search. Finally, the most likely reaction path of hydrodesulfurization of benzothiophene was determined, and it was confirmed with the experimental results.

Experimental system and simulation method
Preparation of catalyst and hydrodesulfurization experiment Preparation of catalyst. 9.70 g of NiCl 2 Á6H 2 O(Sinopharm, AR) and 6.43 g of NaH 2 PO 2 ÁH 2 O (Aladdin, AR) were added to 50 mL of deionized water, and stirred at 293 K until completely dissolved, then the solution was dried at 363 K for 24 h. The dried sample was ground into powder having a particle size of 40-60 mesh. Then the powder was repeatedly washed and filtered, and dried to obtain Ni 2 P catalyst.
Catalyst characterization X-ray diffraction (XRD) analysis. The instrument used was an X'Pert Pro MPD model manufactured by Panalytical, the Netherlands. The measurement conditions were as follows: the diffraction source was Cu-Ka (k ¼ 1.54060), the tube pressure was 45 kV, the tube flow was 40 mA, the scanning speed was 5 /min, the divergence slit was 1 , and the receiving slit was 0.3 mm. The detector was a scintillation counter. The scanning angle was 5 -75 .
As shown in Figure 1, the diffraction peaks were wide, which proved that the prepared Ni 2 P was a crystal structure with small grains, and there were several obvious diffraction peaks at 2h ¼ 40. 7 , 44.6 , 47.3 , 54.2 , 54.9 , which was the same as the standard spectrum of Ni 2 P, indicating the successful preparation of the catalyst.
Catalytic evaluation. Before the hydrodesulfurization experiment, it was necessary to prepare a sulfur-containing model oil with a certain sulfur concentration. The benzothiophene-Noctane (Sinopharm,AR) model oil was selected in this study, and a certain amount of benzothiophene and N-octane were weighed. According to formula (1), 1000 ppm of benzothiophene model oil was prepared by gravimetric method, and the configured model oil was placed in a refrigerator for subsequent experiments. (1) Experimental steps: 1. Put 3.0g of Ni 2 P catalyst and 15mL of benzothiophene model oil into a Hastelloy reactor capable of withstanding high temperature and high pressure, and then put the reaction kettle into a tubular heating furnace. Drain the air in the reactor and check the air tightness of the device. Connect the gas cylinder and open the valve to continuously feed H 2 -N 2 mixed gas into the reactor until the pressure in the reactor reached 1.5MPa (H 2 content 3%, N 2 content 97%). The experimental device is shown in Figure 2. 2. Put the reaction kettle into a tube heating furnace, turn on the tube heating furnace and set the heating rate to 3K/min. A thermocouple for real-time temperature measurement was used to ensure that the reaction temperature was constant. The experimental conditions were selected as follows: reaction time was 5h, reaction pressure was 3.0 MPa, and reaction temperature was 573K. (Through pre-experiments, it was found that the reaction result was the best under this working condition)   3. After the reaction was completed, the tube heating furnace was turned off. When the reaction kettle was naturally cooled, the valve was opened to detect the H 2 S concentration (PN-2000, Shenzhen). During the test, the flow rate was kept stable and the detected gas was collected with an air bag. 4. The gas collected in the air bag was detected by an FID/TCD gas chromatograph analyzer (GC3800, United States). The collected liquid product was analyzed by using an FPD liquid chromatograph (LC, United States).

Hydroxydesulfurization MS simulation of benzothiophene
In this work, hydrodesulfurization of benzothiophene was theoretically explored using GGA-PW91 function based on the Density Functional Theory (DFT). In the calculation, Basis sets were set as DNP, and the electronic structure was analyzed by Density of States (DOS). At each step, the vibration frequency was calculated to ensure the stability of the adsorption configuration. In terms of core treatment, DSPP (DFT Semi-Core Pseudopots) was selected to approximate the electrons in the nucleus of metal atoms. For other atoms, such as hydrogen, carbon, nitrogen, and sulfur atoms, the All Electron Relativistic method was used for calculation. In the calculation process, Smearing was optimized from small to large by anneal function. All calculations used spin polarization. Considering the calculation efficiency and accuracy, the accuracy of all systems was set to Medium. The orbital occupation used the default 0.005 Hartree, and the Density Mixing was set to 0.17. The convergence tolerance of energy was 2 Â 10 À5 Hartree, the maximum gradient was set as 4 Â 10 À3 Hartree/Å , the maximum displacement was set as 5 Â 10 À3 Hartree, and the SCF was set as 1 Â 10 À5 Hartree. The method for DFT-D correction was used to calculate dispersion correction energy.
Determination of the adsorption model. The average bond lengths of Ni-P bond and Ni-Ni bond in nickel phosphide were 2.20 Å and 2.61 Å , respectively, and their unit cell parameters were, The benzene ring structure in the benzothiophene molecule was stable, and the thiophene structure was mainly considered in the simulation calculation. The C atoms around the S atom were named C 1 , C 2 , C 3 , and C 4 , respectively, as shown in Table 2. Since the unit cell structure of Ni 2 P was a close-packed form of ABAB, the distribution of Ni atom and P atom in the adjacent two layers were different. Therefore, there were two kinds of Ni atoms in the Ni 2 P single crystal cell, which were named Ni (1) and Ni (2). The Ni 2 P crystal had a hexagonal crystal structure. The surface density and the interplanar spacing of the (001) surface were the largest in all surfaces, and it was more likely to be broken and exposed on the surface of the catalyst. Therefore, the (001) surface was selected as the adsorption target. Since the Ni 2 P stacking mode was of the ABAB type, there were two (001) surfaces, Ni 3 P 2 surface (B) and Ni 3 P surface (C), and both cut surfaces were displayed in a 3 Â 3 periodic model, as shown in Table 2. The two (001) crystal faces obtained need to be structurally optimized again. Considering the change(i.e., surface relaxation) of surface atoms in the calculation process, the underlying atoms need to be fixed, so that the above three layers of atoms can move freely, as shown in the Table 2.
The surface energy of the two cut surfaces of Ni 3 P 2 and Ni 3 P were calculate. The smaller the surface energy, the more stable the corresponding surface was, then the optimal surface was obtained.
Surface energy is a measure of the destruction of intermolecular chemical bonds when the surface of a substance is created. It is also defined as the extra energy of surface particles relative to internal particles. Its definition is: In the formula, E slab represents the energy of the surface model, eV; E bulk represents the energy of the body model containing the same number of atoms, eV; S is the surface area, Å 2 ; E surf is the surface energy, eV.
After the surface of the Ni 3 P surface and the Ni 3 P 2 surface were relaxed, the structure was further optimized. The optimized lattice parameters were: a ¼ b ¼ 17.6721 Å and c ¼ 15.0979 Å . Table 1 shows the calculation results. The surface energy of the two crystal planes in the Ni 2 P (001) direction was not much different. In this work, the Ni 3 P 2 crystal plane in the Ni 2 P (001) direction was selected as the adsorption surface.
The adsorption sites of the adsorbed molecules on the determined surface can be classified into four types: a top site, a bridge site, and a triple vacancy (hcp site and fcc site). There were two ways to adsorb benzothiophene on the surface of Ni 2 P (001). One was parallel adsorption, that was, the benzothiophene molecule was parallel to the (001) surface; the other was vertical adsorption, that was , the benzothiophene molecule was perpendicular to (001) surface. The difference in adsorption mode determined the different adsorption energies, and the adsorption energies of various adsorption modes were calculated. The lower absolute value of the adsorption energy contributes to more stable adsorption configuration. The Ni 3 P 2 fourlayer surface model was chosen to construct the adsorption configuration. As shown in Table  2, the optimal adsorption configuration was screened by calculating five vertical adsorption configurations and five parallel adsorption configurations. It should be noted that each surface model was only allowed to adsorb one benzothiophene molecule during calculation.
The adsorption energy and structural parameters of benzothiophene on Ni 2 P (001) surface were calculated by simulation as shown in Table 3. In all adsorption configurations, the adsorption energy of hcp was higher than that of top and bridge, and Ni-hcp was the lowest (À1.09 eV). Therefore, the most stable adsorption configuration of benzothiophene on the surface of Ni 2 P was Ni-hcp.
Simulation of hydrodesulfurization reaction. During hydrodesulfurization, benzothiophene adsorbed H atoms in hydrogen. The position of the H atom can be divided into two types near the benzene ring and near the heterocyclic ring. However, considering that the benzene ring did not directly participate in the reaction in the hydrodesulfurization reaction, H atoms were selected to be adsorbed near the heterocyclic ring as the initial conditions in this work.
There were two reaction paths for benzothiophene hydrodesulfurization reaction. For the indirect hydrodesulfurization (HYD), the benzothiophene was first hydrogenated to form  dihydrobenzothiophene, and then the C-S bond was broken to remove the sulfur atom to form ethylbenzene. With regard to the direct hydrodesulfurization (DDS), the C-S bond was directly broken to form styrene, and then further hydrogenated to form ethylbenzene. Moreover, depending on the position of the C-S bond and the position of the hydrogenation, intermediates such as dihydrobenzothiophene and styrene can have different reaction paths, hence the research of specific reaction pathway for the hydrodesulfurization of benzothiophene on the surface of Ni 2 P (001) was very important. Through the transition state search, the activation energy of each reaction pathway was calculated, and the transition state was searched by LST/QST tool in MS software to simulate each possible path in the benzothiophene hydrodesulfurization reaction.
In the initial stage of HYD, benzothiophene was first subjected to hydrogenation saturation of C ¼ C bond to form dihydrobenzothiophene, but in the C-S bond cleavage stage, two different intermediates would appear depending on the position of the hydrogenolysis reaction. Then different paths were produced depending on whether the two intermediates were hydrogenated, and the specific reaction mechanism for HYD was represented by Figure 3(a) to (e) was used to name each reaction, and the four possible reaction mechanisms for HYD of benzothiophene were listed in Table 4. The reactants and products of all reaction were structurally optimized, and the most stable parallel adsorption configuration hcp was selected as the initial adsorption form of the reactants and products. Hydrogen atom was selected as the hydrogen source in the reaction. The transition state search of each reaction was performed to calculate the energy of the reaction.
For DDS, benzothiophene first underwent C-S bond cleavage, and there were two routes depending on the position of the fracture. After desulfurization, styrene was formed and finally hydrogenated to form ethylbenzene, as shown in Figure 4. The reactants and products of all reaction mechanisms in Table 5 were structurally optimized, and the most stable adsorption configuration was selected as the initial adsorption form of the reactants and products, and the transition state search was performed to calculate the reaction energy. By calculating the reaction energy, the reaction path with the lowest activation energy was the most likely reaction path.

Benzothiophene hydrodesulfurization reaction
In order to investigate the possible route of the hydrodesulfurization reaction of benzothiophene, the gas product and the liquid product were examined, and the results are shown in Tables 6 and 7. The gas product was analyzed by gas chromatography. As shown in Table 6, the gas products were mainly hydrogen sulfide, methane, ethane, and ethylene. Since the carrier gas was a H 2 /N 2 mixed gas, the two gases in the product account for the majority. Since the structure of the benzene ring was stable and difficult to be destroyed, the butene and butane in the product may be formed by a polyaddition reaction of a monomer such as methane, ethane and ethylene, or by C ¼ C bond cleavage decomposition of N-octane.
The detection results of the organic products are shown in Table 7, from which the unreacted N-octane and benzothiophene were removed. The analysis believed that the heptane in the product was produced by the dehydrogenation of N-octane. After part of the benzothiophene carbon chain fallen off, benzene and toluene were produced, and with the Table 4. Mechanisms for HYD of benzothiophene(A, B, C, D). Figure 4. DDS reaction pathway of benzothiophene (Wang and Prins, 2008;Yao et al., 2005). Figure 3. HYD reaction pathway of benzothiophene (Wang and Prins, 2008;Yao et al., 2005).

Reaction step Reaction mechanism A Reaction mechanism B Reaction mechanism C Reaction
reaction continued, naphthalene, biphenyl, xylene, ethylbenzene, 1-ethyl and 3-methylbenzene were produced. Part of benzothiophene undergone hydrogenation reaction to generate 3,6-dimethylbenzothiophene. Part of benzothiophene desulfurization produced cyclopentenebenzene, 2-phenylethyl mercaptan, and styrene mercaptan. Among them, 2-phenylethyl mercaptan and styrene mercaptan were intermediate products in the optimal reaction path of HYD and DDS obtained by simulation, respectively. The content of the two intermediate products (2-phenylethyl mercaptan: 9.4%, Styrene mercaptan: 4.0%) proved that HYD was  Table 5. Reaction mechanisms for direct hydrogenation of benzothiophene(E, F, G).

Reaction step
Reaction mechanism E Reaction mechanism F Reaction mechanism G 1 a a b 2 a 1 a 2 c 3 e f d 4 f f more dominant than DDS, which was the main hydrodesulfurization path of benzothiophene. After the reaction, most of the sulfur in benzothiophene was released through the desulfurization reaction to form hydrogen sulfide. Some remained in the intermediate product of the reaction, and a small amount of sulfur remained in 3,6-dimethylbenzothiophene by hydrogenation of the benzothiophene benzene ring.

Hydrogenation reaction path simulation
Indirect hydrogenation. It was found through simulation that the most likely reaction mechanism for HYD is D, and the most likely reaction mechanism for direct hydrogenation is G.
Only the above two reaction paths are described here, and the remaining paths are detailed in Appendix A.
In reaction path D shown in Figure 5, benzothiophene was first hydrogenated to generate dihydrobenzothiophene. The C 4 -S bond in dihydrobenzothiophene was broken to form the intermediate product C 8 H 9 S 2 , and then hydrogenation was continued to generate 2-phenylethyl mercaptan. Finally, sulfur atoms were removed to form ethylbenzene. Figure 6 shows the reaction path D in the HYD process. In reaction a, the H atom, which was originally in stable adsorption, moved closer from the hcp position to the benzothiophene, and the benzothiophene molecule also had a certain displacement. And the original H atom on C 1 atom has also shifted to some extent. In the transition state structure, the original C ¼ C bond became a single bond, the C 1 atom was combined with an external H atom, C 1 and S atoms were slightly shifted. After hydrogenation, the dihydrobenzothiophene molecules rotated to a certain degree. The initial state (IS) of the reaction was benzothiophene, and the final state (FS) of the reaction was dihydrobenzothiophene. The transition state (TS) was selected and calculated using the LST/QST tool.
Reaction c represented the formation process of the intermediate product C 8 H 9 S 2 . During the reaction, H atoms moved closer to C 4 atoms, while S atoms moved closer to Ni atoms. Under the joint action of the two, the C 4 -S bond was broken, and the H atom was finally connected with the C 4 atom to form C 8 H 9 S 2 .
The Figure 6 (reactions c 1 and e) shows the two processes of C 8 H 9 S 2 hydrodesulfurization to ethylbenzene. First, in reaction c 1 , the H atom gradually approached the S atom on C 8 H 9 S 2 , and the two combined to form 2-phenethylthiol. When observing the distance between the S atom and the Ni atom, it was found that the distance increased, and at the same time the entire molecule appeared to a certain degree of rotation. Then in reaction e, the H atom gradually moved closer to the C 1 atom on the 2-phenethyl mercaptan, and the S atom began to approach the Ni atom, which brought the molecule closer to the Ni 2 P (001) surface and extended the C 1 -S bond. Eventually, the -HS structure fell off, adsorbed on the Ni atom, and continued to be hydrogenated to H 2 S, and the H atom combined with C 1 to form ethylbenzene. Table 8 shows the activation energy Ea and the change in reaction energy DE of each reaction obtained from the transition state search. The heat absorption was positive, and the heat release was negative. The most possible path for HYD of benzothiophene can be acquired by comparing the activation energies of the different reaction pathways in Table 8. In HYD process, it can be divided into four reaction paths of A, B, C, and D. Process A included four reactions of a, b, b 1 , and d, and the highest activation energy was 1.98 eV. Process B included three reactions a, b, and b 2 , and the highest activation energy in route B was 1.98 eV. Process C included three reactions, a, c, and c 2 , and the highest activation energy in route C was 1.33 eV. Process D included four reactions a, c, c 1 , and e, the highest activation energy in route D was 1.24 eV. Comparing the highest activation energy of the four reaction paths, the most possible path for HYD of benzothiophene was: BT-DHBT-C 8 H 9 S 2 -PET-EB (reaction mechanism D).
Direct hydrogenation. In the reaction path G shown in Figure 7, the C 4 -S bond in benzothiophene was first broken to form C 8 H 7 S 2 , and then hydrogenated to generate styrene mercaptan, then -HS was dropped to generate styrene, and finally styrene was hydrogenated to generate ethylbenzene.  Reaction b in Figure 8 shows the formation process of the intermediate product C 8 H 7 S 2 . During the reaction, the H atom was closer to the C 4 atom, and the S atom was closer to the Ni atom. Under the combined action of the two, the C 4 -S bond was broken, the H atom and the C 4 atom were connected to form C 8 H 7 S 2 . After the C 4 -S bond was broken, the benzene ring lost the S atom and moved away from the surface of Ni 2 P (001). Reaction c shown in Figure 8 shows the hydrogenation of C 8 H 7 S 2 to styrene. The H atom gradually approached the S atom, and the two combined to form a styrene thiol finally. At the same time, the distance between the S atom and the Ni atom became longer, and the molecule was far away Table 9. Activation energy and reaction energy of each reaction on Ni 2 P (001) plane.

Reaction
Ea  from the Ni 2 P (001) surface. In reaction d, the H atom gradually moved closer to the C 1 atom on the styrene thiol molecule, while the S atom approached the Ni atom, which made the C 1 -S bond stretch, and eventually the -HS structure fell off and was adsorbed on the Ni atom. The H atom combined with C 1 to form styrene. The last step was the hydrogenation of styrene to ethylbenzene, which was reaction f, as shown in Figure 8. Table 9 shows the activation energy Ea and the change in reaction energy DE of each reaction obtained from the transition state search. In the DDS process, it can be divided into three reaction paths: E, F, and G. Process E included four reactions, a, a 1 , e, and f, and the highest activation energy was 2.29 eV. Process F included three reactions a, a 2 , and f, and the highest activation energy in route F was 2.29 eV. Process G included four reactions b, c, d, and f. The highest activation energy in route G was 1.45 eV. Comparing the highest activation energies of the three reaction paths, it was more inclined to reaction mechanism G, and the most possible path for DDS of benzothiophenes was: BT-C 8 H 7 S 2 -CMT-ST-EB (reaction mechanism G).
As shown in Figure 9, the highest activation energy barrier of HYD reaction was smaller than that of DDS (Ea ¼ 1.24 eV<Ea ¼ 1.86 eV), indicating that HYD was more advantageous in the hydrodesulfurization reaction of benzothiophene. Benzothiophene was more likely to converted to dibenzothiophene as an intermediate product and then toward diphenylethanethiol in the reaction. In the HYD, the activation energy barrier of TS2 was the highest, indicating that the rate limiting step in the reaction path was the process of DHBT breaking C 4 -S bond to form C 8 H 9 S 2 . In general, both the DDS reaction and the HYD reaction were exothermic reactions. Figure 10 shows the most likely hydrodesulfurization reaction path of benzothiophene on the Ni 2 P (001) surface.

Conclusion
The specific reaction pathway of benzothiophene hydrodesulfurization is found by experiment and simulation. The surface of Ni 3 P 2 in the (001) plane of the Ni 2 P crystal is as stable as the surface of Ni 3 P within the error range. When benzothiophene is adsorbed in parallel on the Ni-hcp position of Ni 3 P 2 (001) plane, the adsorption energy is lowest and the adsorption is most stable. The HYD of benzothiophene is more advantageous than DDS. The most possible path for HYD is BT-DHBT-C 8 H 9 S 2 -PET-EB. The most possible path for DDS is BT-C 8 H 7 S 2 -CMT-ST-EB. Both of the two processes are exothermic reactions. The accuracy of the simulation results is verified by detecting the substance types and contents in the liquid phase products. 2-phenylethyl mercaptan and styrene mercaptan were intermediate products in the optimal reaction path of HYD and DDS obtained by simulation, respectively. The content of the two intermediate products (2-phenylethyl mercaptan: 9.4%, Styrene mercaptan: 4.0%) proved that HYD was more dominant than DDS, which was the main hydrodesulfurization path of benzothiophene.