The Inﬂuence of Hydrodearomatisation Reaction Kinetics on the Modelling of Sulphur and Aromatics Removal from Diesel Fuel in an Industrial Hydrotreating Process

: Over the years, the hydrotreating process has been considerably improved in order to facilitate the production of environmentally friendly diesel fuels by reducing sulphur and aromatics concentrations, as mandated by contemporary emissions regulations. In this study, different kinetic models for the hydrodearomatisation reaction and the inﬂuence of reaction rate on performance of the industrial trickle bed reactor for hydrotreating of gas oil and light cycle oil fractions were analysed. The impact on reactor temperature, catalyst wetting efﬁciency, and conversion of sulphur and aromatics were determined. The results of simulations were compared with experimental data from an industrial test run and the best model for the observed process is proposed. Reactor performance and overall efﬁciency of the process is strongly dependent on the kinetics of hydrodearomatisaton with respect to aromatics conversion but even more so with respect to the temperature increase in the reactor, which affects all key catalytic reaction parameters, catalyst wetting efﬁciency, and thus the sulphur conversion. Based on the obtained simulation results, it could be concluded that reactor performance is strongly dependent on the hydrodearomatisation reaction. The best predictions of outlet temperature as well as sulphur and aromatic conversion (deviation from the experimental value 0.87 K, 0.01% and 2.57%, respectively) are achieved with the Langmuir–Hinshelwood kinetic model proposed by Owusu-Boakye.


Introduction
Impurities like sulphur and aromatics contained in refined petroleum products lower the quality of fuel and have an important impact on soot emissions with serious effects on human health and potential carcinogenic effects [1][2][3][4]. Several studies have shown that undesired soot emissions from diesel fuel are proportional to the aromatics content of the fuel [2][3][4][5].
The current sulphur concentration limits are very low, and consequently the process of desulphurization is highly optimised to reduce sulphur content, but elimination of aromatic compounds is still difficult, especially for diesel feedstocks with high aromatics content. Nowadays, more attention is paid to the content of aromatics in the product, and aromatics specifications for diesel fuel vary widely (depending on the location and feedstock composition). According to the regulatory standard for fuel specifications, the aromatic hydrocarbon content should not exceed 10 vol%, while at the same time, for small refineries, the reference fuel specification limit is 20 vol% (aromatic hydrocarbon content ≤21.0 wt% and polycyclic aromatic hydrocarbon content ≤3.5 wt% by the D5186-96 test method). The saturation and content of aromatic compounds directly impact the cetane number and the fuel emissions characteristics [3].
Stability and quality of the fuels is improved in the process of hydrotreating as it is the most important process in refineries for sulphur removal and aromatics reduction through the reactions of hydrodesulphurisation (HDS) and hydrodearomatisaton (HDA). In the literature, a broad range of different research papers analysing the hydrotreating process can be found [6][7][8][9][10][11][12] in which hydrodesulphurisation and hydrodearomatisation reactions in trickle-bed reactors are the key step [6,7]. Study and further improvement of hydrotreating processes could be analysed through adequate investigations of kinetics and reactor modelling. The previously developed mathematical model used in this study was developed with input data for complex mixtures of hydrocarbons being processed in industrial units and not from studies of bench-scale and pilot-scale processes, which are often found to be unreliable for industrial-scale applications. A reliable and well-established mathematical model can be crucial for process design, determination of optimum operating parameters, for the understanding of the entire process, and for a realistic representation of the process [6][7][8]. Complexity of the mathematical model implies: detailed analysis of the inlet and outlet mixtures, well-established material and energy balances, the adequate kinetic equations and parameters, the phase equilibrium and thermodynamic calculations for multicomponent mixtures, the rate equations of reactions that are taking place in very complex reaction mixtures in liquid and vapour phases, the evaporation rate as well as the effectiveness factors, and catalyst wetting efficiency. In such a complex reaction network, the reactivity of components is affected by many factors [6,7,9].
Comparisons of the reactivity of different aromatic compounds have been reported in a review study [10]. Most of the quantitative reactivity data characterise reactions of individual aromatic hydrocarbons with hydrogen and few studies have been reported with mixtures of hydrocarbons chosen so that inhibition effects were negligible. However, inhibition effects due to competitive adsorption can affect reactivity significantly. Several review studies about HDA reaction networks and kinetics can be found in the literature [11][12][13].
Hydrodearomatisation reactions are reactions in which aromatic rings are saturated with hydrogen. For polycondensed aromatic hydrocarbons, the hydrogenation of the first ring is in general the fastest step, and the rate of hydrogenation for subsequent rings tend to be lower-rate steps, with the last ring being the least reactive [12,13]. The rate of hydrogenation of the last ring is significantly lower than that of the first one. Overall, aromatic conversion is therefore often limited by thermodynamic equilibrium, since inside the liquid-filled pores of the catalyst particles diffusion of hydrogen from the gas phase through a liquid film and into the particle has to occur before a reaction can take place. Phenanthrene is a good model compound to represent polyaromatics. An important aspect of the hydrodearomatisation reactions is that in typical hydrotreating conditions the conversion can be limited by the thermodynamic equilibrium and by mass transport of the reacting species depending on the process conditions.
There have been several attempts toward a thorough description of the HDA reaction, which can be found in the literature [14][15][16][17][18][19][20][21][22][23][24][25][26]. Most of the HDA kinetic models account for one sum reaction and very few are developed to accommodate three classes of reacting aromatic compounds, mono-, di-, and triaromatics. Kinetic models for the HDA reaction in a gas/oil system mainly assume that hydrogenation and dehydrogenation reactions occur according to the Langmuir-Hinshelwood mechanisms and the HDA reaction is represented as a first-order reversible reaction. Cheng et al. [27] investigated the performance of a fixed-bed reactor in concurrent and counter-current flows to remove sulphur and aromatics in diesel fuel. The model presented by this group is a one-dimensional heterogeneous model that accounts for the HDS and HDA reactions to simulate the concentration profiles of the reactants and products in the gas, liquid, and solid phases.
A complex and reliable mathematical model for hydrotreating of straight run gas oil (SRGO) blended with fluid catalytic cracking naphtha and light cycle oil (FCC-N-LCO) in an industrial trickle bed reactor, where several types of sulphur and aromatic compounds coexist and several classes of reactions occur simultaneously, was previously developed and verified [6,7]. In this study, a previously developed model [6] was used for the investigation of different kinetic models for the HDA reaction. Different kinetic models for the HDA reaction were incorporated in the model [6] and the results of the simulations were analysed in order to obtain the most accurate prediction of reactor performance. Concentrations of sulphur and aromatic compounds in products along with temperature difference in the reactor were verified by comparison with industrial test run data.

Industrial Test Run Data
The industrial test run was carried out in a catalytic hydrotreater reactor, which is normally a steady-state operating tubular flow reactor in which an inlet stream of a straight run gas oil stream was blended with a stream of light cycle oil (LCO) and fluid catalytic cracking naphtha (FCC N) upstream of the reactor. The blended inlet stream was hydrotreated in the presence of a commercial Co-Mo/γ-Al 2 O 3 catalyst. The inlet stream was analysed using GC-MS and different classes of sulphur compounds were found to be present in the stream: substituted benzothiophenes (C1-C4-BTs), dibenzothiophene and naphthothiophene (DBT and NT), and substituted DBTs and NTs (C1-C3-DBTs and -NTs). Properties of the feedstock and output stream used in the test run relevant for the model development and validation are given in the Table 1. All chemical characterisations and measurements were obtained using methods specified by the relevant technical standards, Table 1 in [7]. The composition of the initial feed into the catalytic reactor for this experiment was 87 vol% SRGO and 13 vol% FCC N-LCO. Pressure in the reactor was 40 bar, inlet temperature 606 K, and amount of the catalyst 18,000 kg. A more detailed description of the reactor, catalyst, and the test run can be found in the literature [6,7]. Experimental points used for the model validation in this study were: overall conversion of sulphur of 99.49 wt%, overall conversion of aromatics of 70 wt%, and the reactor temperature difference of 11.5 K.

The Reactor Model
The developed mathematical model used for the simulation of the industrial hydrotreater is deterministic and one-dimensional. The model takes into consideration the volatility of the reaction mixture and therefore the HDS and HDA reactions in the liquid and in the vapour phases. The composition of the reaction mixture was approximated using pseudocomponents according to the distillation curve of the reaction mixture. Pseudocomponents for thermodynamic representation in the vapour liquid equilibrium of different classes of sulphur compounds (BT, DBT1, DBT2, DBT3) were: n-dodecane, n-hexadecane, and n-octadecane, respectively. Aromatics were represented by toluene, tetraline, and phenanthrene as pseudo components. Two additional components in this model were methane and hydrogen as the main constituents of the gaseous phase. Toluene represents the monoaromatics and tetraline stands for diaromatics typically present in diesel fractions. Phenantrene was included to provide information on the successive hydrogenation steps of polyaromatic compounds.
The complex reaction network is even more complicated when vapour-liquid equilibrium (VLE) is included. VLE is crucial in the hydrotreater simulation due to the volatile nature of the oil and high solubility of hydrogen in the oil [6,9]. The presence of reactants in both phases creates a complex reaction system. In this study, vapour-liquid equilibria through the catalyst bed were simulated using the Peng-Robinson equation of state (PR-EOS) with thermodynamic properties and components distribution calculated by UniSim and AspenPlus software. The model assumes trickle flow in the adiabatic catalytic reactor with negligible axial dispersion, while catalyst deactivation during the industrial test run was neglected. Material balance equations in the trickle-bed reactor (TBR) for all reacting components in vapour and in liquid phases are described by the set of ordinary differential equations (ODEs). The material and overall energy balance equations, calculation of catalyst wetting efficiency, and overall catalyst effectiveness factors for vapour and liquid phases are incorporated in the model and described in a previously published paper [6]. Below are the summarised main equations used in the model ( Table 2). All equations and systems of equations were solved in a program developed in MATLAB ® & Simulink ® Release 2010b. Table 2. Equations used in the model.

Parameter: Equation:
Equilibrium constants Liquid-to-vapour flow ratios Densities of the vapour and liquid phases

Heat capacities
Catalyst wetting efficiency Overall catalyst effectiveness; the internal effectiveness factor; generalized Thiele modulus Material balance equations for reacting components in the process of hydrodearomatisation The overall energy balance

Hydrodearomatisation Reactions
Kinetic equations for different classes of sulphur compounds involved in the HDS reaction network are represented by the Langmuir-Hinshelwood rate equations proposed in the literature [19][20][21][22]. These equations were developed on a commercial CoMo/Al 2 O 3 catalyst under operating conditions significant to industrial applications. The kinetic equations for hydrogenation of aromatics used in this study are summarized within Table 3.
Hydrogenation of aromatic rings is reversible, highly exothermic, and a very important reaction in hydroprocessing, and therefore development of HDA reaction rate equations requires great attention. Several researchers have proposed models based on a simple first-order reversible reaction. Typical industrial feeds contain mixtures of tri-, di-, and monoaromatics, and hydrogenation proceeds via consecutive reversible reactions; however. such an approach has certain limitations. The HDA reaction is equilibrium-limited at high temperatures, under which the reverse reaction of naphthalene dehydrogenation occurs. Since HDA reactions are reversible, therefore, an equilibrium conversion is the upper limit for conversion. Cheng and co-authors [27] found that 380 • C was the upper temperature limit that can be practically employed.
Kinetic expressions used in this study are proposed in Table 3 and briefly described below.
Chowdhury [14] experimentally investigated desulphurisation and hydrogenation of aromatics in diesel oil in an isothermal lab-scale TBR where three groups of aromatics were considered: mono-, di-, and polyaromatics. All polyaromatics (tri-, tetra-, and pentaaromatics) behave like triaromatics. In this model, HDA reactions are first-order reversible with the assumption of a completely wetted catalyst bed. The heat of HDA reactions is set at 67 kJ/mol. The temperature dependencies of all the HDA reaction rate constants have been described by the Arrhenius law, where activation energies and the frequency factors have been determined via regression analysis.
Avraam [19] developed a model for trickle-bed reactors for the hydroprocessing of light oil feedstocks containing volatile compounds where reactions of desulphurisation, denitrogenation, saturation of olefins, and hydrogenation of aromatics were considered. Hydrogenation of aromatics was treated as a reversible reaction where A1 is hydrogenated into a totally or partially saturated compound, A2, with a substantial heat of reaction of 180 kJ/mol. This model was validated using pilot-scale data.
In the work of Owusu-Boakye [22], experiments in a lab-scale TBR for single-and two-stage hydrogenation of aromatic compounds in light gas oil were performed. Using the single-site mechanism form of the Langmuir-Hinshelwood rate of reaction, a kinetic model for the HDA reaction was developed. The apparent kinetic parameters of the single-and two-stage processes were obtained using the nonlinear least squares approach. Apparent activation energy and heats of adsorption for calculation of kinetic constants were determined directly from the slopes of the Arrhenius and Van't Hoff plots. Heat of reaction used in this model is 85 kJ/mol.
Yui [23] developed a simple power law kinetic model for HDA reactions based on first-order reversible reactions. Kinetic parameters are calculated using published catalytic aromatics hydrogenation data from a variety of sources. Both forward-and reverse-rate constants are calculated, and heat of reaction used is 255 kJ/mol. De Oliveira [24] brought a novel approach for kinetic modelling of hydrotreating reactions. In this study, hydrotreating of LCO gas oil was simulated using a stochastic simulation method. In view of the fact that the HDA is a monomolecular reaction, the stochastic rate constants of forward and reverse HDA reactions were calculated using quantitative structure/reactivity correlations (QS/RCs). Heat of HDA used in the de Oliveira model is 51 kJ/mol.
Tri ·η Tri −r Tri = k a3 ·C n a3 Tri ·η Tri − k a4 ·C n a4 Tetra ·η Tetra −r Tetra = k a4 ·C n a4 Tetra ·η Tetra A dynamic mathematical model for the commercial hydrotreating process, including hydrodesulphurisation, hydrodenitrogenation, and hydrodearomatisation, was developed by Liu and co-workers [25]. The hydrogenation of each aromatic group (mono-, di-, polyaromatics) is presented as power-low reaction kinetics. The reaction rate constants are calculated according to the Arrhenius law and heat of reaction used is 125 kJ/mol.
Values of calculated kinetic parameters for the inlet temperature in a reactor are given in Table 4. Table 4. Calculated kinetic parameters and rate constants for proposed models.

Model
Kinetic Constant Ea, kJ/mol Value Simulations using different kinetic models were carried out in a user-friendly program developed in MATLAB and Simulink Release 2010b. The system of ordinary differential equations was solved using a fourth-order Runge-Kutta routine with variable step size. Results of simulations were compared mutually and with experimental data.

Results and Discussion
As the result of simulations, molar flows of different classes of sulphur and aromatic components in liquid and in vapor phases were calculated, as well as temperature changes through the reactor, catalyst wetting efficiency, and effectiveness factors. All variables were calculated through the reactor as a function of catalyst weight. Key results to be observed were overall conversions of aromatics and sulphur compounds, since the final goal of the hydrotreating process is the reduction of aromatics and removal of sulphur in the final product. Figure 1 shows the conversions of aromatics calculated using different HDA kinetic models.
Experimental data of total aromatics conversion from the industrial test run was compared to all proposed kinetic models. The general conclusion is that there is a strong difference in the results of simulations when Langmuir-Hinshelwood kinetic models are used when compared to the power law kinetic models.
In Figure 1, conversions of total aromatics can be compared. The reaction rate of HDA is higher as the number of aromatic rings increases (comparing kinetic constants in Models 1 and 6, Table 4). Therefore, the HDA reaction of monoaromatic molecules was the slowest, the rate of hydrogenation of diaromatic molecules was faster, and tri-and polyaromatic molecules were the most reactive compounds. In addition, many studies of HDA reactions have shown that no partially hydrogenated ring compounds were observed in the product analyses, which proves that HDA occurs in a ring-by-ring manner [8,28].  [14], Avraam [19], Owusu-Boakye [22], Yui [23], de Oliveira [24], Liu [25], and experimental data at the reactor outlet for total aromatics conversion.
Experimental data of total aromatics conversion from the industrial test run was compared to all proposed kinetic models. The general conclusion is that there is a strong difference in the results of simulations when Langmuir-Hinshelwood kinetic models are used when compared to the power law kinetic models.
In Figure 1, conversions of total aromatics can be compared. The reaction rate of HDA is higher as the number of aromatic rings increases (comparing kinetic constants in Models 1 and 6, Table 4). Therefore, the HDA reaction of monoaromatic molecules was the slowest, the rate of hydrogenation of diaromatic molecules was faster, and tri-and polyaromatic molecules were the most reactive compounds. In addition, many studies of HDA reactions have shown that no partially hydrogenated ring compounds were observed in the product analyses, which proves that HDA occurs in a ring-by-ring manner [8,28].
The best prediction of total aromatics conversion in this study was achieved using the Langmuir-Hinshelwood kinetic model for HDA reactions proposed by Owusu Boakye (Figure 1). The kinetic model proposed by Chowdhury showed the highest deviation from the experimental point due to the lowest energy activation and consequently the low temperature increase in the reactor (Table 4 and Figures 1 and 2). The power law kinetic model proposed by Yui, as shown in Figure 1, resulted in considerably higher deviation from experimental point.
Since HDS and HDA reactions occur simultaneously and in parallel, in any hydrotreater reactor it is very important to monitor both sulphur and aromatics conversion, as well as temperature increase in the reactor (impacted by both exothermic reactions). Overall sulphur conversion along the catalyst bed is shown in Figure 2. Almost all of the applied HDA kinetic models predicted correctly the reduction of sulphur content in the reaction mixture along the catalyst bed. It can be concluded that all HDA kinetic models except those proposed by Liu and Avraam match very well with the experimental data for sulphur conversion (Figure 2a).  [14], Avraam [19], Owusu-Boakye [22], Yui [23], de Oliveira [24], Liu [25], and experimental data at the reactor outlet for total aromatics conversion.
The best prediction of total aromatics conversion in this study was achieved using the Langmuir-Hinshelwood kinetic model for HDA reactions proposed by Owusu Boakye (Figure 1). The kinetic model proposed by Chowdhury showed the highest deviation from the experimental point due to the lowest energy activation and consequently the low temperature increase in the reactor (Table 4 and Figures 1 and 2). The power law kinetic model proposed by Yui, as shown in Figure 1, resulted in considerably higher deviation from experimental point.
Since HDS and HDA reactions occur simultaneously and in parallel, in any hydrotreater reactor it is very important to monitor both sulphur and aromatics conversion, as well as temperature increase in the reactor (impacted by both exothermic reactions). Overall sulphur conversion along the catalyst bed is shown in Figure 2. Almost all of the applied HDA kinetic models predicted correctly the reduction of sulphur content in the reaction mixture along the catalyst bed. It can be concluded that all HDA kinetic models except those proposed by Liu and Avraam match very well with the experimental data for sulphur conversion (Figure 2a).
The temperature increase within the catalyst bed is shown in Figure 2b. Since the reactions occurring in the reactor (HDA and HDS) are exothermic, a temperature increase is to be expected. The degree of the temperature increase is very dependent on the applied kinetic model. The outlet temperature ranged from 611.7 K to 627.1 K and the best agreement with the experimental data was achieved using the Owusu-Boakye model.
The temperature increase and deviation from the experimental value was the highest for the Avraam and Yui kinetic models, which correspond to the highest proposed value of reaction heat for HDA reaction heat within the kinetic model. The calculated temperature increase was lowest for the de Oliveira kinetic model. Temperature change in the reactor had crucial impact on conversion and also on vapour-liquid equilibrium since equilibrium constants and liquid-to-vapour flow ratios are polynomial functions dependent on the tem-perature. At lower temperatures, more liquid is present in the reactor; as the temperature increases through the catalyst bed, liquid evaporates until catalyst particles are dried out. The temperature increase within the catalyst bed is shown in Figure 2b. Since the reactions occurring in the reactor (HDA and HDS) are exothermic, a temperature increase is to be expected. The degree of the temperature increase is very dependent on the applied kinetic model. The outlet temperature ranged from 611.7 K to 627.1 K and the best agreement with the experimental data was achieved using the Owusu-Boakye model.
The temperature increase and deviation from the experimental value was the highest for the Avraam and Yui kinetic models, which correspond to the highest proposed value of reaction heat for HDA reaction heat within the kinetic model. The calculated temperature increase was lowest for the de Oliveira kinetic model. Temperature change in the reactor had crucial impact on conversion and also on vapour-liquid equilibrium since equilibrium constants and liquid-to-vapour flow ratios are polynomial functions dependent on the temperature. At lower temperatures, more liquid is present in the reactor; as the temperature increases through the catalyst bed, liquid evaporates until catalyst particles are dried out.
The difference in catalyst wetting efficiency between the predictions of the models, using different kinetic equations for the HDA reaction, is very significant (Figure 3). Wetting of the catalyst bed is calculated depending on the reactor design and also the catalyst particle shape and size, but mostly on the liquid flow rate and velocity. The degree of liquid evaporation is dependent on the degree of temperature increase, while temperature increase is dependent on reaction rate and heat released by chemical reactions. The lowest wetting of the catalyst bed was predicted by simulations with the Yui and Avraam HDA models, where liquid evaporates almost instantly. For these models the temperature increase was the highest. The de Oliveira model predicted the lowest temperature increase (Figure 2), therefore wetting efficiency was very high with a substantial quantity of liquid through the entire reactor length. The difference in catalyst wetting efficiency between the predictions of the models, using different kinetic equations for the HDA reaction, is very significant (Figure 3). Wetting of the catalyst bed is calculated depending on the reactor design and also the catalyst particle shape and size, but mostly on the liquid flow rate and velocity. The degree of liquid evaporation is dependent on the degree of temperature increase, while temperature increase is dependent on reaction rate and heat released by chemical reactions. The lowest wetting of the catalyst bed was predicted by simulations with the Yui and Avraam HDA models, where liquid evaporates almost instantly. For these models the temperature increase was the highest. The de Oliveira model predicted the lowest temperature increase (Figure 2), therefore wetting efficiency was very high with a substantial quantity of liquid through the entire reactor length. Experimental data for the wetting efficiency in the simulated industrial reactor during the industrial test run were not available. It was formerly reported that very poor wetting efficiency is to be expected in the industrial-scale reactors that operate at low Experimental data for the wetting efficiency in the simulated industrial reactor during the industrial test run were not available. It was formerly reported that very poor wetting efficiency is to be expected in the industrial-scale reactors that operate at low LHSV [29]. A different study showed that the degree of the wetting efficiency of the catalyst bed affects the reaction rates greatly [30] and decreasing the amount of the liquid phase increases reaction rates and an intensification of the whole process can be achieved [31][32][33][34][35].
The thermodynamic model predicted catalyst bed wetting efficiency in the reactor, taking into account the flows of vapour and liquid phase, which are calculated through VLE calculations based on the PR-EOS. At the beginning of the process, the catalyst bed was partially wetted, while along the catalyst bed the evaporation of the liquid phase occurred. As this model is based on the inlet mixture that contains considerable amounts of volatile compounds, this phenomenon is likely and expected to occur (high methane content in the inlet hydrogen stream).
Overall effectiveness factors in the reactor for all models in vapour and in liquid phases are shown in Figure 4. By calculation of effectiveness factors, the impact of mass transfer limitations of the process is included in the model.  Overall effectiveness factors decreased along the reactor due to the temperature increase, and therefore reaction rates increased, resulting in a growing influence of the mass transfer resistance. In liquid phase, values of effectiveness factors were approaching zero due to the diminishing amounts of liquid in the catalyst bed, while in the gas phase effectiveness factors were higher. The values of calculated effectiveness factors are in agreement with the literature data [28,29].
Experimental data obtained in the industrial reactor were: total sulphur conversion of 99.49 wt%, conversion of aromatic compounds of 70 wt%, and temperature increase in the reactor of 11.5 K. Shown in Table 5 are the results of the total sulphur conversion, temperature increase in the reactor, and conversion of aromatics obtained by the six models and the comparison with experimental results.  Overall effectiveness factors decreased along the reactor due to the temperature increase, and therefore reaction rates increased, resulting in a growing influence of the mass transfer resistance. In liquid phase, values of effectiveness factors were approaching zero due to the diminishing amounts of liquid in the catalyst bed, while in the gas phase effectiveness factors were higher. The values of calculated effectiveness factors are in agreement with the literature data [28,29].
Experimental data obtained in the industrial reactor were: total sulphur conversion of 99.49 wt%, conversion of aromatic compounds of 70 wt%, and temperature increase in the reactor of 11.5 K. Shown in Table 5 are the results of the total sulphur conversion, temperature increase in the reactor, and conversion of aromatics obtained by the six models and the comparison with experimental results.
Values of overall sulphur conversions in Table 5 are slightly higher or lower than the experimental values. Model 6 calculated the lowest and Model 4 the highest values of sulphur conversion. Apart from the reaction kinetics, overall sulphur conversions were following the values of effectiveness factors, which were applied in material balance equations of each model, for calculation of the amounts of all specific sulphur and aromatic compounds through the reactor (Figures 3 and 4). All kinetic models except the Owusu-Boakye (Model 3) and de Oliveira (Model 5) models predicted poorly the conversion of aromatics in the simulations of the industrialscale reactor. The enormous deviation of outlet temperature in Model 4 occurred due to the very high value of the suggested HDA reaction heat, which was applied in the energy balance equation of the model. In addition, a large deviation from the temperature experimental data was observed for the Avraam and Liu models, which again proposed very high values for heat of the reaction.
Keeping in mind the experimental data, Model 3 could describe this catalytic reactor with the highest accuracy, or lowest deviation from the experimental data, among the investigated models. Model 3 described the system quite well, with very good predictions of sulphur and aromatics conversion, as well as the reactor temperature increase. Model 5 predicted well the conversion of aromatics, but it underestimated the outlet temperature and overestimated the sulphur conversion.
Novel research focused on blending different pyrolysis oils with light cycle oil [36,37] could be interesting to apply to this model for evaluating such technology in future research.

Conclusions
In this study, six different kinetic models for hydrodearomatisation reactions were compared through simulations of the industrial catalytic reactor for hydrotreating of straight run gas oil-containing FCC naphtha and light cycle oil fractions. The influence of hydrodearomatisation reaction kinetics on the conversion of sulphur, temperature increase, and wetting efficiency in the reactor were investigated. Experimental data from the industrial test run was used for validation of different models. Based on the obtained simulation results, it could be concluded that reactor performance is strongly dependent on the hydrodearomatisation reaction with respect to aromatics conversion but even more so with respect to the temperature increase in the reactor, which affects all key catalytic reaction parameters and thus also the sulphur conversion. Calculated total conversions of aromatic compounds differed considerably between models used, depending on the type of hydrodearomatisation kinetics total aromatic conversion ranges from to 13% to 63%. Different hydrodearomatisation kinetics do not affect conversion of the sulphur significantly when expressed as the percentage deviation; maximal absolute deviation from the experimental point was 1.2% conversion, but this relatively low deviation is critical for the quality of prediction since real conversions are expected within the narrow range of around 99.5 +/− 0.5%. Depending on the used kinetic model, the outlet temperature ranged from 607.7 K to 625 K, and the best prediction of experimental temperature was achieved using the Owusu-Boakye model, while the Avraam and Yui models overestimated outlet temperature due to very high values of the reaction heat proposed by these models. Catalyst wetting efficiency is an important parameter in the reactor and affects the reaction rates and overall efficiency of the process, but it cannot be evidenced directly in the industrial reactor. Good predictions of total sulphur and aromatics conversions, along with accurate prediction of temperature increases in the reactor, confirmed the relevance of overall efficiency and phase equilibrium calculations within the model. The best predictor of outlet temperature, sulphur, and aromatics conversion was achieved with the Langmuir-Hinshelwood kinetic model proposed by Owusu-Boakye.

Conflicts of Interest:
The authors declare no conflict of interest. Thiele modulus ε porosity of the catalyst particle ε b porosity of the catalyst bed ρ p density of the catalyst pellet, kg/m 3