Hydrocracking is a catalytic process in which the hydrocarbon molecules with longer chains break into lighter hydrocarbons with shorter chains. Hydrocracking units in the petroleum refineries usually feed with the heavy fractions of residual oil that are not commercially valuable and instead produces more valuable lighter fractions of hydrocarbons like gasoline and diesel. The hydrocracking process is especially important due to maximizing the use of crude oil as its resources are substantially reducing. Turning to use unconventional oil and gas reservoirs in recent years is an evidence of this importance    . Hydrocracking is a catalytic reaction at high-temperature and high-pressure conditions in the presence of hydrogen molecules  . Due to coke formation, poison deposition, and solid state transformation, catalyst deactivation occurs during the process lifetime causing a reduction in the production yield. To compensate the catalyst deactivation’s effects during the catalyst lifetime, reactor operating temperature should be adjusted (mainly increase) to maintain the desired production yield  . In this point of view, simulation is an essential tool to adjust the operating temperature while keeping the production’s yields in an acceptable range. Process simulation is also necessary for reactor design to achieve selective intermediate distillate production  .
Hydrocarbons cover a broad range of molecular weights with different types of organic compounds that makes it difficult to monitor the consumption/ production rates in the mixture of reactions in the hydrocracking process. Therefore, process engineers usually lump a group of organic compounds regardless of their molecular shape but based on their true boiling point (TBP) temperature having similar physio-chemical properties. Accordingly, kinetic modeling of heavy oil hydrocracking can be done based on the discrete or continuous lumping approach  while each lump is indexed with its TBP temperature range.
The discrete lumps kinetic model considers each lump characterized by its TBP temperature as one reactive species     . For example, Sedighi et al.  studied a 6-lump kinetic model by considering the vacuum gas oil (VGO) having TBP temperature greater than 380˚C, diesel (260˚C - 380˚C), kerosene (150˚C - 260˚C), heavy naphtha (90˚C - 150˚C), light naphtha (40˚C - 90˚C), and gases (<40˚C) as discrete lumps of the process. In the discrete lumps approach, the heavier products associated with a higher TBP temperature can partially convert to the lighter products in a parallel/series reaction chain while each reaction path has its own reaction kinetic rate. Therefore, the number of model parameters is directly proportional to the number of accounted lumps. When the number of model parameters is higher, the parameter estimation algorithm needs a higher amount of measured data in order to accurately estimate the parameter values. Balasubramanian et al.  considered a 5-lump model based on carbon number and ended up to estimate 40 unknown parameters for their model.
The continuous lumping kinetic model     can overcome this drawback where the number of parameters is independent of the number of lumped fractions in the model. In the continuous model approach, a distribution of reaction products over a range of TBP temperatures represents the model’s output. Integration over a discrete boiling point represents the concentration of desired lump fractions  . Therefore, the continuous model has an advantage over the discrete lumped model, because any fraction of lumped products can be calculated from the same continuous distribution curve.
Both the discrete and the continuous kinetic models have unknown parameters that need to be estimated in order to calibrate the model. Regression techniques have been successfully used to minimize the misfit between the modeled and measured data. Sadeghi et al.  and Elizalde et al.  applied the continuous lumping model over different sets of measured data to minimize the least square error between the modeled and measured points and obtained a point estimate of the model parameters. Kumar et al.  applied hybrid particle swarm optimization to estimate the continuous lumping parameter values.
However, there are some uncertainties associated with the value of the parameters that are not considered in the point estimation methods. The measurement errors, the model structural error (due to model simplifications), and errors associated with the operating conditions (like isotherm reactor assumption) are the main sources of uncertainty in the hydrocracking kinetic models that can affect the value of the estimated parameters. To address the uncertainties in the parameter estimation, a probabilistic approach can be considered. Different techniques have been applied in scientific fields to deal with uncertainties. Commonly, Monte Carlo simulations over a pre-assumed range of parameters    , the autoregressive moving average for uncertainty analysis associated with the time series    , the generalized likelihood uncertainty estimation (GLUE) method  , and the Bayesian approach     are used uncertainty quantification. Albrecht  studied four types of reaction models (with increasing complexity) including linearized second-order reaction, single elementary reaction, single elementary reaction coupled with temperature dependency, and catalytic reaction cycle aiming at evaluating the confidence interval associated with each model parameters. He applied different techniques to regress the model parameters to experimental observations with artificially added noise; and concluded that for highly non-linear models, the Markov Chain Monte Carlo (MCMC) algorithms utilizing a Bayesian approach accurately estimates uncertainty.
Parameter estimation and uncertainty analysis by using Bayesian approach are widely used in different fields of science. Alikhani et al.  applied Bayesian inference to estimate the confidence interval of groundwater residence time distributions by using multiple groundwater age tracers as observed data points. Alikhani et al.  evaluated the information content of long-term measured data from a municipal wastewater treatment plant to evaluate the confidence interval of activated sludge model parameters. They assessed the level of obtained information associated with the measured data points by comparing the entropy of the posterior distribution of parameters to their prior distributions. They concluded that in multi-dimensional parameter estimation, the level of information obtained from a set of measured data is different for each parameter and therefore suggested to perform sensitivity analysis to select a subset of parameters to be estimated. Sun et al.  applied the Bayesian approach to estimate the isomer kinetic decay rate in the phytate (IP6) degradation pathway. Their reaction modeling network is conceptually similar to the discrete lumps kinetic models where higher order molecules in a hierarchical degradation pathway convert to lower order molecules. Therefore, the Bayesian inference approach is assumed to be suitable for this study to obtain the confidence interval of hydrocracking model parameters.
Nonetheless, Bayesian parameter estimation framework needs to be solved by applying MCMC algorithms that usually requires a relatively larger number of model simulations. Therefore, an efficient numerical algorithm  should be selected to reduce the overall computation time of each simulation run.
In this study, the continuous lumping kinetic model is slightly modified to take into consideration the temperature dependency of parameters. The Bayesian parameter estimation approach is applied to measured data obtained from a hydrocracking unit of a petroleum refinery to obtain the posterior credible intervals of the model parameters. Finally, the Monte Carlo simulation is performed by taking samples over the posterior range of estimated parameters to evaluate the confidence interval of the model output concentrations in different operating conditions.
2. Materials and Methods
2.1. Description of the Model
The continuous lumping kinetic model in this study is explained in detail in  and  based on the original model presented in Laxminarasimhan et al.  . The model is defined based on a dimensionless temperature (θ) valued between 0 to 1 showing the range of TBPs between the lowest () to highest () TBP values in the hydrocracking process:
The first-order reactivity (k) of each component can be related to its by:
where is the reactivity of the component with the highest TBP (), and is the model parameter in the operating temperature T. The mass balance equation is then can be obtained as:
where the right-hand side shows the consumption (the first term) and production (the second term) rate of the as the concentration of the component with reactivity k. is the species-type distribution function that transfers the system with N discrete components to a continuous lumping space with the following equation:
The model is called continues because a cumulative distribution is assigned to the fraction yield of the species. represents the production yield of component with reactivity k from cracking of component with reactivity K and is given by:
where A and B are defined as:
where, , and are the model parameters. is the condition that satisfies the, and can be obtained as:
Solving the continuous model results in a distribution of the component concentrations. To obtain the weight fraction of each discrete fraction, the following integration can be performed:
where is the concentration of a specific fraction between the and.
To take into consideration the temperature effect on the model parameters, an Arrhenius-type relationship  is adopted into the model parameters as:
where and are the parameter values at temperature and refetrence temperature, respectively. is the temperature dependency coefficient for any model parameter. By this modification, the model can be applied to a broader range of operating temperatures with a unique set of parameter values. The continuous kinetic model is theoretically applicable for a broad range of hydrocracking products and operating conditions if the model outputs are in a good agreement with the experimental results.
2.2. Bayesian Parameter Estimation
The main idea of the Bayesian inference application in the parameter estimation is extracted from the work by Alikhani et al.  . In this approach, a prior distribution for each parameter should be assigned first. The prior range can be obtained by using the values reported in other studies for the similar operating conditions. The confidence about the prior range can be improved by applying the Bayes’ theorem given the measured data points. The Bayes’ theorem can be shown as:
The results of applying Bayes’ theorem would be obtaining the posterior distribution of model parameters given a measured data set of; while represents the prior distribution. A normal distribution is considered for all the parameters in this study.
In Equation (2), represents the likelihood function. By assuming Gaussian distribution of errors, the likelihood function can be shown as:
where n is the total number of measured data points, is a single measured point while is its corresponding model output for a certain parameter set of. is the likelihood standard deviation, representing the quantified uncertainty associated with the measured data points. represents the five kinetic parameters (, , , , and) at reference temperature plus five temperature dependency coefficients (, , , , and).
The Bayesian parameter estimation is particularly useful for the system of unknown parameters that are highly dependent on the operating conditions in which the calibrated parameter value of other studies is not suitable to be used in another study. Nonetheless, the information about the parameter values in other studies can be used to construct the prior distributions; and Bayesian approach can extract the information in the observed data set to enhance our confidence about the parameter values  . To obtain the posterior distribution by using the Bayesian approach, the mechanistic model needs to be solved by using Markov Chain techniques in a random walk approach  . Metropolis-Hasting algorithm is applied in this study to sample parameter sets of from posterior distribution  .
2.3. Measured Data Points
Measured data points obtained from the case study introduced in  . The measured data were collected from a hydrocracking unit of a petroleum refinery where VGO feeds in 395˚C and 20 Mpa into 2 reactors in series. Reactors are filled with a NiO-WO3/SiO2-Al2O3 catalyst having the density of 790 kg/m3, specific surface area of 235 m2/g, and porosity of 0.5. The hydrocracking products were collected (Table 1) in the form of light petroleum gas, LPG, (<39˚C), naphtha (39˚C - 150˚C), kerosene (150˚C - 250˚C), diesel (250˚C - 380˚C), and VGO (>380˚C). The measured data points at four different reactor temperatures (390˚C, 410˚C, 430˚C, and 450˚C) and six different residence times are shown in Table 1. Part of the measured data points is also shown in Figure 1 and Figure 2.
3. Results and Discussion
Genetic algorithm (GA) used to obtain the point estimates of the modified continuous kinetic model parameters by maximizing the likelihood function (Equation (12)). In total, 10 model parameters plus the likelihood standard deviation treated as unknown parameters to be estimated by applying the observed measured data. The results from GA for two operating temperatures at 390˚C and 450˚C are shown in Figure 1. The results from GA optimization indicate that the continuous lumping model is reasonably able to capture the trend and magnitude of the measured data. Moreover, results in Figure 1 show that applying the temperature dependency of parameters into the kinetic model works very well for the same set of parameter values at all the four reactor temperatures. In this study, 390˚C is considered as the reference temperature and the five main parameters of the continuous model (, , , , and) are estimated in this temperature. The observed data points at the other temperatures are applied to estimate the temperature dependency values of the parameters.
The parameter values obtained by GA optimization were used as starting
Table 1. Measured data points of hydrocracking process used in this study  .
Figure 1. Modeled (lines) vs. measured (symbols) for the reactor temperature at left) 390˚C, and, right) 450˚C. The GA’s point estimate of the model parameters is used to obtain the modeled data.
Figure 2. 95% model output confidence interval (floating bars) vs. measured data (circles) at four different operating temperatures with three different residence times.
points of 10 parallel Markov chains aiming at sampling 100,000 posterior set of parameter values in the Bayesian parameter estimation approach. The first 10,000 of total samples discarded due to burn-in period  . Statistical analysis performed on the remaining 90,000 samples to obtain the 95% credible intervals of the parameters.
Values of 2.5th, 50th and 97.5th percentile plus the standard deviation of each model parameter are shown in Table 2. The negative values of and show that the value of and decrease when operating temperature increase.
Table 2. 95% Bayesian credible interval and the standard deviation of the modified continuous lumping kinetic model parameters.
This finding is in good agreement with the parameter-temperature correlation presented in Elizalde et al.  . The small value for the likelihood standard deviation () shows that the level of uncertainty in the presented system is low; indicating that the errors associated with the measured data and model structural errors are low.
In Table 2, the posterior range of the parameters can dynamically be reevaluated by introducing a new measured data set. In this way, the current posterior range can be introduced as a given prior range to the Bayesian approach and the modified posterior range can be obtained by applying new observed data.
To evaluate the ability of the model to meet the observed data, the Monte Carlo simulation is performed and 5000 realization parameter sets randomly sampled from the 95% posterior range. The model outputs were statistically analyzed and the 95% confidence interval of weight fraction of hydrocracking products are obtained and illustrated in Figure 2.
Results in Figure 2 show that most of the measured data points are fall inside the 95% model output range (shown as floating bars). The results also show that for the lightest fraction (LPG), continuous lumping model weakly predicted the observed data. Nonetheless, the model is predicting the moderate and heavy fractions with a good agreement.
The model confidence intervals reflect the uncertainty level associated with the parameter values and can be used in decision-making step, obtaining factors of safety in designing step, and increasing the level of accuracy in the data measuring (and sampling) step; and generally, can enhance the modeling and simulation efficiency.
In this study, the parameter estimation of the hydrocracking process model is assessed. One system consists of five fractions (LPG, naphtha, kerosene, diesel, and VGO) is modeled by using continuous lumping approach. The continuous model is modified by considering the effect of reactor temperature in the parameter values.
The Genetic Algorithm and the Bayesian parameter estimation approach are applied to obtain the point estimate and the credible interval of the model parameters. A good agreement between the calibrated model output and the measured data is observed showing that the continuous lumping model is able to simulate the presented hydrocracking process. The results also show that the modified continuous lumping model considering the temperature dependency of the parameters extends the ability of the model on different operating temperature. Applying the Bayesian approach resulted in the 95% credible interval of model parameters reflecting the uncertainty associated with parameter values.
A probabilistic simulation is also performed by using the posterior range of the parameters to obtain the confidence interval of model outputs. The results show that for the studied process unit, the uncertainty associated with the measured data and the model structural error is quantitatively low. The Bayesian approach based on the Markov Chain Monte Carlo simulation is shown to be efficient to quantify the uncertainty associated with the parameter values of the continuous lumping model.
The following general conclusions obtained from the presented study:
1) Continuous lumping kinetic model was able to simulate the hydrocracking of VGOs in the range of 390˚C to 450˚C, and the residence time of up to 2 hr.
2) The model parameters were estimated by maximizing the likelihood function between the model outputs and measured data points.
3) The temperature dependency of the model parameters was successfully embedded into the continuous lumping model and the temperature dependency coefficients were estimated.
4) The uncertainty associated with the parameter values was evaluated by applying Bayesian theorem, and MCMC technique and the posterior range of parameters were obtained.
5) Monte Carlo simulations were performed to evaluate the confidence interval of hydrocracking products’ yield.
 Boosari, S.S.H., Aybar, U. and Eshkalak, M.O. (2015) Carbon Dioxide Storage and Sequestration in Unconventional Shale Reservoirs. Journal of Geoscience and Environment Protection, 3, 7-15.
 Kazemi, M. and Takbiri-Borujeni, A. (2016) Flow of Gases in Slit Shaped Organic Nanopores of Shale: A Boundary-Driven Molecular Simulation Study. SPE Western Regional Meeting, Anchorage, 23-26 May 2016, SPE-180441-MS.
 Barkhordari, A., Fatemi, S. and Daneshpayeh, M. (2009) Kinetic Modeling of Industrial VGO Hydrocracking in a Life Term of Catalyst. Proceedings of the 8th World Congress of Chemical Engineering (WCCE8), city, date, page.
 Elizalde, I., Rodríguez, M.A. and Ancheyta, J. (2009) Application of Continuous Kinetic Lumping Modeling to Moderate Hydrocracking of Heavy Oil. Applied Catalysis A: General, 365, 237-242.
 Ali, M.A., Tatsumi, T. and Masuda, T. (2002) Development of Heavy Oil Hydrocracking Catalysts Using Amorphous Silica-Alumina and Zeolites as Catalyst Supports. Applied Catalysis A: General, 233, 77-90.
 Sadighi, S., Ahmad, A., Mohaddecy, S. and Reza, S. (2010) 6-Lump Kinetic Model for a Commercial Vacuum Gas Oil Hydrocracker. International Journal of Chemical Reactor Engineering, 8, p
 Balasubramanian, P. and Pushpavanam, S. (2008) Model Discrimination in Hydrocracking of Vacuum Gas Oil Using Discrete Lumped Kinetics. Fuel, 87, 1660-1672.
 Aris, R. and Gavalas, G.R. (1966) On the Theory of Reactions in Continuous Mixtures. Philosophical Transactions of the Royal Society of London A: Mathematical, Physical and Engineering Sciences, 260, 351-393.
 Sin, G., Gernaey, K.V., Neumann, M.B., van Loosdrecht, M.C. and Gujer, W. (2009) Uncertainty Analysis in WWTP Model Applications: A Critical Discussion Using an Example from Design. Water Research, 43, 2894-2906.
 Neurock, M., Libanati, C., Nigam, A. and Klein, M.T. (1990) Monte Carlo Simulation of Complex Reaction Systems: Molecular Structure and Reactivity in Modelling Heavy Oils. Chemical Engineering Science, 45, 2083-2088.
 Alikhani, J., Stewart, H.A., Takacs, I., Omari, A.A., Murthy, S. and Massoudieh, A. (2015) Time Series Analysis for Uncertainty Evaluation of Influent Fluctuation of a Wastewater Treatment Plant. Proceedings of the Water Environment Federation, 2015, 3519-3526.
 Constantinescu, E.M., Chai, T., Sandu, A. and Carmichael, G.R. (2007) Autoregressive Models of Background Errors for Chemical Data Assimilation. Journal of Geophysical Research: Atmospheres, 112, p
 Beran, J. (1995) Maximum Likelihood Estimation of the Differencing Parameter for Invertible Short and Long Memory Autoregressive Integrated Moving Average Models. Journal of the Royal Statistical Society. Series B (Methodological), 57, 659-672.
 Sathyamoorthy, S., Vogel, R.M., Chapra, S.C. and Ramsburg, C.A. (2014) Uncertainty and Sensitivity Analyses Using GLUE When Modeling Inhibition and Pharmaceutical Cometabolism During Nitrification. Environmental Modelling & Software, 60, 219-227.
 Alikhani, J., Takacs, I., Omari, A. A., Murthy, S., Mokhayeri, Y. and Massoudieh, A. (2014) Inverse Modeling of Nitrification-Denitrification Processes: A Case Study on the Blue Plains Wastewater Treatment Plant in Washington, DC. Proceedings of the Water Environment Federation, 19, 4556-4567.
 Prager, J., Najm, H.N., Sargsyan, K., Safta, C. and Pitz, W.J. (2013) Uncertainty Quantification of Reaction Mechanisms Accounting for Correlations Introduced by Rate Rules and Fitted Arrhenius Parameters. Combustion and Flame, 160, 1583-1593.
 Sohn, M.D., Small, M.J. and Pantazidou, M. (2000) Reducing Uncertainty in Site Characterization Using Bayes Monte Carlo Methods. Journal of Environmental Engineering, 126, 893-902.
 Alikhani, J., Deinhart, A.L., Visser, A., Bibby, R.K., Purtschert, R., Moran, J.E., Massoudieh, A. and Esser, B.K. (2016) Nitrate Vulnerability Projections from Bayesian Inference of Multiple Groundwater Age Tracers. Journal of Hydrology, 543, 167-181.
 Alikhani, J., Takacs, I., Al Omari, A., Murthy, S. and Massoudieh, A. (2017) Evaluation of the Information Content of Long-Term Wastewater Characteristics Data in Relation to Activated Sludge Model Parameters. Water Science and Technology, 75, 1370-1389.
 Sun, M., Alikhani, J., Massoudieh, A., Greiner, R. and Jaisi, D.P. (2017) Phytate Degradation by Different Phosphohydrolase Enzymes: Contrasting Kinetics, Decay Rates, Pathways, and Isotope Effects. Soil Science Society of America Journal, 81, 61-75.
 Alikhani, J., Shoghli, B., Bhowmik, U.K. and Massoudieh, A. (2016) An Adaptive Time-Step Backward Differentiation Algorithm to Solve Stiff Ordinary Differential Equations: Application to Solve Activated Sludge Models. American Journal of Computational Mathematics, 6, 298-312.
 Alikhani, J., Al-Omari, A., De Clippeleir, H., Murthy, S., Takacs, I. and Massoudieh, A. (2017) Assessment of the Endogenous Respiration Rate and the Observed Biomass Yield for Methanol-Fed Denitrifying Bacteria under Anoxic and Aerobic Conditions. Water Science and Technology, 75, 48-56.