Kinetic Modeling of Hydrotreating for Enhanced ... - ACS Publications

Jun 19, 2019 - has been done by means of a statistical significance test based on Fisher's method. ...... concentration of polyaromatics with temperat...
0 downloads 0 Views 1MB Size
Article Cite This: Ind. Eng. Chem. Res. 2019, 58, 13064−13075

pubs.acs.org/IECR

Kinetic Modeling of Hydrotreating for Enhanced Upgrading of Light Cycle Oil Roberto Palos,†,§ Alazne Gutieŕ rez,*,† Idoia Hita,†,⊥ Pedro Castaño,†,⊥ Joris W. Thybaut,‡ Jose ́ M. Arandes,† and Javier Bilbao† †

Department of Chemical Engineering, University of the Basque Country UPV/EHU, P.O. Box 644, 48080 Bilbao, Spain Laboratory for Chemical Technology, Ghent University, Technologiepark 125, B-9052 Ghent, Belgium



Downloaded via BUFFALO STATE on July 31, 2019 at 10:07:22 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.

S Supporting Information *

ABSTRACT: The hydrotreating of light cycle oil (LCO) into high-quality fuels has been investigated experimentally and kinetically, developing a model that accounts for the main and simultaneous reaction pathways: hydrodesulfurization (HDS), hydrodearomatization (HDA), and hydrocracking (HC). The experiments have been carried out in a fixed-bed reactor, NiMo/SiO2−Al2O3 commercial catalyst, 320−400 °C; 80 bar; space time, 0−0.5 gcat h gLCO−1; and H2/LCO volumetric ratio of 1000 Ncm3 cm−3. The proposed kinetic model contains multiple lumps, species, and pathways, leading to the faithful prediction of hydrotreatment products from different viewpoints. The computed kinetic parameters have allowed for simulating the process and seeking the optimal operating conditions. This way, the maximum values obtained for the conversions of HDS, HDA, and HC have been of 90%, 20%, and 65%, respectively; whereas a good compromise between the different hydrotreating goals has been obtained in the 385−400 °C range for a space time of 0.2 gcat h gLCO−1. Finally, the obtained optimal operating conditions have been compared with those optimized in the literature.

1. INTRODUCTION Oil refiners have been facing continuous changes in light of economic crisis, changes in policies, environmental concern, or wars, among others.1 Along these lines, refiners must update their activities to reduce wastes and the environmental impact. Quality standards of automotive fuels have been continuously adjusted to technological and environmental requirements in terms of the content of heteroatoms and of aromatics. Nowadays, the concentration of these undesired substances is limited in the EU to 10 ppm of sulfur and to 35 vol % of aromatics.2 Furthermore, taking into account that extracted oil is increasingly heavier and sourer and that only one-third of the remaining petroleum reserves can be considered as conventional crude oil,3 the upgrading of secondary refinery streams is becoming a key factor for reducing wastes. Thus, among the different refinery streams, light cycle oil (LCO) can be considered as an interesting source of gasoline and middle distillates.4,5 However, its high content of aromatics (up to 88 wt %) and sulfur (up to 4 wt %) appear as their main drawbacks to be directly used in the blending of automotive fuels.6 Therefore, the most suitable strategy for the valorization of LCO and other secondary streams consists of their hydrotreating in two consecutive stages: the first stage using transition metals-based catalysts primarily for the removal of heteroatoms from the LCO and the second stage with noble metal-based catalysts to adapt the composition of the hydrotreated LCO to comply with environmental regulations.7 © 2019 American Chemical Society

Kinetic models that describe the chemical reactions that take place during hydrotreating are essential tools for a rational reactor design, control, and optimization of the commercial plants.8 Therefore, the economic and strategic relevance of hydrotreating processes justifies the effort done to establish kinetic models useful to determine both catalyst and reactor design that conditions the economic viability of the industrial units. Several works in the literature have investigated the kinetic modeling of the secondary streams hydrotreatment.9−12 However, the hydrotreating of real and complex refinery feeds is remarkably different due to the large number of reactions and the complicated reaction network. Hence, modeling reaction kinetics of hydroprocessing units faces multiple challenges because carrying a detailed component identification is a complex task. Consequently, the complexity and heterogeneity of the commonly hydrotreated feedstocks have driven the attention of kinetic models toward model compounds.13−15 On the other hand, many models found in the literature10,16,17 gather the thousands of components present in a typical petroleum fraction into a tractable number of lumps with a simplified reaction network. Therefore, establishing a simplified reaction network that accurately describes the process while keeping a proper balance between a Received: Revised: Accepted: Published: 13064

April 17, 2019 June 11, 2019 June 19, 2019 June 19, 2019 DOI: 10.1021/acs.iecr.9b02095 Ind. Eng. Chem. Res. 2019, 58, 13064−13075

Article

Industrial & Engineering Chemistry Research

(ICP-AES). The textural properties of the catalyst have been determined by means of N2 adsorption−desorption isotherms obtaining a specific surface area of 208 m2 g−1 and a mesoporous volume of 1.90 × 10−2 m3 g−1. In order to determine its acidic properties, tert-butylamine (t-BA) adsorption followed by a temperature-programmed desorption (TPD) has been performed, obtaining a total acidity of 208 mmolt‑BA g−1, as well as infrared spectroscopy measurements using pyridine as a probe molecule to obtain the Brønsted/ Lewis acidic site ratio of 2.25. A more extended and detailed explanation of the employed equipment and methods for determining the properties of the catalyst can be found elsewhere.19 2.2. Hydrotreating of LCO. The runs have been performed in a three-phase laboratory scale fixed bed reactor working in the trickle-bed regime. The schematic representation of the setup and the operating procedure of employed unit can be found elsewhere.7 Employed reaction conditions have been: temperature, 320−400 °C; pressure, 80 bar; space time (τ), 0−0.5 gcat h gLCO−1; H2/LCO volumetric ratio, 1000 Ncm3 cm−3; and, time on stream (TOS), 0−8 h. Prior to reactions, the catalyst (diameter between 150 and 300 μm) has been sulfided in situ for 4 h at 400 °C under a flow rate of 50 cm3 min−1 of a commercial H2S/H2 mixture (Air Liquide, 10% vol H2S). Results have been studied for a TOS of 8 h, after checking that for TOS > 6 h a pseudostable state with constant activity is reached, as a consequence of attaining a stable trickle-bed regime and stable sulfurization state of the metallic function.20 The products distribution, as well as the distillation range distribution of the liquid products, is collected in Table S1 in the Supporting Information. Different reaction indexes have been defined for measuring the extent of hydrodesulfurization (HDS), hydrodearomatization (HDA), and hydrocracking (HC) reactions, based on the concentration of sulfur, total aromatics mass fraction, and mass fraction of gasoil, respectively. Those indexes are

detailed or a simplified kinetic approach presents many advantages. Some advantages are (i) shorter calculation and computation times, (ii) extra reactions that can be easily added, and (iii) a kinetic model that can be easily implemented in process simulation software. In this context, Hita et al.18 made a notorious advance to make kinetic models more closely representing reality by separately considering the kinetics of the three key reactions in the hydrotreating of scrap tires pyrolysis oil. In the present work, we have studied experimentally the first stage of the dual hydrotreating strategy of LCO in a fixed-bed reactor using a NiMo/SiO2−Al2O3 commercial catalyst. Furthermore, with the aim of proposing a kinetic model based on multiple lumps, species, and pathways, a wide range of operating conditions has been tested. Our goal has been to individually quantify the extent of the main reactions: hydrodesulfurization (HDS), hydrodearomatization (HDA), and hydrocracking (HC). The kinetic model reproduces the product composition in terms of sulfur content, chemical nature (paraffins and isoparaffins, naphthenes, 1-ring and 2+ring aromatics), and boiling point range (naphtha, diesel, and gasoil).

2. EXPERIMENTAL SECTION 2.1. Materials. Light Cycle Oil (LCO) is a byproduct of the FCC unit of Petronor S.A. refinery located in Muskiz (Biscay, Spain). The main physicochemical properties of LCO are shown in Table 1. LCO is a quite light (API gravity: 27.7°) Table 1. Main Physicochemical Properties of Light Cycle Oil (LCO) property gravity (°API) density, at 15 °C (kg m−3) kinematic viscosity, at 40 °C (cSt) elemental analysis (wt %) carbon hydrogen sulfur nitrogen simulated distillation (°C) IBP−FBP T50−T95 distillation fractions (wt %) naphtha (350 °C) chemical composition (wt %) paraffins and isoparaffins (P+iP) naphthenes (N) 1-ring aromatics (A1) 2+-ring aromatics (A2+)

value 27.7 887 2.71 89.56 9.37 0.102 0.067

Hydrodesulfurization conversion

XHDS =

(S)LCO − (S)Prod (S)LCO

(1)

where (S)LCO and (S)Prod are the amounts of sulfur (in ppm) in the LCO and in the liquid product stream, respectively

95.7−438.0 288.6−381.3

Hydrodearomatization conversion

XHDA =

(xArom)LCO − (xArom)Prod (xArom)LCO

(2) 11.0 69.7 19.3

where (xArom)LCO and (xArom)Prod are the mass fractions of total aromatics in the LCO and in the liquid product stream, respectively

34.25 3.65 31.10 31.0

Hydrocracking conversion

XHC =

(xGasoil)LCO − (xGasoil)Prod (xGasoil)LCO

(3)

where (xGasoil)LCO and (xGasoil)Prod are the mass fractions of gasoil (TB > 350 °C) in the LCO and in the liquid product stream, respectively. 2.3. Data Treatment. The feedstock and reaction products have been analyzed by means of comprehensive two-dimensional gas chromatography with mass selective detection spectrometry (GC×GC-MSD), enabling an easier identification of the chemical nature of the individual components and, hence, the entire stream. The GC×GC-MSD apparatus comprises an Agilent Technologies 5975C inert XL mass spectrometer coupled in line with an Agilent Technologies 7890A gas chromatograph equipped with two capillary

and high sulfur content (10212 ppm) secondary refinery stream having a diesel content (boiling point, TB 489−623 K) of 69.7 wt %. Besides, it is mainly composed by aromatics (62.1 wt %) with an almost equal distribution of mono- and polyaromatics (31.1 and 31.0 wt %, respectively). A commercial NiMo/SiO2−Al2O3 catalyst supplied by Akzo Nobel has been used for the LCO hydroprocessing. This catalyst contains 2.98 wt % Ni and 7.31 wt % Mo, measured by inductively coupled plasma atomic emission spectroscopy 13065

DOI: 10.1021/acs.iecr.9b02095 Ind. Eng. Chem. Res. 2019, 58, 13064−13075

Article

Industrial & Engineering Chemistry Research

that the improvement of model B in comparison with model A is significant when the following condition is fulfilled25

columns connected serially by means of a valve modulator. The first column is a nonpolar DB-5 ms J&W 122-5532 (30 m × 0.25 mm), whereas the second one is a polar HPINNOWAX (5 m × 0.25 mm). The sulfur content has been measured by a GC system (Agilent Technologies 7890A GC System) equipped with a specific sulfur detector (PFPD) and a HP-PONA (50 m × 0.20 mm) capillary column. Simulated distillation analyses have been performed on a GC system (Agilent Technologies 6890 GC System) according to the ASTM D-2887 Standard and using a DB-2887 (10 m × 0.53 mm) semicapillary column. 2.4. Modeling Methodology. Considering the systematic methodology for chemical kinetics modeling established by Toch et al.,21 the proposed reaction schemes for hydrodesulfurization, hydrodearomatization, and hydrocracking reactions are composed of j reactions in which i lumps are involved. Therefore, the kinetic equations of the formation of each i lump are defined as ri =

FA − B =

i xi* − xi yz zz zz * x i̅ k {

∑ jjjjj nl

i=1

(6)

SSE = (xi* − xi)2

(7)

ν = nlne − n p

(8)

where nl, ne, and np are the number of lumps, experiments, and parameters of the model, respectively. The test based on the Fisher’s method has been applied to distinguish the different reaction schemes proposed for the hydrodearomatization route (HDA). Additionally, the use of a sensitivity analysis to estimate the best set of kinetic parameters is of great interest in heterogeneous kinetic models.22 This strategy consists of assessing how perturbations of the studied kinetic parameters, while maintaining constant the remaining ones, affect the SSE function.

(4)

where ri and rj are the formation rate of lump i and and rate of reaction j, respectively, and (υi)j is the stoichiometric coefficient of the lump i in each j reaction. The system of kinetic equations considered for each reaction scheme has been approached by developing a MATLAB routine. First, it reads the experimental data and the initial values assigned to the apparent kinetic parameters and apparent activation energies. According to the criteria proposed by Alcázar and Ancheyta,22 there have been chosen as initial values for the apparent kinetic parameters and activation energies those previously reported in the literature.18 Afterward, the routine determines the numeric values of the kinetic parameters that minimize the sum of errors in regression between the predicted and the experimental data. Thus, the parameters that minimize the following objective function are calculated OF =

> F1 − α(νA − νB , νB)

where F1−α (vA − vB,vB) is the critical value for the Fisher’s distribution for a level of significance of 95% (α = 0.05). If the condition is not fulfilled, the improvement is not statistically significant. Therefore, if the compared models have different number of parameters, the model with the lower number of parameters is preferred. To carry out the statistical test, the sum of squared errors (SSE) and degrees of freedom (ν) of the kinetic models are defined as follows

∑ (υi)j rj j

SSEA − SSE B SSE B νA − νB νB

3. RESULTS AND DISCUSSION 3.1. Kinetic Models. Based on the established sulfur species, product lumps, and composition fractions, a lumpbased kinetic model that individually quantifies the extent of HDS, HDA, and HC pathways, considering a specific kinetic expression for each route, has been proposed. Accounting for these three different goals simultaneously is a useful strategy to properly simulate the behavior of the reactor. 3.1.1. Hydrodesulfurization (HDS). A proposed reaction scheme for hydrodesulfurization consists of a set of independent reactions, each of them accounting for the conversion of each sulfur-containing lump to obtain H2S and different sulfur-free hydrocarbons as products. The lumps of sulfur-containing species that have been considered in the kinetic model have been the following: benzothiophene (BT), methyl-substituted benzothiophenes (MxBTs), dibenzothiophene (DBT), and methyl-substituted dibenzothiophenes (MxDBTs). More detailed information about the HDS reactions of the aforementioned sulfur-containing species can be found elsewhere.26 A Langmuir−Hinshelwood type of mechanism, postulating that the surface reaction takes place between two adjacently adsorbed reactants, is the most common form of expression for HDS kinetics.27−29 Moreover, the H2S produced from reactive sulfur compounds in the early stage of the reaction is one of the main inhibitors for HDS of the unreactive species.30 It should be considered in kinetic equations by including a competitive adsorption term (KH2S). In order to establish the kinetic equations, the assumptions were as follows: (i) HDS reactions are irreversible; (ii) there is no interaction between the different sulfur species; (iii) due to the large H2 excess, its concentration is effectively constant

(5)

where xi* is the mass fraction of lump i, xi is the corresponding predicted value, and x*i is the arithmetic mean of experimental mass fractions of lump i. It should be mentioned that in hydrodesulfurization the values of xi* and xi correspond to sulfur concentrations (in ppm) in each family (lump) of sulfur compounds. Furthermore, it is required to determine the concentration of the dissolved H2 in LCO, which depends on fugacity, the activity coefficient, and H2 solubility in the LCO. These parameters have been calculated following the empirical correlations described by Riazi and Roomi,23 taking into account the aromatic nature of the LCO, that conditions the solubility of the H2.24 2.5. Test of Significance. Once kinetic parameters have been obtained and in the case that more than one kinetic scheme is proposed, discrimination between different models has been done by means of a statistical significance test based on Fisher’s method. As proposed kinetic models have different degrees of freedom, Fisher’s exact test for the comparison of two models A and B (νA ≠ νB), when model B shows a smaller sum of squared errors than model A (SSEB < SSEA), considers 13066

DOI: 10.1021/acs.iecr.9b02095 Ind. Eng. Chem. Res. 2019, 58, 13064−13075

Article

Industrial & Engineering Chemistry Research relative to the concentration of any of the sulfur-containing lumps; therefore, the kinetic equation can be collapsed into a pseudo-first-order one; and (iv) the equilibrium constant for H2S adsorption is considered independent of reaction temperature. With these premises, the rate of disappearance of each sulfur lump has been defined as −

kjxix H2 k′jxi dxi = = 2 dτ (1 + K H2Sx H2S) (1 + K H2Sx H2S)2

(9)

gLCO−1),

where τ is the space time (gcat h kj is the intrinsic kinetic parameter of the j reaction, k′j is the apparent kinetic parameter of the j reaction, xi is the sulfur concentration (ppm) of the i lump, and xH2 and xH2S are the mass fractions of H2 and H2S in the reaction medium, respectively. It should be mentioned that kinetic modeling of the reactions of HDS, HDA, and HC has been done considering the concentration of H2 (xH2) constant, hence computed kinetic parameters for the three different routes have been the apparent ones (k′j). Even though this strategy is commonly followed in kinetic studies of hydroprocessing, determining the consumption of hydrogen, by means of experimental data either by a hydrogen balance in gas streams or with hydrogen content in the liquid feed and products, is crucial as it will determine the economy of the industrial process.31 On the other hand, the power 2 in the denominator in eq 9 is in concordance with general principles of the LH mechanism, i.e., revealing that two active sites are involved in the rate-determining step.32,33 Additionally, in order to avoid too pronounced correlation between the apparent kinetic parameter and the apparent activation energy during regression, the kinetic parameter has been expressed as a function of temperature according to the reparameterized Arrhenius equation ÅÄÅ −E i ÑÉ ÅÅ j j 1 1 yzzÑÑÑÑ j Å ′ ′ kj = kj , Tr expÅÅ j − zzzÑÑ ÅÅ R jj T Tr {ÑÑÑÖ (10) k ÅÇ

Figure 1. Proposed reaction schemes for HDA reactions of LCO.

reactions are of first order with respect to the products.37 Hence, the following kinetic equations are obtained (referred to mass fractions of each lump in the reaction scheme) for the most complex reaction scheme, i.e., scheme C dx A 2 + dτ

ij x A1x H2 yzz zz − k5x A x H = −k1jjjjx A 2+x H2 − 2+ 2 j Keq,A 2+ − A1 zz k {

(11)

ij i x A1x H2 yzz x x y zz − k 2jjjjx A x H − N H2 zzzz − k4x A x H = k1jjjjx A2 +x H2 − z j 1 2 1 2 z j j Keq,A2 +− A1 Keq,A1− N zz dτ { k { k

dx A1

where k′j,Tr is the apparent kinetic constant of the j reaction at the reference temperature, Ej is the apparent activation energy of the j reaction, R is the universal gas constant, T is the temperature, and Tr is the reference temperature. 3.1.2. Hydrodearomatization (HDA). In HDA reactions the relation between the lumps corresponds to a complex system of individual reactions as it has been previously reported in literature.3,34−36 Hence, in order to describe the HDA reactions accurately, three different reaction schemes (Figure 1) involving four lumps have been considered. Scheme A, the simplest one, considers that different lumps are formed by reversible reactions in series, meaning that 2+- and 1-ring aromatics (A2+ and A1, respectively) are sequentially hydrogenated on the metallic sites of the catalysts, leading to naphthenes (N) which are, in turn, converted in linear and/or branched paraffins (P+iP) via ring opening reactions. In the following schemes possible alkyl cleavages from 1-ring (schemes B and C) and 2+-ring aromatics (scheme C) leading to the formation of paraffins are added. In view of the above, the kinetic model for HDA reactions has been proposed assuming that (i) each of the aromatic lumps acts as reactant and product in reversible hydrogenation reactions; (ii) due to the large excess of H2 used in the experiments, hydrogenation reactions are of first order with respect to the aromatic concentration; and (iii) reverse

ij ij x Nx H2 yzz x x yz dx N zz − k 3jjjx Nx H − P + iP H2 zzz = k 2jjjjx A1x H2 − z j 2 j j Keq,A1− N z dτ Keq,N − P + iP zz k { k {

(12)

x P + iPx H2 zyz dx P + iP ji zz + k4 x x H + k5x A x H = k 3jjjjx Nx H2 − 2+ 2 A1 2 j dτ Keq,N − P + iP zz k {

(13)

(14) −1

where τ is the space time (gcat h gLCO ), kj is the intrinsic kinetic parameter for the j reaction, xi is the mass fraction of the i lump, and xH2 is the mass fraction of H2 in the reaction medium, respectively. Note that in this case i could be P+iP, N, A1, or A2+ which correspond to paraffins and isoparaffins, naphthenes, 1-ring aromatics, and 2+-ring aromatics, respectively. It should be mentioned that, to obtain the equations that describe scheme A, kinetic parameters k4 and k5 in eqs 11, 12, and 14 must be equal to zero. On the other hand, in order to obtain the equations of scheme B, just the kinetic parameter k5 in eqs 11 and 14 must be equal to zero. Once again, the temperature dependence of the kinetic parameters has been expressed according to the reparameterized Arrhenius equation (eq 10). Due to the complexity of the feedstock and the lumps considered in this approach, equilibrium parameters have not been determined from thermodynamic properties of the components but from the ratio of the experimental data at the highest value of space time (τ = 0.5 gcat h gLCO−1) as an apparent equilibrium state is reached at space times higher 13067

DOI: 10.1021/acs.iecr.9b02095 Ind. Eng. Chem. Res. 2019, 58, 13064−13075

Article

Industrial & Engineering Chemistry Research than 0.2 gcat h gLCO−1. This phenomenon has been already observed before18,38 and can be attributed to the reversibility of implied individual hydrogenation reactions. Then, equilibrium parameters for each j stage of the kinetic scheme have been defined as the ratio of the concentration between the products and the reactants involved in each reaction Keq, j =

k′j k −′ j

=

[Prod]i [React]i

distillation fractions. The interconversion takes place through first-order reactions, in such a way that the components of a lump are transformed into other components of a lighter lump. The proposed kinetic scheme is based on the following statements: (i) as it is widely known,43,44 commercial hydrotreating catalysts take on average 4−6 years to be replaced, a time that can be shorter (1−2 years) if heavy oils are hydrotreated and longer if naphtha or diesel are hydrotreated (up to 8 years). Hence, given the composition of the LCO chosen in this work (Table 1), it can be assumed that the catalyst used (commercial NiMo/SiO2Al2O3) could last ca. 5 years in an industrial unit. On the other hand, various strategies can be found in the literature to accelerate the precoking of the hydrotreating catalysts, which consist of increasing LHSV, temperature elevation, use of feeds with high metal content, or a combination of these factors.45−48 However, none of the accelerated aging strategies are capable of fully representing practical deactivation. Therefore, as longterm coke deposition does not influence the reaction scheme and the aim of this work is to establish a scheme capable of reproducing the hydrocracking route, the inclusion of a coke lump in the model is not required. Additionally, to ratify this statement a test of 100 h long has been performed, in which no deactivation of the catalyst has been observed during this period; (ii) a certain degree of conversion from gasoil to naphtha can be assumed; and (iii) it has been verified (Table S1) that the formation of lighter compounds than naphtha (liquefied petroleum gases) is minimum for tested catalyst (designed to avoid overcracking) and within the range of investigated experimental conditions (working at temperatures below 400 °C), it is not necessary to include a lump of these species. In view of the above, the kinetic model for HC has been proposed assuming that (i) reaction stages could be reversible; (ii) hydrocracking reactions are of first order with respect to reactant lump due to the large excess of H2 used in the experiments; and (iii) the activity of the catalyst remains constant, meaning that no catalyst deactivation is observed. Therefore, the equations that describe the evolution with space time of the concentration of each lump are

(15)

where k′j and k′−j are the direct and reverse apparent kinetic parameters, respectively, of reaction j. Then, based on the prior definition of the equilibrium parameter (eq 15), the following temperature dependent empiric correlations have been obtained for the different reversible reactions Keq,A 2+ − A1 = 1.798·10−3T 2 − 2.157T + 6.494· 10−2

(16)

Keq,A1− N = −2.414·10−6T 2 + 3.132· 10−3T − 8.885· 10−1 (17) −4 2

−1

1

Keq,N − P + iP = 1.962·10 T − 2.563· 10 T + 9.129· 10

(18)

where Keq,A2+−A1 corresponds to the equilibrium between lumps A2+ and A1, Keq,A1−N corresponds to the equilibrium between lumps A1 and N, and Keq,N−P+iP corresponds to the equilibrium between lumps N and P+iP. For the gas-phase hydrogenation from benzene to cyclohexane at gas phase conditions, Carrero-Mantilla and LlanoRestrepo39 have obtained an empiric correlation to determine the equilibrium parameter, obtaining similar values to those calculated by eq 17. As expected, the equilibrium parameter Keq,A1−N reaches its maximum value at 376 °C, indicating that at higher temperatures, as the reaction is highly exothermic, the reverse reaction is favored. This thermodynamic limitation is well established in the literature, and it has been observed by several authors at ca. 380 °C in the hydroprocessing of, e.g., heavy gasoil40 or scrap tires pyrolysis oil.41 Furthermore, the thermodynamic limitation of the hydrogenation reaction of aromatics has been experimentally observed (Figure S1). This figure depicts the evolution of the concentration of polyaromatics with the temperature at different values of space time and the tendencies that they will follow, showing that their concentration reaches a minimum at ca. 380 °C in all the cases, demonstrating once again the thermodynamic limitations of these hydrogenation processes. 3.1.3. Hydrocracking (HC). The hydrocracking activity of the catalyst is quantified following the evolution of the different fractions of products, based on the following boiling point criteria: naphtha (Na, TB < 489 K), diesel (D, 489 < TB < 623 K), and gasoil (G, TB > 623 K). The lack of LCO hydrocracking kinetic models in the literature is notorious, and the kinetic schemes for the hydrocracking of heavier streams have been taken as reference36,42 to propose the kinetic scheme shown in Figure 2. This scheme considers the interconversion of distinguished fractions, based on traditional

ij x Dx H2 yzz dxG zz − k 3xGx H = −k1jjjjxGx H2 − 2 j dτ Keq,G − D zz k {

(19)

ij ij x x yz x Dx H2 yzz dx D zz − k 2jjjx Dx H − Na H2 zzz = k1jjjjxGx H2 − 2 j jj dτ Keq,D − Na zz Keq,G − D zz k { k {

x Nax H2 zyz dx Na ji zz + k 3xGx H = k 2jjjjx Dx H2 − 2 j dτ Keq,D − Na zz k {

(20)

(21)

where Keq,G−D corresponds to the equilibrium between gasoil and diesel lumps, and Keq,D−N corresponds to the equilibrium between diesel and naphtha lumps. Furthermore, as previously explained for the kinetic modeling of HDS and HDA, the temperature dependence of the kinetic parameters has been expressed according to the reparameterized Arrhenius equation (eq 10). Exactly to the HDA model, temperature dependence correlations for the calculation of equilibrium parameters in the HC models have been obtained from experimental data obtained in the pseudoequilibrium state at high space time

Figure 2. Proposed reaction scheme for HC reactions of LCO. 13068

DOI: 10.1021/acs.iecr.9b02095 Ind. Eng. Chem. Res. 2019, 58, 13064−13075

Article

Industrial & Engineering Chemistry Research values (τ = 0.5 gcat h gLCO−1). The obtained empirical correlations have been Keq,G − D = 8.365·10−6T 2 − 1.142· 10−2T + 4.005

(22)

Keq,D − Na = 1.401·10−5T 2 − 2.027· 10−2T + 8.679

(23)

3.2. Kinetic Model Fitting and Validation. 3.2.1. Hydrodesulfurization (HDS). The apparent HDS rate constants for each reaction lump at reference temperature of 360 °C, as well as apparent activation energies, are listed in Table 2. As Table 2. Apparent Kinetic Parameters for the HDS Model of LCO lump a

M2BT M3BTb DBTc M1DBTd M2DBTe M3DBTf

k′j,360 °C (gLCO h−1 gcat−1) 88 83 24 21 16 35

± ± ± ± ± ±

4 3 7 5 3 6

Ej (kJ mol−1) 77 120 123 84 89 94

± ± ± ± ± ±

10 9 1 2 1 4

a

Dimethyl benzothiophenes. bTrimethyl benzothiophenes. cDibenzothiophene. dMethyl dibenzothiophenes. eDimethyl dibenzothiophenes. fTrimethyl dibenzothiophenes.

previously explained, we have assumed that there is not any parallel reaction involving two sulfur species and there is an inhibition effect of H2S (considered in the Langmuir− Hinshelwood equations). It should be mentioned that BT or M1BT have been undetected in the reaction products along the investigated operational range. Therefore, their kinetics have been assumed to be extremely fast, and the estimation of their kinetic parameters has not been possible. From the values reported in Table 2, it can be confirmed that hydrodesulfurization of M2BT and M3BT lumps is faster than that of refractory species, with kinetic parameters 4 times higher. Sulfur molecules with a higher number of C atoms in their substituents are more resistant toward HDS than less heavily substituted DBTs.7 However, comparing the reactivity of refractory sulfur species, M3DBT-type species exhibit a higher reaction rate than dibenzothiophene or the other alkyl derivates. This result shows how determining the position of alkyl substituents is, as substituents located close C−S bonds hinder their HDS activity and change the ratio between rates of the two HDS routes, i.e., direct and indirect desulfurization routes.49,50 The lower reactivity of M1DBT and M2DBT families can be explained because their alkyl substituents are mainly located in 4- and 4,6-positions, respectively. These species are the most refractory sulfur compounds, and both steric hindrance and electronic factors are claimed to be responsible for their low reactivity.15 With respect to apparent activation energies, values between 76 and 122 kJ mol−1 have been obtained, which are in concordance with the values reported in the literature for the removal of sulfur-containing molecules in petroleum-derived feedstocks.51,52 The quality and accuracy of the fitting can be evaluated from the comparison between the prediction of the model (lines) and the experimental data (symbols) plotted in Figure 3. It should be mentioned that, due to the different concentration of the sulfur species, the scale of ordinate axis is different for each graph. As can be observed, an increasing space time (higher catalyst loading) significantly reduces the sulfur concentration. Cao et al.53 and Wang et al.54 reported a very similar effect of

Figure 3. Comparison between the experimental data (symbols) and predicted data (lines) for the evolution of M2BT (a), M3BT (b), DBT (c), M1DBT (d), M2DBT (e), and M3DBT (f) lumps.

space time for NiMo catalysts supported on SBA-16 and mesoporous Al2O3, respectively, in the hydrodesulfurization of pure sulfur compounds. Analyzing the temperature effect in the results shown in Figure 3, it can be seen that the increase in HDS conversion for more refractory compounds (Figure 3d− f) is more pronounced in the range from 320 to 360 °C than from this latter temperature to 400 °C, which is not fully predicted by the model. This result is in concordance with the results obtained by Mann et al.55 on the hydrotreatment of heavy gasoil with a NiMo/γ-Al2O3 catalyst, proving that a higher increase is obtained in HDS conversion raising the temperature in the range 300−350 °C than in the range 350− 400 °C. In order to improve the goodness of fit of results depicted in Figure 3, the procedure has been repeated considering the equilibrium parameter for the H2S adsorption temperature dependent. However, the obtained results (not shown) have been worse than those previously depicted. The deviation of the fitting can be attributed to (i) the higher internal diffusional limitations of refractory compounds on the catalyst particles at higher temperatures, so that real activation energies are higher than the apparent ones calculated and displayed in Table 2; (ii) it has not been taken into account that desulfurization of highly refractory sulfur-containing molecules commonly follows a saturation pathway instead of a direct desulfurization route;56 and (iii) to the limitations of the employed kinetic equation to model the HDS, as probably 13069

DOI: 10.1021/acs.iecr.9b02095 Ind. Eng. Chem. Res. 2019, 58, 13064−13075

Article

Industrial & Engineering Chemistry Research

reactions being the hydrogenation of 1-ring aromatics and ring opening of naphthenes. It can, hence, be assumed that almost every 1-ring aromatic formed by partial hydrogenation of a 2+ring aromatics is further converted into a paraffin.58 These results indicate that aromatics with the “simplest” structure are transformed faster and that highly substituted and structurally more complex aromatics have a lower reactivity as a consequence of (i) the steric hindrance between the catalyst and the molecule during the adsorption and (ii) the inductive effect of bigger substituents in the benzene ring which increase its stability.37 On the other hand, according to the obtained apparent activation energies, both hydrogenation steps have a lower apparent activation energy, meaning that these reactions are less favored at higher temperatures, while ring opening reactions are enhanced at higher temperatures. This result is a consequence of the favored dehydrogenation and hydrogentransfer reactions at higher temperatures.5,59 Figure 4 offers a detailed evaluation of the capability of the kinetic model to reproduce the distribution of lumps in the

accounting for an adsorption term of the sulfur species would improve the fitting. A similar procedure for modeling the HDS of middle distillates has been followed by Albazzaz et al.57 as they have grouped all the sulfur-containing molecules into four different groups according to their reactivities. Furthermore, they have also taken into account inhibition effects caused by aromatics and N-containing molecules, concluding that each sulfur group has a different reactivity following a first-order kinetics. A simpler model based on multiple lumps has been developed by Yang et al.14 by grouping sulfur-containing molecules into benzothiophene, dibenzothiophene, and benzonaphthothiophene and their corresponding alkyl-homologues, also achieving a proper prediction of the HDS routes. 3.2.2. Hydrodearomatization (HDA). As explained above, a model discrimination must be performed between the three proposed kinetic schemes for HDA to identify the best performing one. To this purpose the most relevant statistical parameters related to the experimental data fitting for the three different kinetic models (Figure 1) are listed in Table 3. Table 3. Statistic Comparison of Three Kinetics Models (A, B, and C) Proposed for HDA of LCO statistical parameters

A

B

C

ne nl np ν SSE F F1−α

18 4 6 66 1.093 × 10−1

18 4 8 64 1.002 × 10−1 2.892 (FA−B) 3.140

18 4 10 62 1.009 × 10−1 1.285 (FA−C) 3.145

Regarding the sum of squared errors (SSE), model B exhibits the best fitting with a value of 1.002 × 10−1, while for models A and C values of 1.093 × 10−1 and 1.009 × 10−1 were obtained, respectively. On the other hand, it should be taken into account that both the number of experimental data points and lumps used for the experimental data fitting are the same for all three cases, even though the number of parameters increased from model A to C: np,A < np,B < np,C. Consequently, keeping in mind the definition of degrees of freedom (eq 8), νA > νB > νC. Hence, taking kinetic scheme A as the simpler one and looking for its improvement with models B and C, it has been obtained that FA−B < F1−α (2.892 < 3.140) and FA−C < F1−α (1.285 < 3.145), which indicates that neither model B nor model C is statistically better than model A, even though model A shows the worst fitting to experimental data (highest value of SSE). Table 4 lists the apparent kinetic parameter estimates for the HDA model (kinetic scheme A). As it can be seen, the parameters k2 and k3 are higher than k1, meaning that, 2+-ring aromatics are highly reduced with the subsequent formation of 1-ring aromatics and paraffins and isoparaffins, while the concentration of naphthenes remained almost constant. The proposed kinetic model reproduces these results, the fastest

Figure 4. Comparison between the experimental data (symbols) and predicted data (lines) for the evolution with space time of product composition on the hydroprocessing of LCO at 320 °C (a), 360 °C (b), and 400 °C (c).

liquid products. An accurate model (lines) to experimental data (symbols) fitting was achieved with the model also reproducing the marked temperature effect over lump concentration. Additionally, a lower amount of total paraffinic compounds is obtained upon increasing temperature, together with a higher amount of 1-ring aromatics, while naphthenes concentration does not present important variations, highlighting their role as intermediates in the HDA reaction scheme. For a given space time, a similar trend for aromatics was reported by Tang et al.60 for the low-temperature mild hydrotreatment of a coal distillate as product composition was seriously influenced by the shift of thermodynamic equilibrium. Several authors18,61 have reported the effect of space time on aromatic hydrogenation kinetics, predicting higher

Table 4. Apparent Kinetic Parameters for the HDA Kinetic Model of LCO (Scheme A in Figure 1) reaction

k′j,360 °C (gLCO h−1 gcat−1)

Ej (kJ mol−1)

A2+ → A1, k1 A1 → N, k2 N → P+iP, k3

(1.3 ± 1.0)·101 (5.0 ± 0.6)·102 (5.0 ± 0.5)·102

39.3 ± 8.8 10.9 ± 6.2 94.3 ± 0.7 13070

DOI: 10.1021/acs.iecr.9b02095 Ind. Eng. Chem. Res. 2019, 58, 13064−13075

Article

Industrial & Engineering Chemistry Research

increases to approximately 24−28 wt %. Therefore, the highest naphtha yield is obtained at 400 °C and 0.5 gcat h gLCO−1 (reaching a value of 40 wt %), coinciding with the lowest gasoil concentration in the products (6 wt %). It should be highlighted that for space times higher than 0.2 gcat h gLCO−1, the concentrations of the different fractions remain steady, with minimal variations with respect to the values obtained at 0.5 gcat h gLCO−1. The trend is analogue to that observed in Figure 4 when discussing HDA reaction pathways of the LCO. 3.3. Optimizing LCO Hydrotreating. The apparent kinetic parameters obtained for the three studied kinetic routes, namely (i) HDS, (ii) HDA, and (iii) HC, have allowed for the simulation of the process in the studied space time (0− 0.5 gcat h gLCO−1) and temperature (320−400 °C) ranges. The obtained results are depicted in Figure 6 using conversion

hydrogenation rates at higher space time values, in agreement with the results presented in Figure 4. 3.2.3. Hydrocracking (HC). Apparent kinetic parameter estimated at reference temperature (360 °C) and corresponding apparent activation energies obtained for the HC model are summarized in Table 5. As it can be observed, the gasoil-toTable 5. Apparent Kinetic Parameters for the HC Model of LCO reaction

k′j,360 °C (gLCO h−1 gcat−1)

Ej (kJ mol−1)

gasoil → diesel, k1 diesel → naphtha, k2 gasoil → naphtha, k3

(6.7 ± 0.3)·10 (1.1 ± 0.6)·101 (1.0 ± 0.1)·10−2

23.0 ± 2.8 23.2 ± 8.1 62.1 ± 5.5

1

diesel conversion rate is higher than the other rates, e.g., it is 5 times faster than the formation of naphtha from diesel. In addition, the extent of the gasoil-to-naphtha fraction conversion is lesser, as its direct transformation rate is various orders of magnitude lower (k3 ≪ k2 < k1). Based on this very low value for k3, the direct route of naphtha formation from gasoil can be deemed negligible. The statistical insignificance of this kinetic parameter (k3) has been ratified following the sensitivity analysis proposed by Alcázar and Ancheyta22 checking that the residual sum of squares is constant after perturbations of the parameter, which means that this additional reaction stage does not improve the model. A highly accurate HC model fitting (lines) and experimental data (symbols) was achieved, as observed in Figure 5, in the

Figure 6. Conversion maps obtained for hydrodesulfurization (a), total hydrodearomatization (b), polyaromatics hydrodearomatization (c), hydrocracking kinetic models (d), and zoomed area of hydrocracking conversion map (e).

Figure 5. Comparison between the experimental data (symbols) and predicted data (lines) for the evolution with space time of product lumps at 320 °C (a), 360 °C (b), and 400 °C (c).

maps. Note that in the case of HC an extra graph has been included (Figure 6d) where it has been zoomed in the range of 0−0.08 gcat h gLCO−1 and 320−400 °C in order to clear up the evolution followed by the HC conversion. All three simulated routes evidence the enhancing effect of both space time and temperature for the different hydroprocessing goals, favoring the removal of sulfur and aromatics components, as well as the conversion of gasoil toward lighter lumps. This way, at the highest investigated space time (0.5 gcat

whole range of investigated conditions. At all temperatures, diesel is the most abundant fraction, as it was in the LCO (Table 1). For space times below 0.2 gcat h gLCO−1, significant variations can be seen, especially regarding naphtha and diesel yields. Hence, when both temperature and space time are increased, diesel and gasoil fractions decreased in 14−15 and 10−13 wt %, respectively. Consequently, the naphtha fraction 13071

DOI: 10.1021/acs.iecr.9b02095 Ind. Eng. Chem. Res. 2019, 58, 13064−13075

Article

Industrial & Engineering Chemistry Research h gLCO−1), HDS conversions higher than 90% are obtained above 340 °C, while at 0.25 gcat h gLCO−1, temperatures higher than 354 °C are required to reach the same conversion level (Figure 6a). It is also noteworthy that in the range of 380−400 °C sulfur removal higher than 90% can be obtained for quite low values of space time (ca. 0.1 gcat h gLCO−1). In the case of HDA (Figure 6b), a local conversion minimum is observed at temperatures ca. 340 °C, caused by the thermodynamic control of hydrogenation processes, positively affected at temperatures above 340 °C. Attending to the effect of space time, an increase in the 0−0.15 gcat h gLCO−1 range causes a notable effect, especially at higher temperatures, although for higher space time values this influence becomes less relevant. Moreover, attending to Figure 6c it can be seen that the partial hydrogenation of polyaromatics occurs in a higher extent than that of monoaromatics, and, consequently, their conversion rates are much higher than the total ones. Hence, in order to maximize HDA reactions, temperature should preferably be in the 387−400 °C range, while the space time must be established in values higher than 0.25 gcat h gLCO−1. Finally, regarding HC conversion (Figure 6d), significantly high conversions can be obtained in almost the whole range of investigated operating conditions. Interestingly, at space times exceeding 0.1 gcat h gLCO−1 the effect of space time becomes negligible and that of the temperature becomes increasingly relevant. This way, for a given space time, the conversion that can be reached ranges from 50 to 65% by increasing the temperature from 320 to 400 °C. Having said that, it can be seen in the zoomed area in Figure 6e the strong influence of the temperature over the achieved conversion. Thus, for a fixed temperature of 360 °C a small increase in space time from 0.01 to 0.08 gcat h gLCO−1 results in an increase in conversion from 15 to 57%. According to these results, it can be concluded that for optimizing the operating conditions in the hydrotreating of LCO, a compromise between HDS, HDA, and HC conversions should be established. Figure 6 shows that severe conditions enhance all the reaction routes as they favor not only the removal of sulfur and aromatics but also the conversion of heavier gasoil molecules. However, not necessarily highly severe hydroprocessing conditions are required in order to obtain satisfactory results. In this case, in order to maximize the extent of reactions within the three goals, temperature should be established in the range 385−400 °C, but a space time amounting to 0.2 gcat h gLCO−1 could be established as an increase up to 0.5 gcat h gLCO−1 has no significant effect in terms of reached conversion. However, keeping in mind a future plant scale, operating with a bed loaded with a higher amount of catalyst would help to extend the life cycle of the catalyst bed while also keeping desired levels of conversion for longer times.61 Laredo et al.62 have performed an intensive review of the operating conditions employed in the hydrotreating of LCO. Along the lines of this review and with the results of this work, a two-stage strategy seems more plausible using non-noble metal-based catalysts in the first stage and precious metalbased catalysts in the second one. Focusing on the first stage, optimal temperatures and pressures for maximizing highquality fuels production are in the 350−425 °C and 41−100 bar ranges, respectively, which is in concordance with the obtained results in this study (385−400 °C temperature and 80 bar pressure). Regarding space time, values between 0.5 and 1 gcat h gLCO−1 have been selected in the literature as the

optimal ones. In our case, our proposed optimal space time value of 0.2 gcat h gLCO−1 allows for the use of a lower amount of catalyst, which would positively impact the process economics. Nevertheless, keeping in mind the strategy previously mentioned of extending the lifespan of the catalytic bed, our space time could be raised up to 0.5 gcat h gLCO−1 entering within the proposed range. Finally, the volumetric H2/LCO ratio used in this work has been of 1000 vol %, which is also located inside the interval tested in the literature (400− 1478 vol %). In summary, we can conclude that the proposed kinetic model for LCO hydrotreating (HDS, HDA, and HC) corroborates the previous experimental results.

4. CONCLUSIONS A kinetic model based on multiple lumps, species, and pathways has been proposed for the hydrotreating of LCO on a NiMo/SiO2−Al2O3 commercial catalyst. This model, which considers three different hydroprocessing goals (hydrodesulfurization, HDS; hydrodearomatization, HDA; and, hydrocracking, HC), proves to be appropriate for accurately reproducing the products distribution obtained. This kinetic model is of special interest for the design of large-scale hydroprocessing reactors for secondary refinery streams, in order to enhance the refining margin, boost carbon usage, and drop emissions. The HDS kinetic modeling of LCO, based on a Langmuir− Hinshelwood type reaction mechanism, allows for the simulation of the behavior of different sulfur-containing components and lumps. Overall, a good fitting of the model to experimental data has been obtained with light deviations of the most refractory compounds at high space time values, which are attributable to diffusional hindrances for heavy (condensed) molecules within the catalytic structure. For the HDA of the LCO, a sequential scheme formed by the following reactions is suitable: A2+ (2+-ring aromatics) ⇆ A1 (1-ring aromatics) ⇆ N (naphthenes) ⇆ P+iP (paraffins and isoparaffins). The equilibrium parameters of the proposed reverse reactions have been calculated from the pseudoequilibrium data experimentally obtained at high space times, leading to an accurate reproduction of the experimental data. Furthermore, a sequential model of formation, gasoil ⇆ diesel ⇆ naphtha, predicts the transformations between the different fractions in the hydrocracking of LCO, as it has been proved that the additional gasoil → naphtha transformation is not relevant. Employing the kinetic model expressions for the three reaction types, the optimal operating conditions for the first stage of the hydrotreatment of LCO using a NiMo/SiO2− Al2O3 commercial catalyst have been established in the 385− 400 °C temperature range and at a space time of 0.2 gcat h gLCO−1. Operating at these conditions, the hydrotreated LCO would have an optimal composition for a subsequent hydrocracking stage using noble metal-based catalysts, focused on enhancing the hydrocracking route, boosting the production of gasoline and diesel blends.



ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.iecr.9b02095. Products distribution and distillation cuts, evolution of concentration of polyaromatics with temperature (PDF) 13072

DOI: 10.1021/acs.iecr.9b02095 Ind. Eng. Chem. Res. 2019, 58, 13064−13075

Article

Industrial & Engineering Chemistry Research



AUTHOR INFORMATION

XHC, HHDA, and XHDS = conversion in the hydrocracking, hydrodearomatization, and hydrodesulfurization, respectively, reactions xi, xi*, and xi* = mass fraction of the lump i, computed, experimental, and average value of the experimental ones, respectively

Corresponding Author

*Phone: +34 946-015-341. E-mail. [email protected]. ORCID

Roberto Palos: 0000-0002-7716-3428 Alazne Gutiérrez: 0000-0002-5252-6691 Pedro Castaño: 0000-0002-6454-9321 Joris W. Thybaut: 0000-0002-4187-7904 José M. Arandes: 0000-0002-5644-5282

Greek symbols

α = level of significance of the kinetic model νA and νB = degrees of freedom of the models A and B, respectively τ = space time (gcat h gLCO−1) (υi)j = stoichiometric coefficient of lump i in reaction j

Present Addresses

§ Department of Chemical and Environmental Engineering, University of the Basque Country UPV/EHU, Plaza Europa 1, 20018 Donostia-San Sebastian, Spain. ⊥ Multiscale Reaction Engineering, KAUST Catalysis Center (KCC), King Abdullah University of Science and Technology (KAUST), Thuwal, 23955-6900, Saudi Arabia.



REFERENCES

(1) Barfuss, W.; Donges, J. F.; Lade, S. J.; Kurths, J. When Optimization for Governing Human-Environment Tipping Elements Is Neither Sustainable nor Safe. Nat. Commun. 2018, 9, 2354−2363. (2) Directive 2009/30/EC of the European Parliament and of the Council, Pub L 140/88 (Apr 23, 2009). (3) Galadima, A.; Muraza, O. Ring Opening of Hydrocarbons for Diesel and Aromatics Production: Design of Heterogeneous Catalytic Systems. Fuel 2016, 181, 618−629. (4) Xin, L.; Liu, X.; Chen, X.; Feng, X.; Liu, Y.; Yang, C. Efficient Conversion of Light Cycle Oil into High-Octane-Number Gasoline and Light Olefins over a Mesoporous ZSM-5 Catalyst. Energy Fuels 2017, 31, 6968−6976. (5) Escalona, G.; Rai, A.; Betancourt, P.; Sinha, A. K. Selective PolyAromatics Saturation and Ring Opening during Hydroprocessing of Light Cycle Oil over Sulfided Ni-Mo/SiO2-Al2O3 Catalyst. Fuel 2018, 219, 270−278. (6) Laredo, G. C.; Pérez-Romo, P.; Escobar, J.; Garcia-Gutierrez, J. L.; Vega-Merino, P. M. Light Cycle Oil Upgrading to Benzene, Toluene, and Xylenes by Hydrocracking: Studies Using Model Mixtures. Ind. Eng. Chem. Res. 2017, 56, 10939−10948. (7) Palos, R.; Gutiérrez, A.; Arandes, J. M.; Bilbao, J. Catalyst Used in Fluid Catalytic Cracking (FCC) Unit as a Support of NiMoP Catalyst for Light Cycle Oil Hydroprocessing. Fuel 2018, 216, 142− 152. (8) Weitkamp, J. Catalytic Hydrocracking-Mechanisms and Versatility of the Process. ChemCatChem 2012, 4, 292−306. (9) Lopez Abelairas, M.; De Oliveira, L. P.; Verstraete, J. J. Application of Monte Carlo Techniques to LCO Gas Oil Hydrotreating: Molecular Reconstruction and Kinetic Modelling. Catal. Today 2016, 271, 188−198. (10) Calderón, C. J.; Ancheyta, J. Modeling of CSTR and SPR Small-Scale Isothermal Reactors for Heavy Oil Hydrocracking and Hydrotreating. Fuel 2018, 216, 852−860. (11) Alvarez-Majmutov, A.; Chen, J. Stochastic Modeling and Simulation Approach for Industrial Fixed-Bed Hydrocrackers. Ind. Eng. Chem. Res. 2017, 56, 6926−6938. (12) Rodríguez, E.; Félix, G.; Ancheyta, J.; Trejo, F. Modeling of Hydrotreating Catalyst Deactivation for Heavy Oil Hydrocarbons. Fuel 2018, 225, 118−133. (13) Morales-Ortuño, J. C.; Ortega-Domínguez, R. A.; HernándezHipólito, P.; Bokhimi, X.; Klimova, T. E. HDS Performance of NiMo Catalysts Supported on Nanostructured Materials Containing Titania. Catal. Today 2016, 271, 127−139. (14) Yang, Y.; Dai, F.; Li, C.; Xiang, S.; Yaseen, M.; Zhang, S. Kinetic Evaluation of Hydrodesulfurization and Hydrodenitrogenation Reactions via a Lumped Model. Energy Fuels 2017, 31, 5491− 5497. (15) Morales-Valencia, E. M.; Castillo-Araiza, C. O.; Giraldo, S. A.; Baldovino-Medrano, V. G. Kinetic Assessment of the Simultaneous Hydrodesulfurization of Dibenzothiophene and the Hydrogenation of Diverse Polyaromatic Structures. ACS Catal. 2018, 8, 3926−3942. (16) Umana, B.; Zhang, N.; Smith, R. Development of Vacuum Residue Hydrodesulphurization-Hydrocracking Models and Their

Notes

The authors declare no competing financial interest.



ACKNOWLEDGMENTS The financial support of this work was undertaken by the Ministry of Economy and Competitiveness (MINECO) of the Spanish Government (Grant CTQ2015-67425R), the ERDF funds of the European Union, and the Basque Government (Grant IT748-13). I.H. is grateful for her postdoctoral grant awarded by the Department of Education, University and Research of the Basque Government (Grant POS_2015_1_0035). The authors are thankful for the technical and human support provided by SGIker of UPV/ EHU and European funding (ERDF and ESF). The authors also acknowledge Petronor Refinery for providing the feedstock used in this work.



NOMENCLATURE A1, A2+, N, and P+iP = 1-ring aromatics, 2+-ring aromatics, naphthenes, and paraffin and isoparaffins, respectively BT, MxBTs, DBT, and MxDBTs = benzothiophene, methylsubstituted benzothiophenes, dibenzothiophene, and methyl-substituted dibenzothiophenes, respectively Ej = apparent activation energy (kJ mol−1) D, G, and Na = diesel, gasoil, and naphtha, respectively F and F1−α = distribution and critical value of the Fisher’s exact test, respectively HC, HDA, and HDS = hydrocracking, hydrodearomatization, and hydrodesulfurization, respectively Keq and KH2S = equilibrium parameter and H2S adsorption equilibrium parameter, respectively kj and k′j = intrinsic and apparent kinetic coefficient of reaction j, respectively (gLCO h−1 gcat−1) LCO = light cycle oil ne, nl, and np = number of experiments, lumps, and parameters, respectively OF = objective function R = universal gas constant ri and rj = formation rate of lump i and rate of reaction j, respectively SSE = sum of squared errors (S) = sulfur content (ppm) T, TB, and Tref = temperature, boiling point, and reference temperature (K) TOS = time on stream (h) 13073

DOI: 10.1021/acs.iecr.9b02095 Ind. Eng. Chem. Res. 2019, 58, 13064−13075

Article

Industrial & Engineering Chemistry Research Integration with Refinery Hydrogen Networks. Ind. Eng. Chem. Res. 2016, 55, 2391−2406. (17) Srinivas, B. K.; Pant, K. K.; Gupta, S. K.; Saraf, D. N.; Choudhury, I. R.; Sau, M. A Carbon-Number Lump Based Model for Simulation of Industrial Hydrotreaters: Vacuum Gas Oil (VGO). Chem. Eng. J. 2019, 358, 504−519. (18) Hita, I.; Aguayo, A. T.; Olazar, M.; Azkoiti, M. J.; Bilbao, J.; Arandes, J. M.; Castaño, P. Kinetic Modeling of the Hydrotreating and Hydrocracking Stages for Upgrading Scrap Tires Pyrolysis Oil (STPO) toward High-Quality Fuels. Energy Fuels 2015, 29, 7542− 7553. (19) Palos, R.; Gutiérrez, A.; Arandes, J. M.; Bilbao, J. Upgrading of High-Density Polyethylene and Light Cycle Oil Mixtures to Fuels via Hydroprocessing. Catal. Today 2018, 305, 212−219. (20) Gutiérrez, A.; Arandes, J. M.; Castaño, P.; Aguayo, A. T.; Bilbao, J. Role of Acidity in the Deactivation and Steady Hydroconversion of Light Cycle Oil on Noble Metal Supported Catalysts. Energy Fuels 2011, 25, 3389−3399. (21) Toch, K.; Thybaut, J. W.; Marin, G. B. A Systematic Methodology for Kinetic Modeling of Chemical Reactions Applied to n -Hexane Hydroisomerization. AIChE J. 2015, 61, 880−892. (22) Alcázar, L. A.; Ancheyta, J. Sensitivity Analysis Based Methodology to Estimate the Best Set of Parameters for Heterogeneous Kinetic Models. Chem. Eng. J. 2007, 128, 85−93. (23) Riazi, M. R.; Roomi, Y. A. A Method to Predict Solubility of Hydrogen in Hydrocarbons and Their Mixtures. Chem. Eng. Sci. 2007, 62, 6649−6658. (24) Riazi, M. R.; Vera, J. H. Method to Calculate the Solubilities of Light Gases in Petroleum and Coal Liquid Fractions on the Basis of Their P/N/A Composition. Ind. Eng. Chem. Res. 2005, 44, 186−192. (25) Cordero-Lanzac, T.; Aguayo, A. T.; Gayubo, A. G.; Castaño, P.; Bilbao, J. Simultaneous Modeling of the Kinetics for N-Pentane Cracking and the Deactivation of a HZSM-5 Based Catalyst. Chem. Eng. J. 2018, 331, 818−830. (26) Sau, M.; Narasimhan, C. S. L.; Verma, R. P. A Kinetic Model for Hydrodesulfurisation. Stud. Surf. Sci. Catal. 1997, 106, 421−435. (27) García-Martínez, J. C.; Castillo-Araiza, C. O.; De los Reyes Heredia, J. A.; Trejo, E.; Montesinos, A. Kinetics of HDS and of the Inhibitory Effect of Quinoline on HDS of 4,6-DMDBT over a Ni-MoP/Al2O3 Catalyst: Part I. Chem. Eng. J. 2012, 210, 53−62. (28) Yang, L.; Li, X.; Wang, A.; Prins, R.; Chen, Y.; Duan, X. Hydrodesulfurization of Dibenzothiophene, 4,6-Dimethyldibenzothiophene, and Their Hydrogenated Intermediates over Bulk Tungsten Phosphide. J. Catal. 2015, 330, 330−343. (29) Dorneles de Mello, M.; de Almeida Braggio, F.; da Costa Magalhães, B.; Zotin, J. L.; da Silva, M. A. P. Kinetic Modeling of Deep Hydrodesulfurization of Dibenzothiophenes on NiMo/Alumina Catalysts Modified by Phosphorus. Fuel Process. Technol. 2018, 177, 66−74. (30) Rana, M. S.; Ancheyta, J.; Rayo, P.; Maity, S. K. Heavy Oil Hydroprocessing over Supported NiMo Sulfided Catalyst: An Inhibition Effect by Added H2S. Fuel 2007, 86, 1263−1269. (31) Castañeda, L. C.; Muñoz, J. A. D.; Ancheyta, J. Comparison of Approaches to Determine Hydrogen Consumption during Catalytic Hydrotreating of Oil Fractions. Fuel 2011, 90, 3593−3601. (32) Rodríguez, M. A.; Ancheyta, J. Modeling of Hydrodesulfurization (HDS), Hydrodenitrogenation (HDN), and the Hydrogenation of Aromatics (HDA) in a Vacuum Gas Oil Hydrotreater. Energy Fuels 2004, 18, 789−794. (33) Froment, G. F. Modeling in the Development of Hydrotreatment Processes. Catal. Today 2004, 98, 43−54. (34) Cerqueira, H. S.; Mihindou-Koumba, P. C.; Magnoux, P.; Guisnet, M. Methylcyclohexane Transformation over HFAU, HBEA, and HMFI Zeolites: I. Reaction Scheme and Mechanisms. Ind. Eng. Chem. Res. 2001, 40, 1032−1041. (35) Castaño, P.; Arandes, J. M.; Pawelec, B.; Fierro, J. L. G.; Gutiérrez, A.; Bilbao, J. Kinetic Model Discrimination for Toluene Hydrogenation over Noble-Metal-Supported Catalysts. Ind. Eng. Chem. Res. 2007, 46, 7417−7425.

(36) Gutiérrez, A.; Castaño, P.; Azkoiti, M. J.; Bilbao, J.; Arandes, J. M. Modelling Product Distribution of Pyrolysis Gasoline Hydroprocessing on a Pt-Pd/HZSM-5 Catalyst. Chem. Eng. J. 2011, 176− 177, 302−311. (37) Deligny, J.; Germanaud, L.; Dath, J.-P. J.; Brunet, S. Hydrodearomatization of Model Monoaromatics over Ni/Al2O3: Theoretical and Experimental Approaches. Catal. Lett. 2018, 148, 2548−2560. (38) Gutiérrez, A.; Arandes, J. M.; Castaño, P.; Olazar, M.; Bilbao, J. Preliminary Studies on Fuel Production through LCO Hydrocracking on Noble-Metal Supported Catalysts. Fuel 2012, 94, 504−515. (39) Carrero-Mantilla, J.; Llano-Restrepo, M. Vapor-Phase Chemical Equilibrium for the Hydrogenation of Benzene to Cyclohexane from Reaction-Ensemble Molecular Simulation. Fluid Phase Equilib. 2004, 219, 181−193. (40) Mapiour, M.; Sundaramurthy, V.; Dalai, A. K.; Adjaye, J. Effects of the Operating Variables on Hydrotreating of Heavy Gas Oil: Experimental, Modeling, and Kinetic Studies. Fuel 2010, 89, 2536− 2543. (41) Hita, I.; Gutiérrez, A.; Olazar, M.; Bilbao, J.; Arandes, J. M.; Castaño, P. Upgrading Model Compounds and Scrap Tires Pyrolysis Oil (STPO) on Hydrotreating NiMo Catalysts with Tailored Supports. Fuel 2015, 145, 158−169. (42) Ancheyta, J.; Sánchez, S.; Rodríguez, M. A. Kinetic Modeling of Hydrocracking of Heavy Oil Fractions: A Review. Catal. Today 2005, 109, 76−92. (43) Novaes, L. d. R.; Pacheco, M. E.; Martins Salim, V. M.; de Resende, N. S. Accelerated Deactivation Studies of Hydrotreating Catalysts in Pilot Unit. Appl. Catal. A Gen. 2017, 548, 114−121. (44) Kohli, K.; Prajapati, R.; Maity, S. K.; Sau, M.; Sharma, B. K. Accelerated Pre-Coking of NiMo/γ-Al2O3 Catalyst: Effect on the Hydroprocessing Activity of Vacuum Residue. Fuel 2019, 235, 437− 447. (45) Birtill, J. J. But Will It Last until the Shutdown? Deciphering Catalyst Decay! Catal. Today 2003, 81, 531−545. (46) Silvy, R. P. Future Trends in the Refining Catalyst Market. Appl. Catal., A 2004, 261, 247−252. (47) Marafi, A.; Almarri, M.; Stanislaus, A. The Usage of High Metal Feedstock for the Determination of Metal Capacity of ARDS Catalyst System by Accelerated Aging Tests. Catal. Today 2008, 130, 395− 404. (48) Al Bazzaz, H.; Kang, J.-L.; Chehadeh, D.; Bahzad, D.; Wong, D. S.-H.; Jang, S.-S. Robust Predictions of Catalyst Deactivation of Atmospheric Residual Desulfurization. Energy Fuels 2015, 29, 7089− 7100. (49) Zhao, H.; Oyama, S. T.; Freund, H.-J.; Włodarczyk, R.; Sierka, M. Nature of Active Sites in Ni2P Hydrotreating Catalysts as Probed by Iron Substitution. Appl. Catal., B 2015, 164, 204−216. (50) Toutov, A. A.; Salata, M.; Fedorov, A.; Yang, Y.-F.; Liang, Y.; Cariou, R.; Betz, K. N.; Couzijn, E. P. A.; Shabaker, J. W.; Houk, K. N.; et al. A Potassium Tert-Butoxide and Hydrosilane System for Ultra-Deep Desulfurization of Fuels. Nat. Energy 2017, 2, 17008. (51) Sheng, Q.; Wang, G.; Zhang, Q.; Gao, C.; Ren, A.; Duan, M.; Gao, J. Seven-Lump Kinetic Model for Non-Catalytic Hydrogenation of Asphaltene. Energy Fuels 2017, 31, 5037−5045. (52) Novaes, L. da R.; de Resende, N. S.; Salim, V. M. M.; Secchi, A. R. Modeling, Simulation and Kinetic Parameter Estimation for Diesel Hydrotreating. Fuel 2017, 209, 184−193. (53) Cao, Z.; Du, P.; Duan, A.; Guo, R.; Zhao, Z.; Zhang, H.; Zheng, P.; Xu, C.; Chen, Z. Synthesis of Mesoporous Materials SBA-16 with Different Morphologies and Their Application in Dibenzothiophene Hydrodesulfurization. Chem. Eng. Sci. 2016, 155, 141−152. (54) Wang, X.; Zhao, Z.; Zheng, P.; Chen, Z.; Duan, A.; Xu, C.; Jiao, J.; Zhang, H.; Cao, Z.; Ge, B. Synthesis of NiMo Catalysts Supported on Mesoporous Al2O3 with Different Crystal Forms and Superior Catalytic Performance for the Hydrodesulfurization of Dibenzothiophene and 4,6-Dimethyldibenzothiophene. J. Catal. 2016, 344, 680− 691. 13074

DOI: 10.1021/acs.iecr.9b02095 Ind. Eng. Chem. Res. 2019, 58, 13064−13075

Article

Industrial & Engineering Chemistry Research (55) Mann, R. S.; Sambi, I. S.; Khulbe, K. C. Catalytic Hydrofining of Heavy Gas Oil. Ind. Eng. Chem. Res. 1987, 26, 410−414. (56) Ding, S.; Li, A.; Jiang, S.; Zhou, Y.; Wei, Q.; Zhou, W.; Huang, Y.; Yang, Q.; Fan, T. Niobium Modification Effects on Hydrodesulfurization of 4,6-DMDBT Catalyzed on Ni-Mo-S Active Sites: A Combination of Experiments and Theoretical Study. Fuel 2019, 237, 429−441. (57) Albazzaz, H.; Marafi, A. M.; Ma, X.; Ansari, T. Hydrodesulfurization Kinetics of Middle Distillates: A Four-Lumping Model with Consideration of Nitrogen and Aromatics Inhibitions. Energy Fuels 2017, 31, 831−838. (58) Martínez, A.; Arribas, M. A.; Pergher, S. B. C. Bifunctional Noble Metal/Zeolite Catalysts for Upgrading Low-Quality Diesel Fractions via Selective Opening of Naphthenic Rings. Catal. Sci. Technol. 2016, 6, 2528−2542. (59) Lozano, L.; Marin, G. B.; Thybaut, J. W. Analytical Rate Expressions Accounting for the Elementary Steps in Benzene Hydrogenation on Pt. Ind. Eng. Chem. Res. 2017, 56, 12953−12962. (60) Tang, W.; Fang, M.; Wang, H.; Yu, P.; Wang, Q.; Luo, Z. Mild Hydrotreatment of Low Temperature Coal Tar Distillate: Product Composition. Chem. Eng. J. 2014, 236, 529−537. (61) Gutiérrez, A.; Arandes, J. M.; Castaño, P.; Olazar, M.; Barona, A.; Bilbao, J. Effect of Space Velocity on the Hydrocracking of Light Cycle Oil with a Pt-Pd/HY Zeolite Catalyst. Fuel Process. Technol. 2012, 95, 8−15. (62) Laredo, G. C.; Vega Merino, P. M.; Hernández, P. S. Light Cycle Oil Upgrading to High Quality Fuels and Petrochemicals: A Review. Ind. Eng. Chem. Res. 2018, 57, 7315−7321.

13075

DOI: 10.1021/acs.iecr.9b02095 Ind. Eng. Chem. Res. 2019, 58, 13064−13075