A Review of Kinetic Modeling Methodologies for Complex Processes
Une revue de méthodes de modélisation cinétique pour des procédés complexes
IFP Energies nouvelles, Rondpoint de l’échangeur de Solaize, BP 3, 69360
Solaize – France
email: Jan.VERSTRAETE@ifpen.fr
^{*} Corresponding author
In this paper, kinetic modeling techniques for complex chemical processes are reviewed. After a brief historical overview of chemical kinetics, an overview is given of the theoretical background of kinetic modeling of elementary steps and of multistep reactions. Classic lumping techniques are introduced and analyzed. Two examples of lumped kinetic models (atmospheric gasoil hydrotreating and residue hydroprocessing) developed at IFP Energies nouvelles (IFPEN) are presented. The largest part of this review describes advanced kinetic modeling strategies, in which the molecular detail is retained, i.e. the reactions are represented between molecules or even subdivided into elementary steps. To be able to retain this molecular level throughout the kinetic model and the reactor simulations, several hurdles have to be cleared first: (i) the feedstock needs to be described in terms of molecules, (ii) large reaction networks need to be automatically generated, and (iii) a large number of rate equations with their rate parameters need to be derived. For these three obstacles, molecular reconstruction techniques, deterministic or stochastic network generation programs, and singleevent microkinetics and/or linear free energy relationships have been applied at IFPEN, as illustrated by several examples of kinetic models for industrial refining processes.
Résumé
Dans cet article, les techniques de modélisation cinétique des processus chimiques complexes sont examinées. Après un bref aperçu historique de la cinétique chimique, un aperçu des bases théoriques de la modélisation cinétique d’étapes élémentaires et de réactions globales est présenté. Les techniques classiques de regroupement (lumping) sont ensuite présentées et analysées. Deux exemples de modèles cinétiques regroupés (pour l’hydrotraitement de gazole atmosphérique et pour l’hydrotraitement de résidus) développés à IFP Energies nouvelles (IFPEN) sont présentés. La plus grande partie de cette revue décrit des stratégies avancées de modélisation cinétique, dans lesquelles le détail moléculaire est retenu : les réactions entre les molécules sont représentées ou même subdivisées en étapes élémentaires. Pour être en mesure de conserver ce niveau moléculaire à la fois dans le modèle cinétique et dans les simulations de réacteurs, plusieurs obstacles doivent d’abord être éliminés : (i) la charge doit être décrite en termes de molécules, (ii) les grands réseaux réactionnels doivent être générés automatiquement et (iii) un grand nombre d’équations de vitesse avec leurs paramètres de vitesse doit être dérivé. Pour ces trois obstacles, des techniques de reconstruction moléculaire, des programmes de génération de réseaux déterministes ou stochastiques, et des modèles microcinétiques basés sur des événements constitutifs (single events) et/ou des relations linéaires d’énergie libre (Linear Free Energy Relationships) ont été utilisés à IFPEN, comme illustré par plusieurs exemples de modèles cinétiques pour des procédés de raffinage industriels.
© L.P. de Oliveira et al., published by IFP Energies nouvelles, 2016
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
List of Abbreviations
AEBP: Atmospheric Equivalent Boiling Point
ASTM: American Society for Testing Materials
DFRN: Deterministic Full Reaction Network
DFT: Density Functional Theory
DIPPR: Design Institute for Physical PRoperties
DLRN: Deterministic Limited Reaction Network
FID: Flame Ionization Detector
FTICRMS: Fourier Transform – Ion Cyclotron Resonance – Mass Spectrometry
HPLC: HighPerformance Liquid Chromatography
IUPAC: International Union of Pure and Applied Chemistry
LFER: Linear FreeEnergy Relationships
LHSV: Liquid Hourly Space Velocity
LSODE: Livermore Solver for Ordinary Differential Equations
MARI: Most Abundant Reaction Intermediate
MASI: Most Abundant Surface Intermediate
MTHS: Molecular Type and Homologous Series approach
NCD: Nitrogen Chemiluminescence Detector
NIST: National Institute of Standards and Technology
NMR: Nuclear Magnetic Resonance
ODE: Ordinary Differential Equation
PDE: Partial Differential Equation
PDF: Probability Distribution Function
QFER: Quadratic FreeEnergy Relationship
QSAR: Quantitative StructureActivity Relationships
QSRC: Quantitative StructureReactivity Correlations
REM: Reconstruction by Entropy Maximization
SARA: Saturates/Aromatics/Resins/Asphaltenes
SCD: Sulfur Chemiluminescence Detector
SEC: Size Exclusion Chromatography
SLRN: Stochastically Limited Reaction Network
SOL: Structure Oriented Lumping
SSA: Stochastic Simulation Algorithm
TRC: Thermodynamics Research Center
List of Symbols
A_{i} : Preexponential factor for the rate constant of reaction i
c_{i}: Stochastic rate constant for reaction i
C_{j}: Concentration of species j
E : Shannon’s information entropy criterion
E_{i}: Activation energy for reaction i
F_{k}: Mixing rule for analytical property k
h : Planck’s constant (6.626068 × 10^{−34} m^{2} kg s^{−1})
: Intrinsic singleevent rate coefficient of the elementary step
k_{i}: Rate constant for reaction i
k_{B}: Boltzmann’s constant (1.3806503 × 10^{−23} m^{2} kg s^{−2} K^{−1})
K : Number of analytical properties or constraints
K_{c}: Equilibrium constant (concentration basis)
K_{eq}: Overall equilibrium constant
K_{i}: Equilibrium constant for reaction / elementary step i
M_{j}: Molecular weight of component j
M_{AB}: Reduced molecular weight of components A and B (collision theory)
n_{e}: Number of single events
N_{A}: Avogadro constant (6.0221408577 × 10^{23} mol^{−1})
N_{Alkyl}: Number of alkyl substituents in β position of a sulfur atom
N_{AR}: Number of aromatic rings
N_{SR}: Number of saturated rings
N_{TR}: Number of thiophene rings
P_{k}: Value of analytical property k
r_{i}: Reaction rate of reaction i
R : Universal gas constant (8.3144 J mol^{−1} K^{−1})
R_{j}: Rate of production of component j
s : Number of elementary steps in a reaction mechanism
x_{j}: Mole fraction of component j
Greek symbols
: Proportionality constants (in LFER, QFER, QSAR, QSRC equations)
ΔG_{R}: Gibbs free energy of reaction
ΔG^{≠}: Activation Gibbs free energy (Transition State Theory)
ΔH^{≠}: Activation enthalpy (Transition State Theory)
Δn : Difference in molecularity between the reverse and the forward reaction
Δn^{≠}: 1 – molecularity of the reaction that forms the activated complex (TST)
ΔS^{≠}: Activation entropy (Transition State Theory)
λ : Reorganization energy (Quadratic FreeEnergy Relationship)
ν_{j}: Stoichiometric coefficient for component j
ν_{ij}: Stoichiometric coefficient for component j in reaction i
ρ_{H}: Hammett’s sensitivity constant
σ_{i}: Stoichiometric number for reaction/elementary step i
σ_{AB}: Reaction crosssection (collision theory)
σ_{H}: Hammett’s substitution constant
σ^{r}: Symmetry number for the reactant
σ^{≠}: Symmetry number for the activated complex
θ_{i}: Reactivity index of molecule i
θ_{ij}: Reactivity index j of molecule i
Subscripts
eq : At equilibrium conditions
glob : For global symmetry and chirality
overall : For the overall reaction
Superscripts
≠: Transition state or activated complex
Introduction
When dealing with chemical reactions, thermodynamics allows us to know whether a reaction is possible, and to calculate its chemical equilibrium constant. Hence, for a given chemical reaction under given experimental conditions and with a known initial composition of the mixture, thermodynamics is able to predict its equilibrium composition, the final steady state. However, thermodynamics does not say anything concerning the rate at which a reaction proceeds towards equilibrium, nor the pathway that is followed, i.e. the dynamics that lead to that final steady state.
The study of the dynamics of reactions is the field of chemical kinetics, namely the study of the factors that determine the reaction rates and that are responsible for the establishment of chemical equilibria in reversible reactions. The goal of chemical kinetics is therefore twofold:

determining the overall pathway, the reaction network, the reaction mechanisms and the various intermediates, and

deriving quantitative rate equations (with their rate parameters) which predict the rate of reaction as a function of the local conditions (temperature, pressure, phase state, composition, ionic strength, pH, solvent, catalysts, etc.).
The formulation of the reaction mechanism used for the kinetic modeling is always subject to discussion and should be based on all available physical measurements, including surface science data, chemical knowledge, and kinetic experiments (preferably far from equilibrium). Hence, chemical kinetics needs to combine knowledge acquired through various other disciplines, such as analytical chemistry (composition), organic chemistry (reaction types), physical chemistry (structure of matter), classical thermodynamics (equilibrium), statistical thermodynamics and quantum mechanics (rate theories), spectroscopy and computational chemistry (reaction intermediates), and mathematics (parameter estimation and process simulation), among others. By predicting the rates of the various reaction pathways, chemical kinetics allows the prediction of production rates and selectivities, and is therefore a necessary tool in the modeling and design of chemical reactors. Hence, chemical kinetics is one of the pillars of the chemical engineering discipline.
Chemical kinetics is a discipline that quantitatively describes the progress of reactions on a large range of different scales: from interactions between atoms and electrons in chemical bonds, to production rates in chemical reactors. Kinetics on the level of individual molecules describes the modifications that the reactants undergo to form the products. This is often referred to as reaction dynamics or molecular dynamics, and is nowadays an integral part in the domain of chemical kinetics. In these studies (molecular beam experiments, laser spectroscopy, ab initio theoretical chemistry, transition state theory, etc.), the molecular structures that exist along the reaction pathway on the Potential Energy Surface (PES) are analyzed, directly linking reactivity to quantum mechanics. Kinetics on large ensembles of molecules describes the average behavior of the reactants, linking reactivity and the extent of reaction to the operating conditions (temperature, composition, etc.) and thermodynamics.
Chemical kinetics is a relatively young science. Although descriptive kinetics already started just before 1800, with publications from Karl Friedrich Wenzel in 1777 and from Claude Louis Berthollet in 1803, the first known study of quantitative chemical kinetics dates back to 1850, when German physicist Ludwig Ferdinand Wilhelmy published a study on the acid hydrolysis of sucrose into glucose and fructose in a batch reactor, establishing for the first time an empirical relation between the reaction rate and the concentration of the reactants. In his rate equation, the reaction rate was proportional to the concentrations of sugar and acid, while his proportionality constant was an exponential function of the inverse of the absolute temperature. Moreover, he proposed a differential equation to describe the decrease in sugar concentration over time, thereby writing down a differential equation that was probably the first differential equation in chemistry and one that is still used today for elementary bimolecular reactions in closed systems.
In their 1862 publication on the esterification between acetic acid and ethanol, French chemists Pierre Eugène Marcellin Berthelot and Léon Péan de SaintGilles showed that some reactions cannot proceed to total conversion and they had to introduce into their empirical rate equation a limiting factor, which was later shown to be the equilibrium constant. By reexamining these results and interpreting their own experimental data through the laws of classical mechanics, two Norwegian scientists, mathematician Cato Maximilian Guldberg and chemist Peter Waage, made a significant breakthrough by formulating the mass action law, one of the basic laws in chemical kinetics. In their 1864 publication, they introduced the equilibrium mass action law, and formulated the form of the rate equations. They arrived at this form through the laws of classical mechanics: they postulated that molecular collisions have to precede the actual reaction, and proposed that reaction equilibrium is a balance of two opposing ‘affinity forces’, one due to the reactants and one owing to the products. For the esterification reaction, they stated that, at equilibrium, the two opposing ‘affinity forces’ (i.e. the reaction rates) had to be equal, and that each ‘affinity force’, in analogy to the law of gravity, could be expressed as the product of the ‘affinity coefficients’ (i.e. rate constants) and the ‘action masses’ (i.e. concentrations) of the two reacting species. In their later work, they generalized these mass action laws to any number of reactants, both for irreversible and reversible reactions of any molecularity, by explicitly using the stoichiometric coefficients of the reaction in the rate equation, leading to the publication of the mass action law in its final form in 1879.
In the same period, at the University of Oxford, chemist Augustus George Vernon Harcourt sought the help of mathematician William Esson. In their papers from 1865 and 1866, they demonstrated that the reactant concentrations decrease logarithmically with time. Therefore, the rate of reaction does not remain constant, but depends on the concentration of the reactant. They were also the first to treat the kinetics of consecutive reactions. In their 1867 paper, the time dependent evolution was interpreted in terms of differential equations. In their meticulous and systematic study of the concentration dependence of reaction rates, the reaction rate was assumed to be proportional to the original amounts of the reactants and a constant depending on the considered system, which they baptized rate constant k. In their studies on the temperature effect, they also postulated that this rate constant has to become zero at absolute zero, since the reactants no longer move and cannot interact or collide. This, together with their comparison of a chemical reaction to the fall of bodies, where potential energy is transformed into kinetic energy, illustrates that Harcourt and Esson started the introduction of thermodynamics into chemical kinetics.
The next major event arrived in 1884, when Dutch scientist Jacobus Hendricus van ‘t Hoff published the first monograph on chemical kinetics under the title ‘Études de dynamique chimique’. In his book, the author compiled, generalized and further developed the work on chemical kinetics. He described a graphical method for determining the order of a reaction in kinetic studies, studied the influence of the solvent for reactions in solution, introduced the modern concept of chemical affinity to replace the classical mechanical representation, classified reactions in terms of the number of interacting species (which is now referred to as their molecularity), derived kinetic rate equations for mono and bimolecular reactions, and considered polymolecular processes as a sequence of mono and/or bimolecular steps. He derived the van ‘t Hoff relation for the dependence of the equilibrium constant on the absolute temperature, and was also the first to correctly formulate the effect of temperature on the rate constants, although he did not succeed in finding a theoretical interpretation for the various terms. His essential contribution lies in the link he made between the equilibrium constant of a reaction, and the rate constants for the direct and the reverse reaction, thereby formally connecting thermodynamics and kinetics.
In Baltic Germany, Friedrich Wilhelm Ostwald worked on catalysis, electrochemistry, chemical equilibria (Ostwald’s dilution law) and reaction rates. In 1887, Ostwald derived the kinetics of autocatalytic reactions, introduced the terms ‘halflife’ and ‘reaction order’, and worked on ways to determine the reaction order. During the same period, he also coined the word ‘mole’, linked it to ideal gases, and defined it as the molecular weight of a substance in mass grams. Ostwald also wrote the famous twovolume textbook ‘Lehrbuch der Allgemeinen Chemie’, in which he included his developments on catalysis and chemical kinetics. His most important contributions consisted of formulating that each reaction exists and proceeds independently, interacting only on the level of a material balance (principle of independence of chemical reactions), that catalysis was nothing but a kinetic phenomenon, where a foreign compound accelerated the kinetics by lowering the energy barrier, and hence, that a catalyst cannot change the equilibrium of a reaction.
Based on the observation that reactions accelerate when heat energy is added, and on the formulae derived by van ‘t Hoff, Swedish scientist Svante August Arrhenius postulated in 1889 that, in order for a reaction to proceed, an energy barrier has to be overcome before two molecules will react. As did van ‘t Hoff previously, he claimed that the temperature was not the reason for reaction, but that it is responsible for the changes in its rate. In the same work, Arrhenius also had an important additional insight: he postulated that this temperature dependence indicates the existence of an ‘activated molecule’, whose concentration is proportional to the total concentration of reactant molecules, but is exponentially temperature dependent. To satisfy the latter conditions, he proposed an equilibrium between normal and activated reactant molecules. By introducing the concept of an ‘activated state’, Arrhenius can therefore be considered the father of all reaction rate theories. He then expressed the temperature dependence of the reaction rate through the rate constants in his famous Arrhenius equation (while crediting van ‘t Hoff for the origin of his equation), simultaneously coining the term ‘activation energy’ and introducing its symbol, E. Both the preexponential factor (or frequency factor) and the activation energy are considered to be independent of temperature. Arrhenius also provided a molecular level interpretation of the exponential factor as the fraction of collisions with sufficient energy, but his equation remains an empirical equation, however, as the preexponential factor is not explained at all.
In 1905, German physicist Walther Hermann Nernst established what is now called the Third law of thermodynamics (which describes the behavior of matter as the temperature approaches absolute zero, and postulates that the entropy of a perfect crystal at absolute zero is exactly equal to zero). He described the calculation of chemical affinities, provided for the first time a means of determining free energies and hence equilibrium states of chemical reactions from heat measurements, and allowed the unambiguous calculation of entropies.
In Germany, at the end of the 19^{th} century, Max Ernst August Bodenstein studied the bromination of hydrogen to produce hydrogen bromide in conditions where the reaction was not limited by equilibrium. Yet, he found that the rate of the reaction was limited by its product. Without knowing, Bodenstein had discovered a freeradical chain reaction with two propagation steps. Bodenstein also showed that Arrhenius’ law is only applicable to elementary reactions and that overall reactions often show deviations. Later on, in 1913, Bodenstein studied the chlorination of hydrogen and observed that absorption of a single photon could yield up to 10^{6} molecules of HCl, clearly indicating that this reaction could not be a single step reaction. This behavior in various halogenation reactions was explained in 1918 by Walther Hermann Nernst, who conceived the theory for chain reaction mechanisms that was further developed by Bodenstein. In his 1913 paper, Bodenstein also postulated that chemical reactions follow a reaction mechanism consisting of several steps, in which reaction intermediates are present in inferior amounts, that the rate of change of their concentrations can be considered negligible and that the reaction intermediates are therefore in a quasisteady state, an approach which had been proposed six months earlier in 1913 by David Leonard Chapman. Throughout his work and publications, Bodenstein developed, implemented and disseminated this quasisteady state approach, of which he was one of the biggest advocates.
In the same period, chemical kinetics was also extended to systems other than homogeneous gas phase reactions. In 1913, German biochemist Leonor Michaelis and Canadian physician Maud Leonora Menten investigated the kinetics of an enzymatic reaction, in which an invertase catalyzed the hydrolysis of sucrose into glucose and fructose. They proposed a twostep reaction mechanism, in which the first step was at quasiequilibrium, and developed the corresponding rate equation. Reactions on surfaces were studied in a large number of papers, as many elements were found to promote a large range of reactions. The understanding of the elementary mechanisms occurring in heterogeneous catalysis took a great step forward in 1915 when Irvin Langmuir, who worked at the General Electric research laboratory in Schenectady, New York, developed a simple but quantitative theory of adsorption of gases on metallic surfaces that predicted the surface coverage as a function of the gas phase conditions. In 1925, British chemist Hugh Stott Taylor, who worked at Princeton University, distinguished physisorption and chemisorption, and suggested that a catalytic reaction is not catalyzed over the entire solid surface of the catalyst but only at certain locations, which he called active sites. As these chemically active sites can be sparse on the catalyst surface, the reaction can be inhibited by relatively few molecules. Subsequently, in 1926, English chemist Cyril Norman Hinshelwood proposed a general mechanism for reactions on surfaces, and used the Langmuir adsorption isotherm to derive the wellknown LangmuirHinshelwood rate equations for catalytic surface reactions, in which two molecules, both adsorbed on the catalyst surface on adjacent sites, give rise to a bimolecular surface reaction.
Based on the kinetic theory of gases, German chemist Max Trautz and British scientist William Lewis, in 1916 and 1918 respectively, independently derived and published the collision theory, unaware of each other’s work due to World War I. The collision theory starts from the premise that a reaction can only happen when two molecules collide. In this representation, molecules are symbolized by two colliding spheres that undergo reaction if their kinetic energy exceeds the activation energy E, or else simply fly apart after an elastic collision if they do not have enough kinetic energy. Using the Boltzmann distribution for the velocities of the molecules, the reaction rate contains the concentrations of the reactants, the rate of collisions, and an exponential factor containing the activation energy. In this way, the collision theory provides a framework for calculating the preexponential factor of a reaction. This theory gives relatively accurate values for preexponential factors of reactions between monatomic species (as expected from its basic assumptions), but often fails by several orders of magnitude for reactions between molecules with more complicated structures.
In 1922, Frederick Alexander Lindemann proposed a multistep collision model to explain the variation of the apparent reaction order of monomolecular decomposition reactions. He assumed that these monomolecular decomposition reactions consist of three elementary steps, in which molecules were activated by bimolecular collision processes and then either followed by monomolecular dissociation steps or ‘deactivated’ by a bimolecular collision. Soon after that, in 1926, Cyril Norman Hinshelwood developed a more quantitative approach for the Lindemann mechanism by modeling the internal degrees of freedom of the reactant as several equivalent simple harmonic oscillators, and by using statistical methods to determine the probability of the molecule being activated. Further extensions were developed, in 1927 by Oscar Knefler Rice and Herman Carl Ramsperger, and independently by Louis Stevenson Kassel in 1928, by assuming that a molecule is a system of loosely coupled oscillators, leading to the RiceRamspergerKassel (RRK) theory for unimolecular gas reactions.
Around 1930, Michael Polanyi, director of the Haber Laboratory at the Kaiser Wilhelm Institute of Berlin, worked together with a young MexicanAmerican visiting scientist named Henry Eyring on the application of the wave equation for hydrogen to the quantum mechanical description of the atom exchange reaction between ortho and parahydrogen molecules. Their starting point was Fritz London’s approximate solution of the wave function for the atomic transition complex H_{3}. Polanyi and Eyring separated the nuclear and electronic motions, now referred to as the BornOppenheimer approximation, and calculated the energy of the system by considering that the reacting system’s electronic structure would rearrange to an equilibrium state in every nuclear configuration. They refined their calculations by using spectroscopic results to refine their estimates of electronic energies, creating an innovative ‘semiempirical’ method. In their fundamental paper from 1931, they described bimolecular gas reactions between three atoms and were able to exactly calculate, by means of quantum mechanics, the energy of the molecular and atom states, and hence the activation energies, relying both on quantum theory and on experimental data. In this way, they constructed the first Potential Energy Surface by plotting the energy as a function of the distances between the three atoms.
After Henry Eyring returned to the United States to become an instructor at Princeton University in 1931 and Michael Polanyi moved from Berlin to the University of Manchester in 1933, both continued their previously joint work on the energetic description of all configuration states of a chemical system and the practical exploration of PES. In 1935, with only one month difference, both Henry Eyring in Princeton, and Michael Polanyi with Meredith Gwynne Evans in Manchester, independently published a paper on a theory that allowed the description of elementary steps and the calculation their preexponential factors, known as the AbsoluteRate Theory, ActivatedComplex Theory (the term coined by Eyring), or Transition State Theory (the term coined by Polanyi). This theory combines quantum mechanics (to evaluate the PES) and statistical mechanics (to calculate the reaction rate). By evaluating the ratio of the partition functions of the transition state in the PES diagram and the fundamental state of the reactants, they obtained an Arrheniustype equation that allows the calculation of the values of both the preexponential and the exponential factors; the Eyring equation. With this, the last fundamental element of the theoretical framework was added to the classical ‘Chemical Kinetics’ approach.
Since the development of PES and the transition state theory in the early 1930s, many theoretical concepts have been introduced, further developed, refined and extended. Hence, only a short overview will be given. In 1931, Lars Onsager worked on nonequilibrium kinetics and postulated the Onsager regression hypothesis, stating that microscopic fluctuations and macroscopic relaxation are linked, some twenty years before it was finally proven to be true by the socalled fluctuationdissipation theorems. In 1931, Jens Anton Christiansen studied reactions in open sequences, proposed the catalytic cycle, and derived general formulae for reactions on single catalytic sites. Around 1936, he also developed a reaction rate theory by modeling a chemical reaction as an intramolecular diffusion process. This diffusion description was fully systematized in 1940 by Dutch physicist Hendrik Anthony Kramers and is complementary to quantum mechanical transition state theory, which was shown to be a particular case of Kramers’ pure classical method. In 1934, Soviet physicist Nikolay Nikolayevich Semenov was the first to describe the socalled branched chain reaction. He determined the various steps in the mechanisms of chain processes and developed a general and quantitative chain chemical reaction theory that considered both branched and unbranched chain processes. Semenov applied his theory to various types of chain reactions and, more especially, to combustion reactions to study flame spreading, detonation, and the burning of explosives. Several of these aspects were also independently derived by Cyril Norman Hinshelwood, but Semenov extended this work by proposing a theory of degenerate branching, which led to a better understanding of the phenomena associated with the induction periods of oxidation processes. In 1939, Japanese scientist Juro Horiuti, who worked on electrode kinetics and heterogeneous catalysis, introduced the concept of stoichiometric numbers for a reaction mechanism, leading to the derivation of generalized rate expressions. Together with Michael Polanyi, he also proposed the HoriutiPolanyi mechanism for heterogeneouslycatalyzed hydrogenation reactions in 1934. In 1943, Olaf A. Hougen and Kenneth M. Watson extended the LangmuirHinshelwood approach by explicitly considering the reverse reaction. They also used the same approach for reaction mechanisms in which the adsorption or the desorption step was rate limiting. After World War II, Rudolf Arthur Marcus extended the RRK theory for unimolecular gas reactions in 1952 by accounting for the individual vibrational modes of the system, in line with transitionstate theory, thereby enabling the calculation of simple estimates for unimolecular reaction rates from a limited number of characteristics of the PES. In 1954, Dutch scientists Pieter Mars and Dirk Willem van Krevelen described their wellknown multistep reaction mechanism for heterogeneouslycatalyzed oxidation reactions, in which lattice oxygen is used to oxidize the reactant while gasphase oxygen replenishes the lattice vacancies. In 1959, Soviet chemist Boris Pavlovich Belousov published the existence of nonlinear oscillating reactions, now known as the BelousovZhabotinsky reaction, and initiated the field of modern nonlinear chemical dynamics, illustrating that chemical reactions do not have to be dominated by equilibrium thermodynamic behavior. In 1966, Michel Boudart proposed the concept of turnover frequency. He also coined the terms Most Abundant Surface Intermediate (MASI) and structuresensitive reaction. A last trend worth mentioning, but which falls outside the direct scope of chemical kinetics, concerns direct observations and analytical measurements: from the 1950s onwards, kinetic studies were also strongly influenced by the development of various new analytical techniques, and in the case of heterogeneous catalysis, by novel surface study techniques, which culminated in the work of Gerhard Ertl, who provided detailed observations showing how chemical reactions take place on surfaces.
1 Theoretical Background
1.1 Elementary Reactions
An elementary reaction, also called an elementary step or elementary process, is a direct transformation of one or more reactants into one or more products, in which all bonds between the reacting species are broken and formed simultaneously. Figure 1 schematically illustrates the potential energy diagram for the gas phase synthesis of hydrogen iodide. When the collision between a hydrogen molecule and a iodine molecule provides enough energy to create the transition state, this elementary reaction can proceed to form hydrogen iodide.
Figure 1. Illustration of an elementary reaction. 
Most elementary reactions are mono or bimolecular, even though trimolecular elementary reactions have been reported on rare occasions. Higher order reactions cannot be considered as elementary reactions. Indeed, the probability that three or more reactant molecules are at the same time in the same location with the right configuration, conformation and orientation, and with sufficient energy for reaction is extremely low or nonexistent.
Elementary reactions have the following properties:

the reaction proceeds in a single step;

the reaction passes over a single transition state;

no reaction intermediates appear during the transformations of reactants into products;

the reaction can proceed in both directions over the same transition state (principle of microscopic reversibility);

the reaction rate for the forward and the backward elementary step is given by Guldberg and Waage’s mass action law:

the overall order of reaction is given by its molecularity;

the reaction order with respect to each reactant is equal to the stoichiometric coefficient of the reactant;

the rate coefficient is independent of composition and follows the Arrhenius law:

the ratio of the forward to the backward rate constants equals the thermodynamic equilibrium constant:
The latter equation is called the 1^{st} relation between thermodynamics and chemical kinetics. By substituting the Arrhenius law and the van ‘t Hoff relation into this equation, and developing the equations for the different thermodynamic reference states, one also derives the 2^{nd} relation between thermodynamics and chemical kinetics, which relates the heat of reaction (ΔH_{R}) to the difference between the forward and the backward activation energies (E_{+} − E_{}). Further substitution yields the 3^{rd} relation between thermodynamics and chemical kinetics, which relates the reaction entropy (ΔS_{R}) to the ratio of the forward and the backward preexponential factors (A_{+}/A_{}).
Given the abovelisted properties, the kinetic treatment of elementary reactions is relatively simple, since no reaction intermediates appear between reactants and products.
1.2 Multistep Reactions
When not all bonds are broken and formed at the same time, intermediates appear in the reaction pathway going from the reactants to the products. In such a case, the reaction is called an overall or multistep reaction, but many other terminologies are used: global reaction, apparent reaction, operational reaction, complex reaction, composite reaction, multiple step reaction, stepwise reaction, etc. An example of a multistep reaction is illustrated in the schematic energy diagram in Figure 2. As can be seen, in this overall reaction, the reactants A and B are transformed into products R and S in four steps, passing over four transition states and leading to three distinct intermediates. Put differently, this overall reaction can be decomposed into a sequence of four individual substeps or elementary steps. The information on the precise course of reaction in progress, i.e. the set of elementary steps, is called the reaction mechanism.
Figure 2. Illustration of a complex reaction. 
Most chemical reactions are global (i.e. nonelementary) reactions. In many cases, the intermediates are shortlived highly reactive species. Hence, these reactive species are often present in trace concentrations in an analyzed phase, or are species in a nonanalyzed phase (e.g., an adsorbed phase) such as the catalyst surface.
From the thermodynamic side, a number of thermodynamic relations can be written that link the individual elementary steps to the overall reaction. To this end, we will first define the stoichiometric coefficients and the stoichiometric numbers. The stoichiometric coefficients ν_{j} apply to a component in a reaction, and arise from the atomic balance for the reaction. The stoichiometric numbers apply to an elementary step in an overall reaction, and indicate how many times this elementary step is required in the reaction mechanism to arrive at the overall reaction. For an overall reaction with N species whose reaction mechanism consists of s elementary steps, the equilibrium constant of the overall reaction is given by the product of the equilibrium constant of each elementary step to the power of its stoichiometric number:$${K}_{\mathrm{eq},\mathrm{overall}}=\prod _{j}^{N}{\left({C}_{j,\mathrm{eq}}\right)}^{{\nu}_{j}}=\prod _{i}^{s}{\left({K}_{i}\right)}^{{\mathrm{\sigma}}_{i}}$$(4)
If one only looks at the enthalpy term is this equation, one can derive that the overall heat of reaction corresponds to the heat of reaction of each elementary step, multiplied by its stoichiometric number:$${\Delta H}_{R,\mathrm{overall}}=\sum _{i}^{s}{\sigma}_{i}\xb7{\Delta H}_{R,i}$$(5)
An example of these thermodynamic relations is given in Figure 3 for the catalytic synthesis of ammonia from hydrogen and nitrogen.
Figure 3. Illustration of the relations between the thermodynamic parameters of the individual elementary steps and the thermodynamic parameters of a complex reaction. 
On the side of chemical kinetics, the objective is to be able to calculate the rate of the overall reaction. In chemical engineering approaches, empirical expressions such as power laws are often fitted to a set of experimental data without accounting for the reaction mechanism and the various reaction intermediates. By proceeding in such a manner, the physical and chemical background is completely ignored, and such power law models often have a very limited applicability. The rigorous approach consists of applying the mass action law to each individual step, leading to a set of s rate equations. For each species, a net production rate can now be written as a function of the reaction rates:
Depending on the detail that is required and the information that is available, three main solution methods exist to calculate the rate of the overall reaction:

the complete solution,

the quasisteady state approximation,

the ratedetermining step approximation.
The complete solution is obtained by solving the kinetics for the most general case. In this approach, the evolution of the N species is obtained by solving the set of N production rates (containing the s mass action law rate equations of the elementary steps) without any approximations or additional hypotheses. For extremely simple cases, an analytical solution can be found, but most reaction mechanisms require the solution of these equations numerically. The advantage is that this approach is able to describe the transient kinetic behavior of the reaction, including the dynamics of all reaction intermediates. Numerical solutions are in most cases obtained by using the ‘meanfield approach’. The meanfield approximation assumes:

that all species are randomly distributed,

that the interaction of all other individuals on any given individual may be replaced by a single averaged or effective interaction, i.e. the mean field that is due to the neighboring particles.
The meanfield theory is based on the assumption that the fluctuations around the average value of the variables (pressure, temperature, concentrations, etc.) are sufficiently small to be averaged out over large ensembles and hence neglected. While the underlying approximation may easily be fulfilled in gaseous or liquid systems, this is not necessarily true for reactions on solid surfaces, where the interactions between adsorbed species may cause the formation of ordered phases. Additionally, energetic heterogeneities on the surface, surface diffusion and attractive or repulsive interactions between adsorbed species will lead to stronger deviations from the meanfield assumptions. When the meanfield approximation breaks down, the numerical solution has to be calculated by means of MonteCarlo (MC) simulations, which offer in principle an exact description of the interplay between the elementary steps. A major drawback in all of the above solution techniques (analytical solution, meanfield approximation, or MC simulations), however, is that all kinetic parameters of all elementary steps need to be known in order to be able to calculate the complete solution.
The second approach to finding the reaction rates is the QuasiSteady State (QSS) approximation. In order to simplify the set of equations and reduce the number of rate parameters in the system, this approach assumes that there is no net rate of production for the M intermediates in the reaction mechanism. Indeed, highly reactive intermediates will react as soon as they are formed. Hence, after an initial induction period, these intermediates reach a QSS concentration, generally in inferior concentrations. Although an actual steadystate can only be truly reached in open reactors, the QSS approximation is also applied for closed systems (batch reactors) as relaxation times (the time taken to reach the steadystate) are generally shorter than the reciprocal turnover frequency. For each intermediate, its net production rate is therefore set to zero, resulting in a subset of M algebraic equations in M unknown concentrations. Solving this subset of equations makes it possible to obtain the concentrations of the M intermediates as a function of the concentrations of the remaining measured species and the rate parameters. In many cases, analytical expressions can be obtained for the concentrations of the M intermediates, leading to a closed analytical expression for the rate equation of the overall reaction. The main advantages of the QSS approximation are that fewer kinetic parameters are required to obtain the reaction rate, and that an analytical expression is often obtained. In contrast to the complete solution, timedependent phenomena within the full reaction mechanism cannot be accounted for, such as branching chain reactions that lead to explosions.
The third approach is the Rate Determining Step (RDS) approximation, also called the quasiequilibrium approximation. In this approach, one assumes that one of the elementary steps determines the overall reaction rate, while all other elementary steps are sufficiently fast to be considered as being in quasiequilibrium. Hence, for a reaction mechanism with s elementary steps, the assumption of quasiequilibrium for all but one elementary step leads to a subset of s − 1 algebraic equations. This allows the expression of the concentrations of the various intermediates as a function of the concentrations of the measured species and the equilibrium constants of the elementary steps. Equaling the overall reaction rate to that of the ratedetermining step and eliminating the concentrations of the intermediates systematically leads to an analytical expression for the rate equation that depends on the concentrations of the measured species, the rate constant of the ratedetermining step, and the equilibrium constants of the steps in quasiequilibrium. In addition to the latter advantage, the ratedetermining step approach leads to even fewer kinetic parameters which need to be determined. As the assumptions are stronger, the resulting rate equation is a more limited description of the reaction mechanism. In addition to the fact that the timedependent dynamics of the full mechanism cannot be described, as for the QSS approximation, the rate equation from this approach is no longer able to account for a shift in the ratedetermining step as a function of the reaction conditions.
Besides these three main solution methods, other variants in solution techniques or hypotheses can also be used. One can use two (or more) ratedetermining steps to allow the rate equation to shift between ratedetermining steps as a function of reaction conditions, an approach which is intermediate between the QSS and the RDS hypotheses. If one wants to further reduce the complexity of the rate equation and hence the number of rate parameters, one can use the irreversible step approximation by imposing that one or more steps in the reaction mechanism are irreversible, resulting in an even less general rate expression. For catalytic reactions, one can choose to neglect a number of adsorbed reaction intermediates by using either the Most Abundant Reaction Intermediate (MARI) approximation, sometimes called the MASI approximation, which considers that the catalytic surface only consists of empty sites and MARI sites, or the Nearly Empty Surface (NES) approximation, by considering that all intermediates are very weakly bound and that the catalytic surface is almost empty. In all these cases, a choice must be made between the generality of the rate equation and the number of rate parameters that need to be provided.
1.3 Reaction Rate Constants
As already recognized by Arrhenius, rate processes are characterized by rare events, or, formulated differently, statechanging events (reactions) take place on a long time scale when compared to the dynamic time scales involved in the (vibrational, rotational, torsional, translational) motion of molecules. Reaction rate theories aim to understand how an elementary step actually happens and to predict the rate parameters of a reaction, i.e. the preexponential factor and the activation energy, from these microscopic descriptions. Although several rate theories have been developed to calculate these kinetic parameters theoretically, two main approaches have been developed and improved over the years: the collision theory and the transition state theory.
In the collision theory, an elementary reaction happens when suitable reactants collide with sufficient energy and with the correct orientation and conformation. Based on the kinetic theory of gases and the MaxwellBoltzmann distribution for the velocities of the molecules, the following reaction rate is obtained for a bimolecular reaction:$$r={\left({\sigma}_{\mathrm{AB}}\right)}^{2}\xb7{\left[\frac{8\pi \mathrm{RT}}{{M}_{\mathrm{AB}}}\right]}^{1/2}\xb7{N}_{A}\xb7{e}^{\frac{{E}_{a}}{\mathrm{RT}}}\xb7{C}_{A}\xb7{C}_{B}$$(6)
From this expression, the preexponential factor can be isolated and directly calculated on theoretical grounds:$$A={\left({\sigma}_{\mathrm{AB}}\right)}^{2}\xb7{\left[\frac{8\pi \mathrm{RT}}{{M}_{\mathrm{AB}}}\right]}^{1/2}\xb7{N}_{A}$$(7)
The collision theory gives relatively accurate values for preexponential factors of reactions between monatomic species (as expected from its basic assumptions), but often fails by several orders of magnitude for reactions between molecules with more complicated structures. The reason for this is that the reacting molecules have been considered to be spherical and hence able to react in all directions. For actual molecules, this is no longer true, as the orientation of the collisions and the conformation of the molecules is not always suitable for reaction.
The most influential rate theory is the transitionstate theory. It is based on the evaluation of the barriers on the PES between the transition state, or activated complex, and the reactants. This theory combines quantum mechanics (to evaluate the PES) and statistical mechanics (to calculate the reaction rate). By evaluating the ratio of the partition functions of the transition state in the PES diagram, and the fundamental state of the reactants, the Eyring equation is obtained for the rate constant:$$k=\frac{{k}_{B}\xb7T}{h}\xb7{e}^{\frac{\Delta {S}_{c}^{\ne}}{R}}\xb7{e}^{\frac{\Delta {H}_{c}^{\ne}}{\mathrm{RT}}}\xb7{\left({c}_{0}\right)}^{{\Delta n}^{\ne}}$$(8)
The transitionstate theory tries to predict the rate constant from the initial reactive flux of molecules going over the activation barrier, and will therefore always overestimate the rate constant. Although Eyring recognized the overestimation and introduced a ‘transmission coefficient’, it has not been possible to calculate the transmission coefficient within the transitionstate theory itself. Hence, the transitionstate theory provides a clear upper bound for the rate constant. Even with this limitation, additional difficulties appear, first of all in obtaining a precise knowledge of PES, and secondly in finding the correct configuration of the activated complex. Indeed, the most accurate quantum chemical techniques are computationally very demanding and can therefore only be applied for fixedgeometry energy calculations. Cheaper methods, such as DFT, are needed to map out the PES and to calculate the vibrational frequencies of the transition state. Unfortunately, the DFT functionals are currently not sufficiently reliable to provide accurate rate constants. Hence, the transition state theory, an educated guess for the nature of the activated complex, and modern computational chemistry tools often allow the prediction of rate constants for bimolecular reactions only within an order of magnitude, or at best within a factor of three. Alternatively, Benson (1976) gives a list of rules for estimating values for ΔS^{≠} and ΔH^{≠} that are needed in the Eyring equation on the basis of the characteristics of different types of reactions, for reactions both in the gas phase and in solution. Benson also gives tables of ‘group contributions’. For fast estimates of organic reaction rates, Benson’s rules are often the method of choice. The transitionstate theory is therefore primarily used to understand how a chemical reaction takes place, but has been less successful in its original goal of directly calculating rate constants on theoretical grounds.
Since rate parameters can still not be determined with sufficient accuracy from reaction rate theories, the preexponential factors, activation energies, adsorption constants, etc., are in most applications determined based on experimental data. This is generally done by parameter estimation, also called parameter identification, techniques that determine the rate parameters in the rate equations by minimizing the deviations between the model predictions and the experimental results.
The classic approach starts by carrying out a large set of kinetic experiments. From the analysis of these data, one or several reaction mechanisms are then proposed, and their corresponding rate equations are derived. After estimating the parameters of each of the competing kinetic models, the various models are ranked according to a goodnessoffit criterion, and the best kinetic model is finally selected. In this approach, since the model selection is performed on the current data set (operating conditions, feedstock compositions, etc.), the selected reaction mechanism is often oversimplified and one should carefully investigate the reliability and robustness of the resulting kinetic model, especially when extrapolating outside the domain of the data set.
More recently, the parameters of the rate equations can be determined on the basis of a microkinetic analysis. Such a microkinetic analysis uses input from surface science, from computational chemistry, and from in situ or in operando kinetic studies to gain insight into the reaction intermediates and the actual elementary steps, and to independently determine the value of some of the rate parameters from separate experiments, thereby reducing the number of rate parameters that need to be estimated from the kinetic experiments. In this way, a more fundamental kinetic model can be obtained that integrates the knowledge of the various research fields.
2 Description of the Problem
In the sections above, the general framework for chemical kinetic modeling has been described. This theoretical background allows the development of kinetic models not only for simple systems, but also for very complex processes. In the latter case, several difficulties arise very rapidly.
Modeling complex processes usually involves several thousand components and reaction intermediates. In the case of petroleum refining processes, the number of molecules and the number of reactions grows almost exponentially with carbon number. Even the enumeration of the isomers quickly becomes impossible without computer algorithms (Tab. 1). To give some orders of magnitude, processes treating naphtha fractions typically deal with up to a thousand components that contain between 5 and 11 carbon atoms, while for processes treating middle distillates, there are already several hundred thousand components that contain typically between 12 and 24 carbon atoms. For vacuum gas oil fractions that typically contain molecules with 25 to 40 carbon atoms, the number of possible isomers exceeds 10^{14}, even when heteroatoms such as sulfur, nitrogen and oxygen are not accounted for. For Vacuum Residues (VR) containing molecules with more than 40 carbon atoms and with high amounts of heteroatoms, the number of isomers is almost intractable.
Number of possible isomers of Alkanes, Alkenes, Alkynes (Cayley, 1875; Read, 1976; Hendrickson and Parks, 1991; Bytautas and Klein 1998; Faulon et al., 2005)
Between the components of this extremely large set of molecules, an even larger number of reactions can occur. The reaction network will need to describe the chemistry as precisely as needed for the understanding of the various phenomena, but also to improve the prediction of the process performances and the product qualities. Again, full reaction networks can only be computergenerated.
Associated with these large reaction networks are network reduction techniques. Indeed, for complex feedstocks, the reaction network becomes so large that it becomes very difficult to handle during the simulation. Two main classes of network reduction techniques can be distinguished: a priori lumping techniques and a posteriori lumping techniques. In the former, the full network is not created, but a smaller version of it is proposed, grouping components into ‘lumps’ and reactions into ‘lumped reactions’. In a posteriori lumping techniques, the full reaction network is first generated, before it is reduced into a more tractable network by a scientific approach.
For each reaction of the final reaction network, a rate equation needs to be proposed or generated, to which a number of rate parameters will be associated. For large reaction networks, one not only needs to propose ways to automatically generate the rate equations, but also to reduce the number of rate parameters that need to be identified.
Last but not least, in order to be able to simulate the evolution of the composition in the reactor, one needs to provide an appropriate description of the feed, which can either be a molecular level description of the feed when detailed modeling techniques are used, or a lumped description of the feed.
No matter which modeling approach is selected to simulate complex chemical systems, the four following steps need to be clearly defined first:

Feedstock composition: What is the required detail?

Reaction network: What reactions need to be considered?

Rate equations and rate parameters: At what rate is each component formed?

Continuity equations: How is the evolution governed inside the reactor?
In this review, detailed reactor models (including vaporliquid equilibrium thermodynamics, heat and mass transfer between phases or inside catalyst particles, external heat exchangers, etc.) and their corresponding continuity equations will not be discussed.
In what follows, systems with several hundred or several thousand components, reaction intermediates, reactions and elementary steps will be considered. In such cases, two main strategies can be applied: a lumping strategy or a detailed modeling strategy. In the first strategy, the chemical complexity is reduced by grouping chemical compounds into families or ‘lumps’ and grouping individual reactions into lumped reactions between these families. In this case, a meanfield approach together with a ratedetermining step approximation is generally used, as will be shown in the examples below. The corresponding rate equations are then written at a macroscopic level for a limited number of pathways between the analytically observable lumps of components, while the corresponding continuity equations are solved by means of a classic determinist solver for a set of ordinary or Partial Differential Equations (PDE). In the detailed modeling strategy, the chemical detail is generally retained at the molecular level, or even at the level of the reaction intermediates. For these cases, we will describe techniques to provide a detailed feedstock composition, to generate and reduce full reaction networks, to derive the corresponding rate equations, and to decrease the number of rate parameters. Concerning the reaction rates, one can either use a meanfield approach and apply the QSS approximation, or calculate the complete solution by means of MC methods, as will be illustrated below.
3 Approach by Lumping
3.1 General Information
The lumping approach consists of regrouping chemical compounds by similar properties called ‘lumps’. The lumps are then considered as homogeneous ensembles on which a kinetic model normally used for the molecular compounds can be applied. This approach is often used for processes where the molecular characterization of the reactant mixtures is difficult or impossible because the feedstock is too complex, as is the case in the majority of petroleum refining processes (catalytic reforming, hydrotreating, hydroprocessing, catalytic cracking, thermal coking, etc.).
The development of a lumping approach usually proceeds through the following steps:

description of the feedstock by choosing a set of lumps;

description of the relationships between the lumps by building a kinetic network of lumped reactions;

determination of the rate equations and their associated parameters by optimizing the model on experimental data.
The choice of lumps is always a compromise between the capabilities of the analytical techniques to characterize and quantify them on one hand, and the needs of the final user in terms of model prediction and precision on the other. In most cases, the analytical techniques are the limiting step and force the choice of the lumps. In the field of petroleum refining, these analyses can be subdivided into chemical separation techniques (Mass Spectrometry (MS), Liquid Chromatography (LC), Solvent Extraction (SE), etc.) and into physical separation techniques (distillation mainly). Historically, lumped models were developed from the 1960’s onwards, at a time when petroleum analyses were very simple but yielded little information, and when computing power was limited. By way of example, the first lumped model for the Fluid Catalytic Cracking (FCC) process proposed by Weekman and Nace (1970) contained only 3 lumps corresponding to 3 distillation cuts: ‘Feed’, ‘Gasoline’ and ‘Gas+Coke’. Over time, thanks to the development of more efficient separation techniques and the increase in computing power, lumped models became more and more complex with a continuous increase in the number of lumps. For the catalytic cracking process, Jacob et al. (1976) proposed an extended 10lump model in 1976, while in 1994, Pitault et al. (1994) further increased the complexity by introducing an 18lump kinetic network. More recently, at the end of the 1990s, the lumping approaches were deeply modified with the introduction of ‘reconstruction methods’, which allow the definition of several hundred to several thousand chemical lumps from a set of limited analyses and expert knowledge. For example, Christensen et al. (1999) proposed a kinetic model for catalytic cracking with approximately 3000 lumps that is based on the Structure Oriented Lumping (SOL) approach developed by Quann and Jaffe (1992). These novel techniques will be detailed in the second part of this article. Figure 4 illustrates the evolution of the lumped kinetic models for catalytic cracking, as proposed by Weekman and Nace (1970), Jacob et al. (1976) and Pitault et al. (1994).
Figure 4. Illustration of the evolution of the lumped kinetic models for the catalytic cracking process: (a) Weekman and Nace (1970); (b) Jacob et al. (1976); (c) Pitault et al. (1994). 
3.2 Advantages and Drawbacks of Lumped Kinetic Models
Despite their increasing complexity, lumped kinetic models are relatively easy to develop because the number of lumps and the number of reactions remain limited. Moreover, due to the multicompound characteristics of the lumps, the reaction pathways are generally global with no intermediate species and the kinetic rate equations are often simple (pseudoorder reactions, Langmuir approach in heterogeneous kinetics, etc.). Their kinetic parameters (preexponential factors, activation energies, adsorption constants, thermochemical constants, etc.) are often determined by minimizing the deviations between model and experimental data coming from pilot units or industrial plants. This simplicity allows the collection of highspeed kinetic models that require limited computing power, a feature that is very interesting for the optimization and control of petroleum processes. This explains why this type of model has been predominant in petroleum refining for over 50 years.
However, despite this advantage, lumping methods also have some drawbacks that need to be kept in mind. Firstly, the lumps are often defined as ensembles of compounds with similar physicochemical properties (determined by the analytical techniques) but not necessarily with similar reactivities. If the thermodynamic equilibria inside the lumps are not maintained by fast intralump reactions, the properties of the lumps will consequently be modified during the reaction due to the transformation of these internal compounds as a function of their own reactivities. In this case, the base hypothesis that considers that the properties of the lumps are homogeneous and constant during reaction is clearly wrong (Li and Rabitz, 1989; Wei and Kuo, 1969).
The second drawback of lumped models is that they are not associated to a molecular kinetic theory but are directly derived from experimental data coming from pilot units in most cases. In petroleum refining, these experiments are long and expensive. Consequently, even if the theoretical formulation of lumped kinetic models is relatively quick and cheap, their parameter estimation requires more time and money. Indeed, in order for lumped models to be robust and feed independent, a wide variety of experimental data in terms of operating conditions and even more importantly, feedstock composition is needed. Moreover, with the increase in both the complexity of the kinetic models and in the number of lumps, the analyses needed to describe the mixture become more and more important, further increasing the experimental cost for the development of lumped models.
Finally, the last main drawback of lumped kinetic models is related to their relative inability to determine the physicochemical properties of the effluent from its composition and the properties of its lumps. In practice, these models often use correlations in order to obtain an estimation of the desired product properties, but this approach is empirical and very limited. Moreover, because the properties of the lumps can change during the reaction, this estimation method remains somewhat arbitrary.
3.3 Examples of Lumped Kinetic Models Developed at IFPEN
Between 1995 and 2005, IFPEN developed lumped kinetic models for several hydroprocessing processes such as Atmospheric Gas Oil (AGO) hydrotreating or fixedbed residue hydroprocessing.
One of the first lumped models developed at IFPEN for the AGO hydrotreating was a 9lump model (Fig. 5a) accounting for only 1 sulfur lump for dibenzothiophenes (+ H_{2}S), 1 nitrogen lump for carbazoles (+ NH_{3}) and 4 hydrocarbon lumps (+ H_{2}) representing saturates, monoaromatics, diaromatics and triaromatics (Bonnardot, 1998; MagnéDrisch, 1995). The analyses required to define these lumps were very simple: elementary sulfur content for the sulfur lump, elementary nitrogen content for the nitrogen lump and MS for the other lumps. The reaction network contained 6 overall reactions. The reactions for aromatic hydrogenation were considered to be reversible, while the reactions for hydrodesulfurization and hydrodenitrogenation were defined as being irreversible. The rate equations were derived from the following hypotheses: existence of two types of active sites for hydrogenation and for hydrogenolysis, Quasi Steady State Approximation for the intermediate species, LangmuirHinshelwood approach, and equilibrium for adsorption/desorption reactions. The experimental domain was relatively large in terms of H_{2} partial pressure (from P_{base} –22 bar to P_{base} +28 bar), H_{2}S partial pressure (from 0.5 to 5 bar) and Liquid Hourly Space Velocity (LHSV) (from 1 to 3 Sm^{3}/m^{3}/h), but only one temperature (T_{base}) and only one type of feedstock (two full range Light Cycle Oils, denoted LCO, and three LCO aromatic extracts) was used. The kinetic model was implemented in a single phase plugflow reactor model, and contained 15 parameters, which were determined from a database of 90 experimental data points with 5 independent responses each. In 2000, López García improved the prediction of the model in terms of hydrodesulfurization (López García, 2000; López García et al., 2003) by introducing 3 sulfur lumps which represent the 3 most refractory sulfur classes of the AGO: DBT, 4DBT and 46DBT (Fig. 5b). The determination of the 3 sulfur lumps was determined by Gas Chromatography (GC) coupled with a sulfurspecific SCD detector. This extended reaction network contained 9 overall reactions. To reduce the number of rate constants in the rate equations, the author applied the RDS approach for the various reactions on the two types of active sites: hydrogenation sites and hydrogenolysis sites. The experimental domain was extended in terms of operating temperature (from T_{base} −20 °C to T_{base} +10 °C) and Liquid Hourly Space Velocity (from 0.5 to 4 Sm^{3}/m^{3}/h). Although three full range LCO and five LCO aromatic extracts were used, the feedstock variability was still limited, as only one type of feed (FCC derived feeds) was used. Further improvements concerned the introduction of thermodynamic constraints for the reversible reaction, and the use of a twophase plugflow reactor model based on a GraysonStreed flash calculation. By introducing the effect of the temperature, the number of parameters was increased to 31, which were determined from a database of 122 experimental data points with 8 independent responses each. To further improve the prediction capabilities of the kinetic model for the AGO hydrotreating process, the direct lumping approach based on the analytical capabilities was abandoned, and a new kinetic model based on a feedstock reconstruction approach (see Sect. 4.1.2) was developed. For this model, the number of lumps were classified by chemical family and by carbon number (López García et al., 2010), resulting in a total to 597 lumps. The kinetic network was also extended by accounting for 14 reaction families related to hydrogenation/dehydrogenation, hydrodesulfurization and hydrodenitrogenation pathways (Fig. 5c). The corresponding rate equations were derived based on the presence of two types of active sites (hydrogenation and hydrogenolysis) using the RDS approach. The experimental domain encompassed a large domain of operating temperatures (from T_{base} −20 °C to T_{base} +20 °C), of operating pressures (from P_{base} –57 bar to P_{base} +30 bar), and of Liquid Hourly Space Velocities (from 0.3 to 4 Sm^{3}/m^{3}/h). The largest extension concerned the feedstocks, as 24 different types of industrial gas oils were included: straight run gas oils from a large variety of crude oils, LCO, coker gas oils and their mixtures. This feed diversity is particularly important to confer a high degree of robustness on the model.
Figure 5. Evolution of the IFPEN lumped kinetic models for atmospheric gasoil hydrotreating: (a) Bonnardot (1998); (b) López García (2000); (c) López García et al. (2010). 
Concerning fixedbed residue hydroprocessing, IFPEN developed a series of lumped kinetic models based on two main separation analyses, the first by volatility (distillation), the second by polarity (Saturates/Aromatics/Resins/Asphaltenes (SARA) LC) (Haulle, 2002; Le Lannic, 2006). Combining both analyses allowed the definition of 5 petroleum fractions: asphaltenes (ASP), resins (RES), aromatics of the 520 °C+ cut (ARO+), aromatics of the 520 °C− cut (ARO−), and saturates (SAT). For each of these fractions, an elemental analysis was performed in order to determine its composition in terms of carbon, hydrogen, sulfur, nitrogen, oxygen, nickel and vanadium. The first kinetic scheme described the evolution of the quantity of these different fractions, with an additional specific focus for the sulfur, vanadium and nickel elements (Fig. 6a). The total number of lumps for this approach was equal to 19. The reactor was modeled as a singlephase plugflow reactor with a LangmuirHinshelwood formalism to manage the adsorption of the different species on the catalyst. Subsequently, this model was generalized (Le Lannic, 2006; Verstraete et al., 2007) by defining 2 additional fractions (RES+ and RES−) and by decoupling each fraction by chemical element (C, H, O, N, S, Ni, V). This new approach was more rigorous, but the number of lumps was increased to 40 (Fig. 6b). This new kinetic model was coupled with a catalyst morphology model to account for the diffusion of asphaltenes and resins inside the catalyst and for the modification of the catalyst properties (porosity, surface area, etc.) caused by the deposition of nickel sulfides and vanadium sulfides by hydrodemetallization (Toulhoat et al., 2005). A third kinetic and reaction model modified both the lumped kinetic network (Fig. 6c), by introducing additional asphaltenes families, and the reactor model, by introducing diffusion/mass transfer limitations based on the StephanMaxwell equations (Ferreira, 2009; Ferreira et al., 2014, 2010).
Figure 6. Evolution of the IFPEN lumped kinetic models for fixedbed residue hydroprocessing: (a) Haulle (2002); (b) Le Lannic (2006); (c) Ferreira (2009). 
4 Detailed Kinetic Modeling Strategies
Due to the complexity of the feedstock, but also because of the limitations in computer hardware and software, classical kinetic models of industrial processes, such as petroleum refining or petrochemical processes, have traditionally been based on a lumping approach, as illustrated above. However, over the last decades, these industries have been subjected to more stringent environmental legislation and tighter product quality constraints, which are phrased in terms of molecular or atomic composition of feedstocks and process products (Neurock, 1992). Novel kinetic models must therefore be able to predict the performance of processes at the molecular level. This cannot be ensured by a lumping approach due to its limitations in describing the molecular composition throughout the entire reactor.
The limitations of lumped kinetic models therefore motivated the development of more fundamental and more detailed kinetic models. Two other approaches can be distinguished to model the complex process kinetics: a mechanistic approach and a molecular approach. The difference between these approaches resides in the level of detail at which the reaction pathways are described.
Mechanistic models are the more fundamental approach based on a detailed explicit description of the reaction pathways, including its reactants, elementary steps (such as homolytic bond scission, radical recombination, olefin protonation, ion deprotonation, βscission, isomerization, hydride abstraction, etc.), and its reaction intermediates (such as ions and radicals). In this approach, a limited number of a priori assumptions are needed and the rate parameters are more fundamental in nature.
Molecular models result from an intermediate approach between the mechanistic and lumping approaches. Here, a chemical reaction system is modeled at the molecular level without including its reaction intermediates. The reactions are viewed as moleculetomolecule transitions and each reaction is characterized by an overall rate constant. The effects of reaction intermediates are included in the rate equations by imposing assumptions during their derivation.
Both detailed kinetic approaches enable one to overcome the drawbacks related to the lumping approach, since they both retain a molecular description of the reaction system throughout the entire reactor simulation, which allows the creation of feedstockindependent models. However, such models expect a molecular description of the feedstocks, of the reaction pathways, and of the rate equations and rate parameters. In what follows, methods to determine a molecularlevel description of the feedstock are first discussed in detail and illustrated. Subsequently, the simulation of detailed reaction networks is discussed by simultaneously treating the generation and reduction of the detailed reaction networks, and the rate equations, rate parameters and reactor simulation, before illustrating examples of various applications.
4.1 Molecular Description of the Feedstocks
The first step in the simulation of a detailed moleculebased kinetic model is to determine the molecularlevel description of the feedstock. This step is a crucial component of the moleculebased kinetic modeling since the feedstock description provides the major part of the ‘initial conditions’ of a molecular reaction model. The molecular description can be obtained either by using, when possible, advanced analytical characterization techniques, or by numerically creating a molecular representation of the feedstock through composition modeling algorithms, here called molecular reconstruction methods.
In this section, the existing methods of obtaining the molecular description of the feedstock will be reviewed.
4.1.1 Analytical Characterization
Common analytical techniques used to obtain a molecular characterization of the feedstocks are MS and GC. MS consists of transforming molecules, in the gaseous state, into ion fragmentation patterns by electron bombardment. These fragments are then separated and detected by their masstocharge ratio (m/e). The graphical representation of the ion intensity as a function of m/e makes up the mass spectrogram which is unique for each molecule (Wauquier, 1994). The MS technique is suitable for quantifying an individual targeted compound or a small number of compounds (2030 molecules) in a mixture. However, the main application of the MS technique, in particular in petroleum laboratories, is a quantitative analysis of the chemical families of the form C_{ n }H_{2n+Z }X where C and H are carbon and hydrogen, respectively, n is the number of carbon atoms, Z is the hydrogen deficiency and X refers to heteroatoms such as S, N, O (ASTM D2425, 2007; ASTM D2786, 2007; Fafet et al., 1999a,b; Fischer and Fischer, 1974).
GC separates volatile components of mixtures through a chromatographic column with a gaseous mobile phase by physical characteristics, such as diffusivity, adsorption, absorption and volatility. As the GC technique is a separation technique but not an identification method, the chromatography column must be coupled to a specific detector. The most common detectors are:

Flame Ionization Detector (FID) (Adam et al., 2008; ASTM D6730, 2007; Johansen et al., 1983; Vendeuvre et al., 2005),

Sulfur or Nitrogen Chemiluminescence Detector (SCD/NCD) (Adam et al., 2009; Dzidic et al., 1988; López García et al., 2002, 2003; Revellin et al., 2005; Tuan et al., 1995),

Atomic Emission Detector (AED) (Andersson and Sielex, 1996; Depauw and Froment, 1997),

mass spectrometer (Bouyssiere et al., 2004; López García et al., 2002; Teng and Williams, 1994).
GC techniques have a high separation efficient, equivalent to 2.5 10^{5} theoretical plates, which allows them to identify up to 200 compounds in a mixture (Mondello et al., 2002). Due to its high resolution, GC is widely used for quantitative analysis of petroleum gases (Speight, 1991) and gasoline fractions (C_{5}C_{10}) (Johansen et al., 1983; Teng and Williams, 1994).
The MS and GC techniques can identify and quantify several hundred molecules in a mixture. However, the number of hydrocarbon isomers increases exponentially with the number of carbon atoms. For example, the number of possible isomers for a molecule with 5 carbon atoms (C_{5}) is 8, with 10 carbon atoms (C_{10}) it is 474, and with 15 carbon atoms (C_{15}) it is over 40 000 (Hudebine, 2003). This is why, beyond C_{10}, the analytical techniques are not able to obtain a direct identification of all molecular species.
Other advanced analytical methods, such as comprehensive two dimensional Gas Chromatography (GC × GC) (Bertoncini et al., 2013; Dutriez et al., 2009, 2010, 2011; Giddings, 1987; Mahe et al., 2011; Mondello et al., 2002), HighPerformance Liquid Chromatography (HPLC) (Miller et al., 1983), and Fourier Transform – Ion Cyclotron Resonance – Mass Spectrometry (FT–ICR–MS) (Qian et al., 2001), have been developed to provide a molecular description of complex feedstocks. Today, these advanced techniques do not provide quantitative molecular information about the distribution of compounds within a complex mixture, especially for the high carbon number range, but rather structural characteristics of the molecules. Moreover, they are generally very timeconsuming, and the obtained data are often difficult to interpret.
4.1.2 Molecular Reconstruction Methods
In an effort to overcome the limitations of the analytical characterization, current research focuses on the development of composition modeling techniques, also called molecular reconstruction methods. The composition modeling consists of transforming available analytical information, such as the overall properties and molecular structural information, into a molecular representation of the feedstock.
Several molecular reconstruction methods have been developed (Al Halwachi et al., 2012; Allen and Liguras, 1991; AlvarezMajmutov et al., 2014, 2015, 2016; Boek et al., 2009; Neurock et al., 1994; Peng, 1999; Pyl et al., 2011; Quann and Jaffe, 1992; Sheremata et al., 2004; Zhang, 1999). These methods can be classified into three groups according to the approach used to represent the molecular composition: approach by model molecule, deterministic molecular reconstruction approach, and stochastic molecular reconstruction approach.
Approach by Model Molecule
The first reconstruction methods consisted of generating a model molecule to represent the petroleum cuts. The objective of these methods is to generate a single molecule that contains the average characteristics of mixture to be represented.
The first reconstruction method proposed in the literature concerned coal constituents, and was developed by van Krevelen (1952). It consisted of creating a 2D molecule for a given fraction from its carbon content, hydrogen content, refractive index, density, and aromatic carbon content. Williams developed a molecule reconstruction algorithm for the aromatic fraction of an oil sample using molecular weight, elemental analysis, and ^{1}H NMR (Williams, 1958). The algorithm was also applied to four asphalt fractions samples, to a virgin gas oil sample and to a catalytic cycle stock sample. Brown and Ladner (1960) developed a similar approach for coallike material. One of the most wellknown methods was devised by Speight (1970) and creates a 2D molecule for an asphaltene fraction from its elemental analysis, average molecular weight and proton nuclear magnetic resonance (^{1}H NMR). From this analytical information and some assumptions, Speight developed 7 mass balances for different types of carbon atoms in order to determine the average hydrocarbon structure of the petroleum cut. Almost simultaneously, Hirsch and Altgelt (1970) proposed a similar technique that also used the average density as an analytical input. The authors also suggest to replace some assumptions with 5 parameters that had to be estimated by the modeler.
Over time, these methods have been improved by adding additional analytical techniques, such as ^{13}C NMR (Ali et al., 1990, 2006; Altgelt and Boduszynski, 1994; Cantor, 1978; Dickinson, 1980; Knight, 1967; Petrakis and Allen, 1987; Sato, 1997; Sato et al., 1998; Suzuki et al., 1982; Takegami et al., 1980), Size Exclusion Chromatography (SEC) (AlZaid et al., 1998; Gauthier et al., 2008), infrared spectroscopy (Montgomery and Boyd, 1959; Qian et al., 1984), functional groups analysis (Faulon, 1994) or pyrolysis analysis (Artok et al., 1999; Faulon et al., 1990; Kowalewski et al., 1996), which allows for more molecular input for the mixture and reduces the number of assumptions.
From the analytical information, the modeler can easily build a representative model molecule of a mixture. Additionally, most of these methods only need the overall properties which can be obtained without too much experimental effort. However, a model molecule is far from being representative of all molecules in a mixture. Indeed, a model molecule does not provide the effects of polydispersity of complex mixture in terms of physical properties and reactivities, which are essential for the process modeling.
Deterministic Molecular Reconstruction Approach
The deterministic molecular reconstruction approach consists of modeling the mixture composition by means of a set of molecules whose mixture properties are close to those of the mixture to be represented. The synthetic mixture of molecules is generated by estimating the mole fractions of the predefined database of compounds. The development of such an approach is usually done in the following manner:

creation of a predefined database of compounds,

definition of the objective function (analytical constraints),

modification of the mole fractions of the compounds.
The creation of a database of compounds is extremely important for the accuracy of the molecular reconstruction method, because this step chooses which compounds will be considered. Evidently, the absence of the key compounds in the database leads to an inaccurate molecular representation of the mixture. Depending on the type of mixture that is considered, a different molecular library has to be selected. The database is therefore created from qualitative analytical information and expert knowledge. As mentioned in Section 4.1.1, advanced analytical techniques can provide molecularlevel information. This information can be a qualitative identification of the molecules for light mixtures, such as gasoline, or the identification of the molecular structures for more complex feedstocks, such as middle distillates or vacuum gas oils. Thus, the predefined database can be defined as a set of molecules, as proposed by Allen and Liguras (1991) or as a library of molecular structures, such as the SOL (Structure Oriented Lumping) approach proposed by Quann and Jaffe (1992) or the MTHS (Molecular Type and Homologous Series) method proposed by Peng (1999) and Zhang (1999).
Whatever the predefined database, the pure component properties of its compounds are required. For light compounds, pure component properties are often available in thermophysical properties databases, such as the NIST Chemistry WebBook (Linstrom and Mallard, 2015), TRC (Frenkel et al., 2000) or the DIPPR (Daubert and Danner, 1989). If the pure component properties are not available, they can be calculated either by inspection of the structure of the molecule, to assess chemical formula and molecular weight, by correlations, or by group contribution methods, which are frequently used to estimate normal boiling points, densities, and critical properties (Benson, 1976; Benson and Cohen, 1993; Constantinou and Gani, 1994; de Oliveira et al., 2013a; Hudebine and Verstraete, 2011; Joback and Reid, 1987; MarreroMorejón and PardilloFontdevila, 1999).
To calculate the various properties of a mixture of molecules from its composition, mixing rules need to be used that are based on the mole fractions of the compounds and their pure component properties. It is common practice to assume that the mixing rules are linear with respect to either the pure component property or its blending index (Hudebine and Verstraete, 2011), even if some of the properties cannot be blended linearly. For every piece of analytical information, an analytical constraint can be written. This set of analytical constraints is also associated to a number of physical constraints: a material balance to ensure that the sum of mole fractions is equal to 1, and mole fraction constraints that require each mole fraction to be positive and to fall within the range from 0 to 1.
Once the database of compounds is defined and the analytical and physical constraints have been set up, the mole fraction of each compound needs to be determined. This is generally done by minimizing the difference between the available quantitative analytical data and the properties calculated from the composition of the synthetic mixture. Even using linear mixing rules, the estimation of the mole fractions is often an illdefined problem because the number of unknowns (mole fractions) largely exceeds the number of observables (analytical data). To reduce the number of unknowns, several techniques have been proposed in the literature. Allen and Liguras (1991) apply a constrained Simplex method to minimize a weighted criterion using the enthalpies of formation of the molecules as weights, yet respecting analytical constraints. Jaffe et al. (2005) suggest the minimization of a specific criterion using a Lagrange multiplier method. Zhang and coworkers (Ahmad et al., 2011; Aye and Zhang, 2005; Peng, 1999; Wu and Zhang, 2010; Zhang, 1999) developed an approach to estimate the mole fractions of a mixture from their overall properties by interpolating into a database of wellcharacterized mixtures. Recently, Pyl and coworkers (Pyl et al., 2011) use the parameterized probability functions to describe the distribution of mole fractions in the compounds.
The deterministic molecular reconstruction approach is the most efficient of the reconstruction methods. Its conception is based on a large set of analytical data and few assumptions are needed concerning the presence or absence of a compound in a mixture. Moreover, this method is usually quite fast in determining a molecular composition of a mixture. However, this approach suffers from two shortcomings. Firstly, the development of the predefined databases is a very timeconsuming and expensive task because of the advanced analytical tools that are needed to obtain the molecular information. Of course, once these databases are created, the determination of the molecular composition of new mixtures is fast, but the initial analytical work remains an extremely limiting factor. Secondly, this approach cannot be applied to the very complex mixtures for which no analytical technique is currently powerful enough to identify the chemical structure of the most abundant chemical species contained in these mixtures. This is typically the case for petroleum residues or biomass fractions.
Stochastic Molecular Reconstruction Approach
The stochastic molecular reconstruction approach models the molecular composition of a mixture from a set of Probability Distribution Functions (PDF) for molecular structural attributes (number of rings, chains length, number of chains, etc.). The governing principle of this approach is that any molecule can be considered as an assembly of molecular attributes. Each attribute is represented by a PDF which provides the probability of finding this attribute in a molecule of the mixture (Neurock et al., 1994; Trauth et al., 1994).
The concept of using PDF to describe complex mixtures has initially developed in the field of polymers for characterizing the molecular size distribution of polymers (Flory, 1936; Libanati, 1992). The application of the PDF was extended to the petroleum field by Klein’s group (Campbell and Klein, 1997; Neurock et al., 1994; Trauth et al., 1994) based on detailed analysis of heavy fractions done by Boduszynski (Boduszynski, 1988, 1987). Boduszynski showed that the structural properties of a petroleum fraction follow statistical distributions when plotted against the Atmospheric Equivalent Boiling Point (AEBP), the molecular weight or the number of carbon atoms. Therefore, Neurock et al. (1994) suggested the characterization of petroleum fractions by means of a set of PDF for molecular attributes. They developed a Stochastic Reconstruction (SR) method to generate an equimolar set of molecules from a set of PDF for molecular attributes which are characteristic of the mixture to be represented. In practice, each PDF is sampled via a MC procedure to identify the structural attributes of a molecule, as illustrated in Figure 7. The sequence of the PDF sampling steps is oriented by a building diagram, which is defined by the modeler. The selected attributes are randomly assembled to obtain the structure of the molecule. MC sampling of the set of PDF can be repeated N times so as to obtain a mixture of N molecules.
Figure 7. Illustration of the SR method for creating an asphaltene molecule (Neurock et al., 1994). 
The initial SR method imposes as an input a set of PDF parameters which are representative of the mixture to be reconstructed. However, the PDF parameters are either known from prior work, or obtained from detailed analytical methods, which are very time consuming. Therefore, Klein’s group (Trauth et al., 1994) proposed to include an optimization loop to estimate the PDF parameters from the overall properties as in Figure 8. To do this, each PDF is described by a parameterized mathematical law (Gaussian distribution, gamma distribution, exponential distribution, etc.). For a set of PDF parameters, a mixture of molecules is generated, and their overall properties are calculated. The calculated properties are then compared to the experimental data by means of an objective function (Hudebine and Verstraete, 2004). The objective function is then minimized by modifying the PDF parameters. The optimization is performed using a global optimization method, such as simulated annealing, genetic algorithms, or particle swarm optimization.
Figure 8. Flow diagram of the SR method with the optimization loop. (Trauth et al., 1994). 
The SR method is an ingenious technique to create a molecular representation of the complex mixture from the overall properties. To overcome the lack of molecular information, this approach supposes certain assumptions of similarity or uses high internal constraints that cannot be provided by the overall properties. However, the SR method also has disadvantages. The selection of the molecular attributes, the PDF and the building diagram are no trivial tasks to be performed by the modeler. The same mixture represented by different sets of molecular attributes and/or PDF may lead to different solutions. Secondly, this method is computationally very demanding because the optimization loop needs to generate a large number of molecules before reaching an accurate molecular representation (Hudebine and Verstraete, 2004; Van Geem et al., 2007; Wu and Zhang, 2010). Finally, SR is not always capable a fitting all mixture properties simultaneously due to the numerous constraints imposed by the PDF and building diagram (Hudebine and Verstraete, 2004).
4.1.3 Examples of Molecular Reconstruction Methods Developed at IFPEN
Given the limitations of lumped models, IFPEN has developed the molecularbased kinetic models that prompted the development of molecular reconstruction methods.
Reconstruction by Entropy Maximization (REM) (Hudebine and Verstraete, 2011, 2004; Hudebine, 2003) was the first method developed at IFPEN. The REM method is a deterministic molecular reconstruction approach that is based on the principle of maximum Shannon entropy, according to which an information entropy criterion must be maximized in order to obtain the optimal result (Shannon, 1948). The mole fractions of the predefined set of molecules are adjusted by maximizing the Shannon entropy criterion (E), given by Equation (9), subject to the mass balance (Eq. 10) and the analytical constraints (Eq. 11). The maximization of the criterion is performed by combining the Lagrange multiplier method and the conjugate gradient optimization method. This method was successfully applied to LCO gas oils (Hudebine, 2003), FCC gasolines (Hudebine and Verstraete, 2011) and steam cracking naphthas (Van Geem et al., 2007).$$E=\sum _{j=1}^{N}{x}_{j}\mathrm{ln}\left({x}_{j}\right)$$(9)with$$\sum _{j=1}^{N}{x}_{j}=1$$(10)and$${F}_{k}\left(\mathbf{x}\right)={P}_{k}\hspace{1em}k=1,K$$(11)where E represents Shannon’s information entropy, N is the number of molecules present in the predefined database of molecules, x _{ j } represents the mole fraction of molecule j, x is the vector of mole fractions x _{ j }, P _{ k } is the value of analytical property k, F _{ k } is the mixing rule for analytical property k, and K is the number of analytical properties or constraints.
In order to extend the REM method to more complex mixtures, IFPEN developed a twostep reconstruction (SRREM) method in which REM has been coupled to a SR method (Hudebine and Verstraete, 2004; Hudebine, 2003). The SR method is first applied to generate an equimolar database of molecules that are typical for a given mixture. The REM method is then used to improve the mixture properties of the previously generated set of molecules by modifying their molar fractions. The SRREM method overcomes the drawbacks of both reconstruction methods. The SR method is a very constrained method that creates an equimolar mixture, by generating molecules that are always typical of the studied mixture due to the use of specific molecular attributes and building rules based on analytical information and expert knowledge, but at final convergence, it only approaches the analytical data. In contrast, the REM method is a fast technique to obtain a mixture whose properties closely match the analytical data, but it needs a particularly welldefined initial set of molecules (Hudebine and Verstraete, 2004).
Two different SRREM algorithms were developed at IFPEN: a direct algorithm and an indirect SRREM algorithm (de Oliveira, 2013; de Oliveira et al., 2013a,c; Hudebine and Verstraete, 2004). The direct SRREM algorithm (Fig. 9a) consists of applying the SR and REM methods as two consecutives steps to reconstruct a mixture. Thus, an equimolar database of molecules is generated for each studied mixture which is then submitted to REM method to improve the mixture properties. However, as SR step of this reconstruction method uses a MC approach to generate its initial set of molecules, the direct SRREM still remains a computationallydemanding algorithm. The indirect SRREM algorithm (Fig. 9b) was therefore developed and may be applied when a large number of mixtures of a given type need to be represented. The idea behind this algorithm is to generate a large reference database of representative molecules for a type of mixture once and for all. Subsequently, only the REM step is employed, requiring only a few seconds of CPU time. For each mixture, the REM step starts from this reference database to obtain a new set of molecules whose properties are very close to those of the mixture to be represented. All synthetic mixtures will now contain the same molecules as the reference database, but with modified mole fractions.
Figure 9. Flow diagram of the direct SRREM algorithm (a) and indirect SRREM algorithm (b) for reconstruction of petroleum VR. (de Oliveira et al., 2013a). 
The SRREM method was developed for and validated on LCO gas oils (Hudebine and Verstraete, 2004; de Oliveira et al., 2012; López Abelairas et al., 2016) and was later extended to vacuum gas oils (CharonRevellin et al., 2011; Verstraete et al., 2004) and VR (de Oliveira et al., 2014, 2013ac; Verstraete et al., 2010).
As mentioned in Section 3.3, IFPEN developed a feedstock reconstruction method, called statistical reconstruction, in order to improve the AGO hydrotreating kinetic model. Statistical reconstruction (Hudebine et al., 2011; López García et al., 2010) is a deterministic molecular reconstruction approach that provides a matrix of structural lumps or pseudocompounds to represent the composition of middle distillates. The rows of the matrix contain 28 different chemical families, such as paraffins, naphthenes, aromatics, sulfur compounds, while the columns represent the molecule size expressed by carbon atoms (from 1 to 30). Each matrix element contains the mole fraction of the pseudocompounds which are estimated from four chemical analyses: MS, sulfur speciation, nitrogen speciation and simulated distillation. The first three analyses provide the overall mole fraction of the 28 chemical families, which are then distributed among the pseudocompounds based on the simulated distillation data by using two gamma distributions. These statistical distributions represent the variation of the size of molecules within each family by means of an increase in alkyl chain length. The first distribution defines the chain length for acyclic structures (paraffins, mercaptans, amines) and the second distribution describes the length of the side chain grafted on the core of the cyclic compounds. For each chemical family, the distribution is adapted with respect to its minimum and maximum number of carbon atoms and to the initial and final boiling point of the simulated distillation.
4.2 Simulation of Large Complex Reaction Networks
4.2.1 Network Generation and Network Reduction
Generation of reactions and reaction networks fall under the larger concept of ‘computerassisted chemistry’. This vast field of ‘computerassisted chemistry’ comprises a large variety of applications such as kinetics packages, network generation programs, tools for designing experiments, property estimation methods, structure elucidation software, structure visualization applications, computational chemistry, quantum chemistry, force field methods, computerassisted organic synthesis planners, and many more.
Concerning reaction networks, several automated network generation algorithms are available in the literature, both for computeraided organic synthesis and for kinetic and reactor modeling purposes. Automated network generation algorithms have been used in different applications:

hydrocarbon pyrolysis (Clymans and Froment, 1984; Hillewaert et al., 1988; Chinnick et al., 1988; Froment, 1991; Billaud et al., 1991; Dente et al., 1992; Di Maio and Lignola, 1992; Broadbelt et al., 1994, 1995; Susnow et al., 1997; DeWitt et al., 2000; Faulon and Sault, 2001; Kruse et al., 2001, 2002, 2003; Matheu et al., 2001, 2003a,b; Bounaceur et al., 2002; Grenda et al., 2003; Klein et al., 2006; Van Geem et al., 2006, 2008; Levine and Broadbelt, 2009; Vandewiele et al., 2012; Karaba et al., 2013; De Bruycker et al., 2015; Van de Vijver et al., 2015a,b),

hydrocarbon oxidation and combustion (Chevalier et al., 1990, 1992; ValdesPerez, 1992a; Ranzi et al., 1995; Blurock, 1995; Côme et al., 1996; Gaffuri et al., 1997; Glaude et al., 1997, 1998; Iyer et al., 1998; Warth et al., 2000; BattinLeclerc et al., 2000; Green et al., 2001; Heyberger et al., 2001; Matheu et al., 2001, 2003b; Németh et al., 2002; Ratkiewicz and Truong, 2003, 2006; Song et al., 2003; Katare et al., 2004; Buda et al., 2006; Pfaendtner and Broadbelt, 2008a,b; Hognon et al., 2012),

isomerization (Guillaume et al., 2003ac; Surla et al., 2004, 2011),

alkylation (Martinis and Froment, 2006),

olefin oligomerization and alkylation (Prickett and Mavrovouniotis, 1997c; Guillaume, 2006; Shahrouzi et al., 2008; Toch et al., 2015),

propane aromatization (Katare et al., 2004; Bhan et al., 2005),

catalytic reforming (Joshi et al., 1999; Klein et al., 2006; Wei et al., 2008; SoteloBoyás and Froment, 2009; Cochegrue et al., 2011),

catalytic cracking (Feng et al., 1993; Prickett and Mavrovouniotis, 1997a; Joshi et al., 1998; Dewachtere et al., 1999; Christensen et al., 1999; Moustafa and Froment, 2003; Froment, 2005; QuintanaSolórzano et al., 2007a,b, 2010; Xue et al., 2014),

hydrocracking (Baltanas and Froment, 1985; Baltanas et al., 1989; Vynckier and Froment, 1991; Quann and Jaffe, 1992, 1996; Quann, 1998; Martens and Froment, 1999; Mizan and Klein, 1999; Martens et al., 2000, 2001; Thybaut and Marin, 2003; Laxmi Narasimhan et al., 2003b, 2004; ChavarríaHernández et al., 2004; Jaffe et al., 2005; Klein et al., 2006; Kumar and Froment, 2007a,b; ChavarríaHernández et al., 2008; ChavarríaHernández and Ramírez, 2009; Vandegehuchte et al., 2012, 2015),

hydrotreating (Hou and Klein, 1999; Klein et al., 2006),

hydrogenation (Bera et al., 2011, 2012),

hydrogenolysis (ValdesPerez, 1994; ValdesPerez and Zeigarnik, 1997),

methanoltoolefins (Park and Froment, 2001a,b, 2004; Alwahabi and Froment, 2004; Froment, 2005; Kumar et al., 2013a,b),

methanoltogasoline (Jin and Froment, 2013),

FischerTropsch synthesis (Klinke and Broadbelt, 1999; Temkin et al., 2002; Lozano Blanco et al., 2006, 2008, 2011),

coke formation (Moustafa and Froment, 2003; QuintanaSolórzano et al., 2005),

biomass conversion (Rangarajan et al., 2010),

oxidative coupling of methane (Simon et al., 2007),

silicon hydride clustering (Wong et al., 2004),

hydrocarboxylation (Zeigarnik et al., 1997),

carbonylation (Bruk et al., 1998),

organic synthesis (Fontain and Reitsam, 1991),

biological processes (ValdesPerez, 1992b; Klein et al., 2002; Faeder et al., 2003, 2005; Li et al., 2004; Hatzimanikatis et al., 2004, 2005; Blinov et al., 2005a,b; Mayeno et al., 2005; Henry et al., 2007, 2010; Hill et al., 2008; Finley et al., 2009),

and many more.
This illustrates that automated network generation techniques are well suited for any type of reaction, as long as one can define the transformations that are typical of a particular system. Most of these algorithms generate an exhaustive list of the vast number of reactions, thus avoiding errorprone manual reaction network developments. Good overviews of some of these automated network generation methods for different applications can be found in Ugi et al. (1993), Prickett and Mavrovouniotis, (1997b), Tomlin et al. (1997) and Klein et al. (2006).
Ugi et al. (1993) tried to classify network generation algorithms and identified three types:

empirical methods,

semiformal methods

formal methods.
According to their classification, empirical methods are based on preestablished reaction libraries and expert systems. Examples of such databases were already underway in the 1960s (Vleduts, 1963) and continue to be developed (Ihlenfeldt and Gasteiger, 1996). Semiformal techniques are based on heuristic algorithms that generate reactions, either overall reactions or elementary steps, in a network from a few basic reaction types, but do not have selection procedures to remove improbable reactions from the set of conceivable solutions. Formal techniques generate the various elementary steps at the mechanistic level by taking into account the physical and chemical description of the molecule, and select the most probable elementary steps to be retained in the final network from all possible steps.
Once the reaction network is available, it can be used for simulation purposes. For the complex feedstocks, the reaction network may become very large and very difficult to manage along the simulation, however. In such a case, a direct simulation may not be possible without some reduction of the network. Reduction of a reaction network is not so easy. Complex mathematical tools are required that generally use the full knowledge of the reaction network, including the reaction rates in order to classify the reactions according to different scales of reactivity: ‘fast’ reactions and ‘slow’ reactions. During the development of the model, some qualitative information or orders of magnitude may be available, but quantitative information on the kinetics is rarely known. Hence, a priori reduction of a reaction network is usually based on empirical rules, such as limiting the maximum molecule size, excluding molecule types, ignoring reaction families, constraining the total number of species and/or reactions, selecting reactions on their reaction (free) enthalpy or on their presumed reaction rate, etc. A bad choice of these ‘reduction rules’ can easily lead to biased reaction networks, resulting in incorrect predictions and behaviors during simulation.
Rigorous reduction techniques aim to retain as much as possible the detail of the original reaction network. When the kinetic equations between species are known and sets of species are in thermochemical equilibrium, one can group the species in equilibrium into lumps, and write a new reduced kinetic network between lumps that is exactly equivalent to the original network (Wei and Kuo, 1969; Li and Rabitz, 1989; Vynckier and Froment, 1991; Cochegrue et al., 2011). In this sense, this is strictly an exact lumping technique as defined by Wei and Kuo (1969), since the composition of each lump is known through the equilibrium relations. Hence, when experimental data shows that the species of a lump are in thermochemical equilibrium, they proposed an a posteriori lumping scheme, also known as a late lumping scheme, which groups species at thermodynamic equilibrium inside a single lump. It is therefore a rigorous network reduction technique that allows the production of a simplified network, which can then be used in reactor simulations (Vynckier and Froment, 1991). Indeed, using such a rigorous network reduction method avoids the problems related to the diversity of feedstocks, unlike a priori relumping of species, also known as early lumping, which is generally a ‘blind’ lumping approach devised to cope with the lack of componentbased analytical data and where empirical kinetics do not account for the composition of the lumps. Reduction techniques may also use stochastic network generation (Faulon and Sault, 2001; Mizan et al., 1998; Shahrouzi et al., 2008), or seeding and deseeding techniques (Joshi, 1998; Joshi et al., 1999; Watson et al., 1996).
Reaction network reduction techniques can be divided into five categories:

global reduction,

response modeling,

chemical lumping,

statistical lumping,

detailed reduction (Frenklach, 1987).
Global reduction techniques transform a complete reaction scheme into a small number of overall reaction steps. These techniques comprise adhoc methods such as empirical fitting, reduction by approximations, and lumping. Global reduction techniques are problemspecific and cannot be generalized. Response modeling techniques consist of mapping model responses (species concentrations) and model variables (initial boundary conditions, rate parameters, transport properties) through functional relationships, typically combinations of simple algebraic functions such as polynomials, by using either experimental data or computer experiments. As for global reduction techniques, response modeling solutions are problemspecific since they require data to construct the algebraic functions. The chemical lumping method is guided by similarity in the chemical structure or chemical reactivity of the reacting species, i.e. the reactions in the entire network have essentially the same rate parameters, as is the case in polymerization reactions. The statistical lumping technique is used when growth processes can be described by the Smoluchowski equation. Finally, the detailed reduction technique consists of identifying and removing noncontributing reactions. An effective reduction strategy is to compare the individual reaction rates with the rate of a chosen reference reaction, which can be a given ratelimiting step or the fastest reaction. Consequently, the detailed reaction reduction approach is a general technique that is applicable to any reacting system.
4.2.2 Rate Equations and Rate Parameters
The rate equations for mechanistic kinetic models are straightforward. Indeed, due to its fundamental approach that explicitly describes the elementary steps, the rate of each reaction follows Guldberg and Waage’s law of mass action. In molecular kinetic models, the reactions are viewed as moleculetomolecule transitions, without accounting for the detailed reaction mechanism and its reaction intermediates. Hence, the rate equations may need to include the indirect effects of the reaction intermediates by using powerlaw type rate equations or LangmuirHinshelwood type rate equations. This type of rate equation can be derived using the above described QSS approximation, RDS approach, or one of their variants.
The number of rate parameters in a large reaction network will be huge. It is therefore almost impossible to obtain a value for each rate parameter through experimental work or through theoretical calculations. Hence, in almost all detailed kinetic models, rate parameters are typically defined for homologous series of reactions. The idea is that the reaction rate is determined by the transformation that happens at the ‘reacting center’ of the molecule, i.e. by the groups that are directly involved, while the rest of the molecule is only a spectator. This will be considered to be particularly true in large molecules. It should be mentioned, however, that this remains an important simplification of the reality, since actual molecules are continuously vibrating, bending and rotating in such ways that the ‘nonreacting’ moieties may still influence the reactivity.
For homogeneous acidcatalyzed reactions, Brønsted (Brønsted and Pedersen, 1924; Brønsted, 1928) observed that strongly exothermic reactions have a low activation energy and derived from his experimental work that the catalytic activity is directly related to the ionization constant K _{ a }. The Brønsted equation can be written as:$$k={\alpha}_{1}\xb7{\left({K}_{a}\right)}^{{\alpha}_{2}}$$(12)but is generally given in its linearized form:$$\mathrm{ln}\left(k\right)={\alpha}_{1}+{\alpha}_{2}\xb7\mathrm{ln}\left({K}_{a}\right)$$(13)
The Brønsted equation can also be written differently. Introducing into this equation the van ‘t Hoff relation for the ionization equilibrium constant and the Eyring equation for the rate constant leads to:$$\mathrm{ln}\left(\frac{{k}_{B}T}{h}\right)\xb7\frac{{\Delta G}^{\ne}}{\mathrm{RT}}={\alpha}_{1}{\alpha}_{2}\xb7\frac{{\Delta G}_{a}}{\mathrm{RT}}$$(14)
After rearranging, one obtains:$$\u2206{G}^{\ne}=\frac{{\alpha}_{1}\xb7\mathrm{RT}}{\mathrm{ln}\left({k}_{B}T/h\right)}+\frac{{\alpha}_{2}}{\mathrm{ln}\left({k}_{B}Th\right)}\xb7\u2206{G}_{a}$$(15)or,$$\u2206{G}^{\ne}={\beta}_{1}+{\beta}_{2}\xb7\u2206{G}_{a}$$(16)
This form of the Brønsted equation shows that the Gibbs free energy of activation for a reaction ΔG^{≠} is directly related to its Gibbs free energy of reaction ΔG_{a}, i.e. the Gibbs free energy released by the reaction. The Brønsted equation is therefore historically the first example of a ‘linear Gibbs energy relation’ (the IUPAC recommended term), more commonly known as ‘Linear FreeEnergy Relationships’ (LFER).
One of the most famous LFER is the Hammett equation used in organic chemistry. Hammett (1937) found that, for a large set of reactions involving benzoic acid derivatives with meta and parasubstituents, the ratios of the rate constants were related to the ratios of the ionization equilibrium constants:$$\frac{k}{{k}_{0}}={\left(\frac{K}{{K}_{0}}\right)}^{{\rho}_{H}}$$(17)
The Hammett equation is again generally written in a linearized form:$$\mathrm{ln}\left(\frac{k}{{k}_{0}}\right)={\rho}_{H}\xb7\mathrm{ln}\left(\frac{K}{{K}_{0}}\right)={\rho}_{H}\xb7{\sigma}_{H}$$(18)or,$$\mathrm{ln}\left(k\right)=\mathrm{ln}\left({k}_{0}\right)+{\rho}_{H}\xb7{\sigma}_{H}$$(19)where k is the reaction rate constant of a substituted reactant, k_{0} is the reaction rate constant of the unsubstituted reactant, σ_{H} is the substitution constant, which represents the difference of the Gibbs free energies of ionization and accounts for field, inductive, and resonance effects, and ρ_{H} is the sensitivity constant, which only depends on the type of reaction, not on the type of substituent. The original Hammett equation was later refined, leading to the Taft equation, the SwainLupton equation, the GrunwaldWinstein equation, the YukawaTsuno equation, etc.
The origin of such linear freeenergy relationships was first derived by Evans and Polanyi (1935, 1936, 1938), who approximated the PES by straight lines. This approximation leads to a theoretical justification of a LFER, a linear relation between the activation free energy ΔG^{≠} and the Gibbs free energy of the reaction ∆G_{R}:$$\u2206{G}^{\ne}={\beta}_{1}+{\beta}_{2}\xb7\u2206{G}_{R}$$(20)where β_{1} and β_{2} depend on the slopes of the lines and the initial distance between reactants and products. This relation can be intuitively understood for endergonic elementary reactions, since the free energy barrier for activation has to increase when the free enthalpy of reaction increases. Within a series of closely related atomtransfer reactions, Evans and Polanyi (1938) also showed that the activation energies depend linearly on the reaction enthalpies. This ‘Linear Activation Energy–Reaction Energy Relationship’ or linear energy relationship is generally referred to as the EvansPolanyi relationship and has the following form:$${E}_{a}={E}_{0}+\alpha \xb7\u2206{H}_{R}$$(21)
Evans and Polanyi applied this relationship to organic reactions, but the EvansPolanyi relationships were later shown to apply to many different types of atom transfer reactions. Semenov (1954, 1958) studied a large number of atom and radical reactions, extended the EvansPolanyi relationships to chain reactions, and suggested a value of 48 kJ.mol^{−1} for E_{0}, while α = 0.25 for exothermic reactions and α = 0.75 for endothermic reactions. Temkin was the first to apply the EvansPolanyi relationships to heterogeneous catalysis, together with Mochida and Yoneda (1967ac).
Marcus (1968) approximated the PES by two intersecting parabolas whose minima correspond to the free energy levels of the reactants and of the products respectively. The free enthalpy barrier now corresponds to the point of intersection of the two parabolas, and its values is given by:$$\u2206{G}^{\ne}=\frac{\lambda}{4}{\left(1+\frac{\u2206{G}_{R}}{\lambda}\right)}^{2}=\frac{{\left(\lambda +\u2206{G}_{R}\right)}^{2}}{4\lambda}=\frac{\lambda}{4}+\frac{\u2206{G}_{R}}{2}+\frac{{\left(\u2206{G}_{R}\right)}^{2}}{4\lambda}=\u2206{G}_{\mathrm{int}}^{\ne}+\frac{\u2206{G}_{R}}{2}+\frac{{\left(\u2206{G}_{R}\right)}^{2}}{16\u2206{G}_{\mathrm{int}}^{\ne}}$$(22)where λ is called the ‘reorganization energy’ and depends on the curvature of the parabolas and the distance between the energy minima of reactants and products. The group λ/4, or , is called the ‘intrinsic barrier’, i.e. the free energy barrier for a reaction whose free energy of reaction equals zero. Marcus initially developed this equation to correlate solution reaction rates for electron transfer, but this relation has been proven useful for atom, proton, hydride, and group transfers as well. The above Marcus equation is sometimes referred to as a Quadratic FreeEnergy Relationship (QFER). When the reorganization energy λ is much larger than the Gibbs free energy of reaction, the curvature becomes negligible and the Marcus equation reduces to a LFER.
Numerous attempts have been made to extend such linear freeenergy relationships or EvansPolanyi relationships. In many cases, the developed relationships are arrived at by means of a statistical approach that relates the rate constant k or the activation energy E_{a} to various physicochemical properties of the reactants and/or products. Furthermore, many of these empirical or semiempirical equations no longer correlate the activation free energy ΔG^{≠} or activation energy E_{a} to other freeenergy terms or energy terms. Hence, one should no longer refer to such relations as linear freeenergy relationships or linear energy relationships, but rather use the terms Quantitative StructureActivity Relationships (QSAR), Quantitative StructureReactivity Correlations (QSRC) or simply the term Correlation Analysis. The basic premise behind these QSAR equations is that, for a homologous series of reactions, the rate constants can be correlated to some molecular property. The simplest form of a QSAR is the simple linear relationship:$$\mathrm{ln}\left({k}_{i}\right)=\alpha +\beta \xb7{\theta}_{i}$$(23)where k_{i} is the rate constant of molecule i for a given reaction, θ_{i} is the reactivity index of molecule i, and α and β are the correlation coefficients for a given reaction family. In many cases, a multilinear QSAR is developed for each reaction family:$$\mathrm{ln}\left({k}_{i}\right)={\beta}_{0}+\sum _{j=1}^{J}{\beta}_{j}\xb7{\theta}_{\mathrm{ji}}$$(24)where k_{i} is the rate constant of molecule i for a given reaction, θ_{ji} is the reactivity index j for molecule i, and the J+1 parameters β_{j} are the correlation coefficients for a given reaction family. Many other forms of QSAR equations (quadratic, etc.) can be found, and some theoretical basis is given in Wold and Sjöström (1978).
In the literature, many examples of QSAR equations can be found, especially in pharmaceutical research, toxicology, and biological sciences, either with or without a theoretical basis. To illustrate the various levels of theoretical background, some examples of QSAR equations for applications in the refining industry are given here. Neurock (1992) derived a linear QSAR that correlated the rate constant for the dealkylation of alkyl aromatics to the heat of reaction, and a linear QSAR that correlated the rate constant for aromatics hydrogenation to the πelectron density. Korre (1995) developed a multilinear QSAR for the rate constant for aromatics hydrogenation to the heat of reaction, the number of hydrogen atoms that were incorporated, the number of aromatic rings, and the number of naphthenic rings. Froment et al. (1994) developed a type of a groupcontribution method for the rate constants for the hydrodesulfurization of all benzothiophene and dibenzothiophene isomers, by accounting for the steric and electronic effects of the various alkyl substituents and their positions around the sulfur atom. Liguras and Allen (1989) correlated the overall rate constant for catalytic cracking of paraffins to the total number of carbon atoms, and the number of carbon atoms of each type (primary, secondary, tertiary atoms). López García et al. (2010) used, depending on the reaction family, a value for the reaction rate that was either constant or a function of the carbon number. For the hydrocracking of vacuum gas oils, Stangeland (1974) proposed a nonlinear correlation of the rate constant with the boiling point temperature of the fraction or the component.
In conclusion, LFER and QSAR are empirical correlations observed between rate parameters and equilibrium constants, rate coefficients and/or structural parameters, usually among a series of reactions within a given mechanism (Connors, 1990). They are an empirical but quantitative way of ‘systematizing’ the similarity of reactions in a welldefined domain of applicability. They imply that the steric and entropy terms remain almost constant through a homologous series of reactions. These relations are also referred to as ‘extrathermodynamic relationships’, since nothing makes these correlations necessary from a theoretical point of view (Chapman and Shorter, 1978). Despite this, they are relatively common, due to the fact that the structure of a molecule outside of the ‘reacting center’ has only a limited influence on the reaction rate, and are very useful in chemistry and in chemical engineering applications.
4.2.3 Reaction Network Simulation
Inserting a detailed kinetic model into a steady state onedimensional reactor model generally yields a set of coupled Ordinary Differential Equations (ODE) (for ideal plug flow reactors or batch reactors) or a set of coupled algebraic equations (for ideal continuous stirred tank reactors). When the detailed kinetic model is inserted into transient and/or more detailed reactor models, coupled PDE and/or AlgebroDifferential Equations (ADE) arise. Such a set of equations can then be solved using classic numerical tools. This meanfield approach is the classic way to perform reactor simulations.
In some cases, it is not possible to arrive at a complete description of the reaction network. In other cases, the meanfield approximation breaks down, as can happen for surface reactions with energetic heterogeneities on the surface, surface diffusion, and/or attractiverepulsive interactions between adsorbed species. In such cases, a change in point of view is needed. One alternative way to simulate large systems is to change the basic description of the problems. Instead of following a macroscopic description of the mixture (concentrations, pressures, etc.), the evolution of a population of molecules and their transformations are described. By analogy to fluid mechanics, this is similar to the distinction between an Eulerian description (behavior of properties function pressure, specific gravity, velocity, etc.) and a Lagrangian description (individual motion and properties of a molecule).
An important class of models involves the description of random phenomena. The probabilistic character of chemical reactions has historically been used for the description of reacting ideal gases using the notion of colliding molecules and energy levels to determine whether a reaction will occur or not. In these applications, this description was not used for simulation purposes, but for the description of the microscopic events that lead to reaction. The theoretical background of the probabilistic description of reaction events is based on the Chemical Master Equation (CME). This equation provides the timeevolution of a system of countable states where switching between states is treated probabilistically. From a probabilistic point of view, the following set of differential equations is obtained:$$\frac{\partial}{{\partial}_{t}}\mathcal{P}\left(\stackrel{\u20d7}{X},t{\stackrel{\u20d7}{X}}_{0},{t}_{0}\right)=\sum _{j=1}^{M}\left\{\mathcal{P}\left(\stackrel{\u20d7}{X}{\stackrel{\u20d7}{v}}_{j},t{\stackrel{\u20d7}{X}}_{0},{t}_{0}\right)\times {a}_{j}\left(\stackrel{\u20d7}{X}{\stackrel{\u20d7}{v}}_{j}\right)\mathcal{P}\left(\stackrel{\u20d7}{X},t{\stackrel{\u20d7}{X}}_{0},{t}_{0}\right){\times a}_{j}\left(\stackrel{\u20d7}{X}\right)\right\}$$(25)
This equation is difficult to solve and has no analytical solution, except for some very simple academic cases. The development of computers and appropriate algorithms allowed the numerical simulation of those reaction events. This approach was formalized by Gillespie (1976), resulting in an algorithm (Stochastic Simulation Algorithm – SSA) that is as simple as the theoretical description (CME) is complex.
Gillespie (1976) introduce the notion of a propensity function that represents the probability that a reaction occurs. By rewriting the equations, he proved that the evolution between two random events can be decoupled by the random generation of two numbers, the first giving the time of the next reaction, and the second choosing the reaction that will occur according to the propensity function. This algorithm is classically called the SSA and is shown in Figure 10.
Figure 10. Gillespie’s SSA (Gillespie, 1976). 
Each realization of the algorithm gives a particular simulation. The average of a large number of reactions converges to the deterministic simulation of a mixture with the same macroscopic properties (Fig. 11). The size of the initial population has an important effect on the speed of the convergence (Shahrouzi, 2010; de Oliveira et al., 2013a).
Figure 11. Results of stochastic and deterministic simulations for the reaction with c_{1} = 0.5 s^{−1}, c_{2} = 0.1 s^{−1}, starting from X_{A}(0) = 100, X_{B}(0) = 0. 
With an a priori determined reaction network, there are no advantages of using the SSA to integrate the continuity equations. For very large reaction networks, however, the memory requirements for the SSA are generally lower than classic numerical integration packages.
4.2.4 Examples of Complex Reaction Network Simulations Developed at IFPEN
Classic SingleEvent Kinetic Modeling
The singleevent kinetic modeling methodology has been developed at the Ghent University by Froment and coworkers. The ‘single event’ concept links the reactivity of a molecule to its geometry and to a limited number of intrinsic kinetic parameters, which depend only on the type of molecule and the type of reaction. Using the concept of ‘single events’ is an efficient way of modeling reactions without introducing reductive assumptions about the composition of the feed (Froment, 1999, 2005).
The singleevent approach has been applied to many acidcatalyzed processes and to metalcatalyzed processes. This modeling approach has been applied for acidcatalyzed processes such as:

isomerization (Guillaume et al., 2003a,b,c; Surla et al., 2004, 2011),

alkylation (Martinis and Froment, 2006),

olefin oligomerization (Guillaume, 2006; Shahrouzi et al., 2008; Toch et al., 2015),

methanoltoolefins (Park and Froment, 2001a,b, 2004; Alwahabi and Froment, 2004; Froment, 2005; Kumar et al., 2013a,b),

catalytic reforming (SoteloBoyás and Froment, 2009; Cochegrue et al., 2011),

catalytic cracking (Feng et al., 1993; Dewachtere et al., 1999; Beirnaert et al., 2001; Moustafa and Froment, 2003; Froment, 2005; QuintanaSolórzano et al., 2005, 2007a,b, 2010; Xue et al., 2014),

hydrocracking (Baltanas and Froment, 1985; Baltanas et al., 1989; Vynckier and Froment, 1991; Svoboda et al., 1995; Schweitzer et al., 1999; Martens and Froment, 1999; Martens et al., 2000, 2001; Martens and Marin, 2001; Thybaut et al., 2001, 2009; Thybaut and Marin, 2003; Laxmi Narasimhan et al., 2003a,b, 2004, 2006, 2007; ChavarríaHernández et al., 2004, 2008; Valéry et al., 2007; Kumar and Froment, 2007a,b; ChavarríaHernández and Ramírez, 2009; Mitsios et al., 2009; Choudhury et al., 2010; Vandegehuchte et al., 2012, 2015).
Singleevent kinetic modeling has also been applied to FischerTropsch synthesis (Lozano Blanco et al., 2006, 2008, 2011).
The singleevent methodology has already been extensively described in the literature (Vynckier and Froment, 1991). At IFPEN, this methodology has been applied to:

paraffin isomerization (Guillaume et al., 2003a; Surla et al., 2004, 2011),

olefin oligomerization (Guillaume, 2006; Shahrouzi et al., 2008),

catalytic reforming (Cochegrue et al., 2011),

hydrocracking (Schweitzer et al., 1999; Valéry et al., 2007; Mitsios et al., 2009; Guillaume et al., 2011).
In the singleevent methodology, the level of detail is very high, since the transformations on the acid sites concern the elementary steps between carbenium ions. The reactions considered depend on the process. For bifunctional catalysts (isomerization, catalytic reforming, hydrocracking), the metal phase reactions concern the hydrogenation and dehydrogenation steps, while the acid phase reactions concern the protonation of olefins and the deprotonation, Hydride Shift (HS), Methyl Shift (MS), ethyl shift, alkyl shift, intraring alkyl shift, Protonated CycloPropane (PCP) isomerization, Protonated CycloButane (PCB) isomerization, and βscission of carbenium ions. For catalytic cracking, the same acid phase reactions are supplemented with protolytic cracking reactions of paraffins (Fierro et al., 2001, 2002) and hydride transfer between olefins and carbenium ions.
The network generation is based on a semiformal network generation algorithm in which reaction rules describe the reaction pathway occurring on the catalyst as a function of the reactant and its configuration. These reaction rules are iteratively used by an algorithm that handles a numerical description of molecules and of reactions and generates in this way the exhaustive network of reactions that includes all reaction intermediates. The full description of the molecule and reaction coding is described in Surla et al. (2011). Figure 12 shows the reaction network for the catalytic reforming of nhexane to illustrate the level of detail included in this network generation algorithm.
Figure 12. Reaction network for the catalytic reforming of nhexane (without accounting for hydrogenolysis reactions). 
Once the network is generated, a reaction rate is attributed to each reaction. Stability considerations concerning the carbenium ions involved in the reactions, and the further decomposition of the elementary steps involving similar activated complexes allow the writing of reaction rates that are constrained by the concept of ‘single events’. For an elementary step, the activation entropy term in the Eyring equation is separated into its various contributions (Benson, 1976). For a homologous series of reactions in the network, the translational, vibrational and electronic contributions are considered to be the same, but the rotational contributions strongly depend on the symmetry of the molecule (Willems and Froment, 1988a,b). Hence, the entropy contributions due to the symmetry of the molecule and to its chirality are grouped into a factor and set apart from the other entropy contributions, which are grouped into what is termed the intrinsic activation entropy (Eq. 26).$$k=\frac{{k}_{B}T}{h}\mathrm{exp}\left(\frac{T\xb7\u2206{S}^{\circ \ne}\u2206{H}^{\circ \ne}}{\mathrm{RT}}\right)=\frac{{\sigma}_{\mathrm{glob}}^{r}}{{\sigma}_{\mathrm{glob}}^{\ne}}\times \frac{{k}_{B}T}{h}\mathrm{exp}\left(\frac{T\xb7\u2206{S}_{\mathrm{int}}^{\circ \ne}\u2206{H}^{\circ \ne}}{\mathrm{RT}}\right)$$(26)
The first factor is a ratio of global symmetry numbers and is termed the number of single events, while the remaining part is termed the intrinsic rate coefficient:$$k=\frac{{\sigma}_{\mathrm{glob}}^{r}}{{\sigma}_{\mathrm{glob}}^{\ne}}\times \frac{{k}_{B}T}{h}\mathrm{exp}\left(\frac{T\xb7\u2206{S}_{\mathrm{int}}^{\circ \ne}\u2206{H}^{\circ \ne}}{\mathrm{RT}}\right)={n}_{e}\times \stackrel{\u0303}{k}$$(27)where k is the rate coefficient of the elementary step, is the intrinsic rate coefficient of the elementary step, also termed the singleevent rate coefficient, n_{e} is the number of single events.
The intrinsic rate coefficient of the elementary step only depends on the nature of the involved ions, and of the type of reaction. Thus, the reaction rates are parameterized by a limited number of parameters (Surla et al., 2011). For example, the PCP isomerizations are fully described by four kinetic parameters although hundreds of thousands of PCP isomerizations can occur in complex mixtures. Additional thermodynamics constraints and assumptions allow the further reduction of the number of kinetic parameters (Surla et al., 2011; Cochegrue et al., 2011).
The generated reaction network can be simulated without any further assumptions, as shown by Verstraete (1997). Since the singleevent approach generates all possible reaction surface intermediates, the resulting reaction networks are huge. Network reduction techniques can be advantageously applied on singleevent reaction networks in order to reduce the number of kinetic equations of the network and to solve the continuity equations. For some of the processes (catalytic reforming, hydrocracking), the hydrogenation/dehydrogenation equilibrium is easily established on the metal phase. At the same time, experimental data obtained on various feedstocks and at various contact times demonstrated that the paraffins with the same number of carbon atoms and the same number of branches are in thermodynamic equilibrium with each other (Verstraete, 1997; Schweitzer et al., 1999). This leads to the conclusion that isomerization reactions without change in the branching degree (HS, MS, ethyl shift) are very fast compared to isomerization reactions with a change in the number of branches (PCP), cracking reactions and cyclization/ring opening reactions. Cyclic compounds with the same number of carbon atoms and the same number of branches are also in equilibrium with each other (Verstraete, 1997), indicating that ring contraction or expansion reactions are very fast compared to branching isomerization reactions. With this experimental information, it is possible to reduce the reaction network by means of an a posteriori relumping technique without losing any information (Cochegrue et al., 2011), as described by Vynckier and Froment (1991). With this technique, molecules are collected into ‘groups of isomers’ according to their carbon number and their degree of branching. Inside a group of isomers, the molecules are considered to be in equilibrium because of the very fast interconversion reactions, which allows the determination of the composition within each group. The size of the reduced reaction network is much smaller, since it contains between 10 and 1000 times fewer species and fewer reactions, while the introduction of a set of multiplicative ‘lumping coefficients’ allows maintaining the same set of fundamental rate parameters for the relumped reactions in the reduced network.
Using the properties of these lumping coefficients and the recurrence relations between the lumping coefficients even allowed writing these a posteriori relumped kinetic models without the full generation of the network (Valéry, 2002; Valéry et al., 2007; Mitsios et al., 2009; Guillaume et al., 2011). This socalled lateral chain decomposition method divides the various components (molecules, intermediates and activated complexes) into lateral chains and activated zones (Fig. 13). Detailed inspection of these various parts allowed the derivation of a set of recursive formulae that calculate the various structural and thermodynamic properties that are needed to determine the lumping coefficients. It should be stressed that this lateral chain decomposition technique is exactly equivalent to the original calculation method based on the explicitly generated reaction network, as long as Benson’s group contribution method is used to obtain the necessary thermodynamic data in both approaches (Benson and Cohen, 1993). With this lateral chain decomposition method, much larger reaction networks with higher carbon numbers can be generated (up to approximately 50 to 60 carbon atoms with the current version of the algorithm) very easily and almost instantly. The main advantage of these recursive formulae resides in the fact that application of the singleevent modeling methodology to very large reaction networks and very high degrees of branching is highly simplified, both in comparison to the explicit generation of the entire reaction network and to other approaches in the literature (Guillaume et al., 2011). This lateral chain decomposition method is quite flexible and can easily be extended to more general lateral chains, by introducing new insertion patterns and adapting the recursion formulae accordingly.
Figure 13. Reaction pathway and activated complex (a) for a β scission step and (b) for a PCP isomerization step. L_{1} to L_{5} are substituents which are part of the activated zone, while A and B are lateral chains, which do not participate in the elementary step (Guillaume et al., 2011). 
After relumping, the reduced reaction network for catalytic reforming contains almost 100 lumps and about 500 reactions, while the hydrocracking network contains a few thousand reactions between a few thousand lumps, whose evolution needs to be simulated. In both cases, the reactor simulations are performed by numerically integrating the continuity equations, which are ODE, by means of the LSODE package (Hindmarsch, 1980, 1983; Petzold, 1983).
In conclusion, refining processes based on acid or bifunctional catalysts, such as isomerization, catalytic reforming, catalytic cracking, and hydrocracking, involve a tremendous number of species (reactants, intermediates and products) and reactions, both of which may rapidly exceed several thousand for industrial feedstocks. The singleevent method represents an extremely elegant approach for the modeling of such processes, and also offers several advantages. First of all, it is based on a fundamental understanding of the chemistry and uses an exhaustive computergenerated reaction network that lists all reactions and reaction intermediates by means of a set of simple rules. Secondly, by applying the singleevent concept and by introducing detailed thermodynamic constraints between the rate parameters, the huge reaction networks can be simulated using only a limited number of independent fundamental kinetic parameters that are feedindependent. Since these rate coefficients are independent of the number of carbon atoms in the molecules, it is possible to identify them through studies on model molecules and to predict the behavior of complex feedstocks from experimental data that can be acquired more easily. A third characteristic is that the singleevent methodology also comes with a rigorous latelumping network reduction method that uses experimentallyverified and reasonable assumptions on the chemical equilibria to regroup the species at equilibrium into ‘groups of isomers’ or lumps. Under these assumptions, the substantially reduced relumped reaction network is strictly equivalent to the detailed network, while the use of lumping coefficients allows maintaining the same fundamental kinetic parameters to simulate the complete system without any loss of information. Finally, a lateral chain decomposition method has been developed for the singleevent methodology that allows the direct calculation of the lumping coefficients for the relumped reaction network by means of recursive formulae. This allows the application of the singleevent methodology to very large reaction networks and very high degrees of branching without explicitly generating the entire reaction network.
Network Reduction Techniques for SingleEvent Kinetic Modeling of Olefin Oligomerization
At IFPEN, the singleevent methodology was also applied to the oligomerization of light olefins. Two additional problems arose, however. The first one concerns the generation of the reaction network. Unlike the reaction networks for catalytic reforming, hydrocracking, or catalytic cracking, the reaction network for olefin oligomerization will theoretically never end, since light olefins can be added onto any olefin of any size. Hence, no exhaustive reaction network can be generated and, in practice, the reaction network is therefore limited by specifying a maximum number of carbon atoms. The second problem concerns the a posteriori relumping that was used previously. The experimental data for olefin oligomerization shows that thermodynamic equilibrium is not reached within the molecules with the same carbon number and the same degree of branching. It is therefore not possible to reduce the detailed reaction network to a manageable number of ODE to allow a reactor simulation to occur within a reasonable time. However, the concept of ‘single events’ still allows the description of the elementary steps with a small number of intrinsic parameters (Shahrouzi et al., 2008).
To resolve these two problems (no exhaustive reaction network and no tractable set of differential equations), stochastic simulation methods were applied. Several approaches were compared to investigate various network limitation methods. First, the full reference reaction network was generated by allowing species with up to 16 carbon atoms using a deterministic network generation program, which exhaustively enumerates all species and elementary steps. The resulting reaction network contained 279,045 species and 1,729,586 elementary steps. Every additional carbon number approximately leads to a threefold increase in the number of species and in the number of reactions. This C_{16} reaction network was simulated twice, once using a deterministic numerical integration using LSODE and once by applying Gillespie’s SSA. As expected, both approaches gave the same results.
A first network reduction approach consisted of limiting the deterministic network generation algorithm by a rulebased approach, such as putting an upper limit on the number of species and/or reactions. In this example, the network generation algorithm was stopped once the first 100,000 species (about 1/3 of the full reaction network) were generated. The resulting network is the Deterministic Limited Reaction Network (DLRN).
A second network reduction approach consisted of constructing a stochasticallygenerated, but limited reaction network. As in the previous methods, the reactor simulation is performed by the SSA, but starting from an ‘empty’ reaction network. The algorithm therefore has to generate all the possible reactions for the components present in the feed. Each of these reactions is now stored into the ‘initial’ reaction network. The SSA will then select one of these reactions and create the product species. If the product species are new to the current simulation, their possible reactions are generated and stored in the reaction network, while the SSA selects the next reaction. This procedure is repeated until the maximum storage limit of the reaction network is attained, in our case the first 100,000 species. From this point onwards, the simulation is continued on the stored network. For the next simulation, a new network is generated and the new simulation is based on this new network. In conclusion, this approach integrates the network generation into the SSA. Hence, the network is explored along the most probable reaction routes. Moreover, when several simulations are performed on the same feed, the average network combines all random networks and is therefore larger than each individual network. The strength of this network reduction technique resides in the fact that the exploration of the full reaction network is driven by probability considerations (exploring the more probable routes) and by averaging (combining several limited reaction networks).
Figure 14 schematically depicts the two network reduction techniques in comparison to the reference reaction network. Shahrouzi (2010) investigated the impact of these two network reduction methods on the simulation results. All simulations were performed using the SSA, and the Deterministic Full Reaction Network (DFRN) was used as a reference. The stochastic simulation of the DFRN represents the reference simulation, and a single simulation starting from 100,000 feed molecules took approximately 56.6 hours of CPU time on a dualcore CPU running at 2.33 GHz (Shahrouzi, 2010). The first rulebased network reduction technique provided a DLRN that contained the first 100,000 species and their reactions. During its stochastic simulation, the DLRN was loaded into memory, leading to a memory requirement that was approximately 10 times lower. A single stochastic simulation took approximately 4.5 hours of CPU time, which is 12.5 times faster than the stochastic simulation of the DFRN reference network. The second network reduction technique created a Stochastically Limited Reaction Network (SLRN). During the stochastic simulation, the SLRN grows into memory until 100,000 species and their reactions are stored. Again, the memory requirements are approximately 10 times lower than for the DFRN. The stochastic simulation itself also took approximately 4.5 hours of CPU time. The results of the three simulations are compared in Figure 15. For the first reduction technique, the results of the DLRN are quite different compared to those for the DFRN, not only for the reaction intermediates (the carbenium ions), but also for the final products (the olefins). The quality of the simulation with the deterministic reduced reaction networks was shown to highly depend on the reduction rules that are used during the network generation (total number of species, total number of reactions, rank of the products, reaction (free) enthalpies, presumed reaction rate, etc.). Moreover, by selecting a maximum number of species, the resulting network also depends on the processing order of the network generation algorithm. For the second reduction technique, the results for the SLRN agree well with the simulations of the DFRN. A very good agreement is observed for the final products, even though some discrepancies arise for the high carbonnumber carbenium ion intermediates (Shahrouzi, 2010). This illustrates that the SLRN allows the capturing of all the features of the full deterministic reaction network with a tenfold reduction in computer memory and in CPU time. This is mainly due to the fact that the full reaction network is explored based on reaction probability considerations along the more probable reaction pathways.
Figure 14. Representation of the DFRN (with 47 species), the DLRN (with only the first 20 species), and SLRN (with only the 20 species created by the ‘high probability’ pathways). 
Figure 15. Mass percent product distributions of (a) olefins and (b) carbenium ions for the oligomerization of C4 olefins up to C16 at 67% of conversion. Feed composition is 50% isobutene, 25% 1butene and 25% 2butene. The initial molecular population is 100,000. In the two limited reaction networks, the PCP steps are not considered and the protonation/deprotonation equilibrium is also considered in the DLRN. 
The SLRN network reduction technique also has an additional advantage: when several stochastic simulations are performed, the use of a new SLRN during each simulation also allows the coverage of a bigger network than if the first SLRN was reused in all simulations. Hence, an alternative network reduction technique can also be suggested here: if the size of the cumulative reaction network remains acceptable, it can be exported and reused with deterministic simulation tools. This approach would allow to stochastically generate reaction networks based on the reaction probabilities.
In the examples above, the chosen implementation still requires the provision of memory space to store the reaction network. An alternative technique consists of performing stochastic simulations without a pregenerated network by using a Markovian approach. Consider a population of molecules without a predefined reaction network. If one can reconstruct for each molecule the different reactions that can occur and their propensity function, the reactivity of the ensemble of molecules can be stored via a global propensity function. If this information can be updated after every reaction event, the evolution of the population of molecules no longer requires the storage of a huge predefined reaction network. This approach is particularly interesting for processes where the initial mixture is limited to a few different molecules because the initial propensity functions are easily computed. Shahrouzi (2010) applied this approach to olefin oligomerization reactions. The main benefit of this Markovian simulation technique resides in the exploration of different parts of a theoretically quasiinfinite reaction network without having to store it.
Full Stochastic Kinetic Modeling of Hydrotreating and Hydroconversion Processes
For refining processes that treat complex mixtures, such as heavy oil fractions, a full stochastic kinetic modeling approach was developed. In the proposed twostep kinetic modeling approach (Fig. 16), the first step, the ‘composition modeling step’, overcomes the lack of molecular detail of complex petroleum fractions by using a set of molecules that are carefully selected based on available analyses. The available analytical information, which is generally relatively limited, is transformed into a set of molecules that has the same properties as the petroleum fraction. This transformation is carried out through the abovedescribed SR – Entropy Maximization (SRREM) algorithm, which stochastically assembles various structural blocks before adjusting the molar fractions of the resulting molecules. Once the composition of the feedstock is adequately represented, the second step, called the ‘reaction modeling step’, simulates the effect of the reactions on this set of molecules by means of a kinetic MC method. For this simulation, the abovedescribed SSA was selected, thereby retaining the molecular detail of the process throughout the entire reaction simulation. This approach is based on a Markovian process, which tracks the transformation of a discrete population of molecules, event by event. Each event describes the transition from one state to another as the result of a chemical reaction. To perform the reactor simulations, the Markovian technique does not need the full reaction network, but only a list of reaction rules that define the possible reactions for each molecule. The probability of each reaction is determined depending on the reactant molecules and the reaction rate constant. By tracking the evolution of the population of molecules, the reaction network is generated ‘onthefly’. Overall, this twostep methodology therefore models both the feedstock composition and the process reactions at a molecular level.
Figure 16. Twostep stochastic kinetic modeling approach. 
This modeling methodology was first developed and validated on somewhat simpler gas oil feeds (de Oliveira et al., 2012; López Abelairas et al., 2016). The hydroprocessing of industrial LCO gas oils was simulated and compared to experiments carried out in an isothermal fixedbed upflow reactor containing 200 mL of a sulfided commercial NiMo/Al_{2}O_{3} catalyst (López Garcia, 2000; López García et al., 2003) and operating at 320 °C and 70 bar. The LCO gas oils were characterized through elemental analysis (carbon, hydrogen, sulfur), MS, ^{13}C RMN and simulated distillation. The LCO gas oils were represented by 5000 molecules which were reconstructed from these analyses by using the probability distribution functions of the structural attributes. To simulate the hydrogenation of the LCO gas oils, a limited set of reaction types needs to be considered: hydrogenation of aromatic rings, dehydrogenation of cyclohexane rings, and hydrodesulfurization of sulfur species. The evolution of the set of 5000 molecules was assessed by performing 50 stochastic simulations and averaging the properties of the mixture.
For the LCO gas oils, the content of the aromatic families of the reconstructed feed mixture shows a good agreement with the analytical values. The same is true for the other properties and illustrates the representativeness of the generated set of 5000 molecules, thus proving that the LCO gas oil is well reconstructed. The simulations of the LCO hydrotreating process showed that the overall properties of the effluent were well predicted for various feedstocks and operating conditions (Fig. 17). The prediction of the hydrogenation of aromatic compounds is quite accurate. The results illustrate that the methodology is able to provide a good prediction of the behavior of hydrotreating process behavior on complex feeds with different molecular compositions. The methodology is also able to predict molecular properties of the effluent that are not accessible in traditional kinetic models.
Figure 17. Comparison between experimental and simulated effluent weight fractions for the hydrocarbon families (SAT – Saturates, MONO – Monoaromatics, DI – Diaromatics, TRI+ – Polyaromatics with at least three aromatic cycles) and sulfur families (BT – Benzothiophenes, DBT – Dibenzothiophenes) during hydrotreating of two different LCO gas oils. 
To illustrate the full potential of the method, the abovedescribed methodology was also applied to the hydroconversion of an Athabasca VR. Elemental analysis (C, H, S, N, O), average molecular weight, specific gravity, simulated distillation, SARA separation and ^{13}C NMR were used to characterize the VR feed analytically (de Oliveira et al., 2013a). The conversion of the VR was again simulated using the SSA. For each molecule, an algorithm makes an inventory of the possible reactions of the molecule, and calculates the corresponding rate coefficients by means of appropriate QSRC for the hydrogenationdehydrogenation equilibrium constants, the hydrogenation and hydrodesulfurization rate constants, the ring dealkylation, and the ring opening rate constants:$$\mathrm{ln}\left({K}_{\mathrm{eq}}\right)=2.951813.2145\xb7{n}_{H}+0.36999\xb7\left\frac{{\u2206H}_{R}^{0}}{{\mathrm{RT}}^{0}}\right0.7837\xb7{N}_{\mathrm{SR}}\frac{\u2206{H}_{R}}{R}\xb7\left(\frac{1}{T}\frac{1}{{T}_{\mathrm{Ref}}}\right)$$(28) $$\mathrm{ln}\left({c}_{\mathrm{Hydro}}\right)=9.85297.0607\xb7{n}_{H}^{1.18}+76.751\xb7{10}^{3}\xb7\left\frac{{\u2206H}_{R}^{0}}{{\mathrm{RT}}^{0}}\right+0.61632\xb7{N}_{\mathrm{AR}}+0.33026\xb7{N}_{\mathrm{SR}}2.8000\xb7{N}_{\mathrm{TR}}\frac{{E}_{a,\mathrm{Hydro}}}{R}\xb7\left(\frac{1}{T}\frac{1}{{T}_{\mathrm{Ref}}}\right)$$(29) $$\mathrm{ln}\left({c}_{\mathrm{HDS}}\right)=2.7500.4699\xb7{N}_{\mathrm{AR}}^{2.291}5.076\xb7{N}_{\mathrm{Alkyl}}^{0.540}\frac{{E}_{a,\mathrm{HDS}}}{R}\xb7\left(\frac{1}{T}\frac{1}{{T}_{\mathrm{Ref}}}\right)$$(30) $$\mathrm{ln}\left({c}_{\mathrm{Dealk}}\right)=0.2331\mathrm{ln}\left(1+\mathrm{exp}(0.50\xb7({N}_{\mathrm{atoms}}8.00)\right)\frac{{E}_{a,\mathrm{Dealk}}}{R}\xb7\left(\frac{1}{T}\frac{1}{{T}_{\mathrm{Ref}}}\right)$$(31) $$\mathrm{ln}\left({c}_{\mathrm{RingOp}}\right)=0.3285\frac{{E}_{a,\mathrm{RingO}}}{R}\xb7\left(\frac{1}{T}\frac{1}{{T}_{\mathrm{Ref}}}\right)$$(32)
In the above equations, the reaction enthalpies were calculated by means of Benson’s group contribution method; while all activation energies were directly taken from literature: the hydrogenation and desulfurization activation energies were taken from López García (2000), and the dealkylation activation energies from taken from Gauthier et al. (2007) and DanialFortain (2010).$${E}_{a,\mathrm{Hydro}}=50900\mathrm{J}\mathrm{mol}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\mathrm{with}{T}_{\mathrm{ref}}=598.15\mathrm{K}$$(33) $${E}_{a,\mathrm{HDS}}=50000+20000\xb7{N}_{\mathrm{AR}}^{0.807}+28350\xb7{N}_{\mathrm{Alkyl}}^{0.321}\hspace{1em}\mathrm{with}{T}_{\mathrm{ref}}=598.15\mathrm{K}$$(34) $${E}_{a,\mathrm{Dealk}}=217000\mathrm{J}/\mathrm{mol}\hspace{1em}\mathrm{with}{T}_{\mathrm{ref}}=688.15\mathrm{K}$$(35) $${E}_{a,\mathrm{RingO}}=50000\mathrm{J}/\mathrm{mol}\hspace{1em}\mathrm{with}{T}_{\mathrm{ref}}=688.15\mathrm{K}$$(36)
The VR was represented by a set of 5000 molecules in order to ensure an optimal balance between the required CPU time and the accuracy of the feedstock representation. The mixture properties obtained after the SR step show a good agreement with most of the analytical data. These results illustrate that the reconstructed feedstock composed of 5000 molecules represents quite accurately the Athabasca VR, although VR fractions contain more complex molecules than LCO gas oils. For the hydroconversion of VR fractions, the evolution of the synthetic mixture of 5000 molecules was obtained by performing 100 stochastic simulations and averaging the properties of the populations of molecules. The overall conversion is well predicted and the overall product yields and product property predictions follow the correct trends with the operating conditions (de Oliveira et al., 2014). The reactor simulations compared favorably (Fig. 18) to experiments carried out in a batch reactor at 410 °C and 150 bar (DanialFortain et al., 2010). The effect of temperature was validated by comparing the simulation results to experimental data on a Ural VR at 395 °C and 410 °C.
Figure 18. Comparison between experimental data and simulation results for (a) residue conversion, (b) sulfur removal, and (c) product yields during hydroconversion of Ural VR at 395 °C and 410 °C. 
In conclusion, this twostep kinetic modeling strategy has been applied to the hydroprocessing of oil fractions: hydrotreating of LCO and hydroconversion of VR. The simulation results were compared to the experimental data, illustrating that the combined ‘molecular reconstruction – stochastic simulation’ approach nicely simulates the process performances. The proposed approach has several advantages. First of all, the complex feedstock is represented by means of a synthetic mixture of molecules. This representation allows the retention of a molecular description of the reaction system throughout the entire reactor simulation. Secondly, the Markovian approach of the SSA allows the tracking of the evolution of the population event by event, without a predefined reaction network. To this end, only a limited number of reaction types need to be defined, and the reaction network is generated ‘onthefly’ as the reactions proceed. Finally, by following the modifications to the synthetic mixture throughout the reactor simulation, the evolution of all product yields and product properties can be assessed by evaluating the yields, composition and characteristics of subgroups of molecules. This means that this methodology is also able to predict various molecular properties of the effluent that are not accessible in traditional lumped kinetic models.
Conclusions
Chemical kinetics studies the rates of chemical reactions and addresses how the reaction rates depend on concentrations, temperature, pH, solvents, nature of the catalysts, to mention a few reaction conditions.
In this review paper, an overview was given of kinetic modeling techniques for complex processes. When systems with several hundred or several thousand components, reaction intermediates, reactions, and elementary steps are being considered, two mean modeling strategies can applied: a lumping strategy or a detailed modeling strategy. In the first strategy, the chemical complexity is reduced by grouping chemical compounds into families or ‘lumps’ and grouping individual reactions into lumped reactions between these families. In this case, a meanfield approach together with a RDS approximation is generally used. Two examples of lumped kinetic models (atmospheric gasoil hydrotreating and residue hydroprocessing) developed at IFPEN have been presented.
The largest part of this review described detailed kinetic modeling approaches. In these detailed modeling strategies, the chemical detail is generally retained at the molecular level, or even at the level of the reaction intermediates throughout the entire reactor simulation, which allows the creation of the feedstockindependent models. This implies, however, that a molecular description of the feedstock can be provided to the model. When advanced analytical techniques are not able to provide a detailed characterization of the feedstock, a molecular representation of the feedstock can be created through molecular reconstruction algorithms, also called composition modeling techniques. Another hurdle to be taken concerns the generation of large size reaction networks. Indeed, as the number of species and the number of reactions increases almost exponentially with increasing carbon number, the reaction network becomes so large that it needs to be computer generated and managed along the simulation. In this review, deterministic and stochastic network generation and simulation approaches have been discussed, together with examples of kinetic models developed at IFPEN for several refining process applications.
All these examples given illustrate that one can acquire a more fundamental understanding of very complex processes by applying advanced techniques for composition, kinetic and reactor modeling.
References
 Adam F., Bertoncini F., Coupard V., Charon N., Thiébaut D., Espinat D., Hennion M.C. (2008) Using Comprehensive TwoDimensional Gas Chromatography for the Analysis of Oxygenates in Middle Distillates I. Determination of the Nature of Biodiesels Blends in Diesel Fuel, J. Chromatogr. A 1186, 236–244. [CrossRef] [PubMed] [Google Scholar]
 Adam F., Bertoncini F., Dartiguelongue C., Marchand K., Thiébaut D., Hennion M.C. (2009) Comprehensive TwoDimensional Gas Chromatography for Basic and Neutral Nitrogen Speciation in Middle Distillates, Fuel 88, 938–946. [CrossRef] [Google Scholar]
 Ahmad M.I., Zhang N., Jobson M. (2011) Molecular ComponentsBased Representation of Petroleum Fractions, Chem. Eng. Res. Des. 89, 410–420. [CrossRef] [Google Scholar]
 Al Halwachi H.K., Yakovlev D.S., Boek E.S. (2012) Systematic Optimization of Asphaltene Molecular Structure and Molecular Weight Using the Quantitative Molecular Representation Approach, Energy & Fuels 26, 10, 6177–6185. [CrossRef] [Google Scholar]
 Alwahabi S.M., Froment G.F. (2004) Single Event Kinetic Modeling of the MethanoltoOlefins Process on SAPO34, Ind. Eng. Chem. Res. 43, 17, 5098–5111. [CrossRef] [Google Scholar]
 AlZaid K., Khan Z.H., Hauser A., AlRabiah H. (1998) Composition of high boiling petroleum distillates of Kuwait crude oils, Fuel 77, 5, 453–458. [Google Scholar]
 Ali F., Ghaloum N., Hauser A. (2006) Structure representation of asphaltene GPC fractions derived from Kuwaiti residual oils, Energy & Fuels 19, 231–238. [CrossRef] [Google Scholar]
 Ali L.H., AlGhannam K.A., AlRawi J.M. (1990) Chemical Structure of Asphaltenes in Heavy Crude Oils Investigated by N.M.R, Fuel 69, 519–521. [CrossRef] [Google Scholar]
 Allen D.T., Liguras D.K. (1991) Structural Models of Catalytic Cracking Chemistry: A Case Study of Group Contribution Approach to Lumped Kinetic Modeling, in Chemical Reactions in Complex Mixtures: The Mobil Workshop, Sapre A.M., Krambeck F.J. (eds), Springer, pp. 101–125. [CrossRef] [Google Scholar]
 Altgelt H.K., Boduszynski M.M. (1994) Composition and Analysis of Heavy Petroleum Fractions, Marcel Dekker, Inc. [Google Scholar]
 AlvarezMajmutov A., Chen J., Gieleciak R., Hager D., Heshka N., Salmon S. (2014) Deriving the Molecular Composition of Middle Distillates by Integrating Statistical Modeling with Advanced Hydrocarbon Characterization, Energy & Fuels 28, 12, 7385–7393. [CrossRef] [Google Scholar]
 AlvarezMajmutov A., Gieleciak R., Chen J. (2015) Deriving the Molecular Composition of Vacuum Distillates by Integrating Statistical Modeling and Detailed Hydrocarbon Characterization, Energy & Fuels 29, 12, 7931–7940. [CrossRef] [Google Scholar]
 AlvarezMajmutov A., Chen J., Gieleciak R. (2016) MolecularLevel Modeling and Simulation of Vacuum Gas Oil Hydrocracking, Energy & Fuels 30, 1, 138–148. [CrossRef] [Google Scholar]
 Andersson J.T., Sielex K. (1996) Dimethylbenzothiophenes and Methyldibenzothiophenes in Crude Oils from Different Sources, J. High Resolut. Chromatogr. 19, 49–53. [CrossRef] [Google Scholar]
 Artok L., Su Y., Hirose Y., Hosokawa M., Murata S., Nomura M. (1999) Structure and Reactivity of PetroleumDerived Asphaltene, Energy & Fuels 13, 287–296. [CrossRef] [Google Scholar]
 ASTM D2425 (2007) Standard Test Method for Hydrocarbon Types in Middle Distillates by Mass Spectrometry, Annu. B. ASTM Stand. [Google Scholar]
 ASTM D2786 (2007) Standard Test Method for Hydrocarbon Types Analysis of Gas Oil Saturates Fractions by High Ionizing Voltage Mass Spectrometry, Annu. B. ASTM Stand. [Google Scholar]
 ASTM D6730 (2007) Standard Test Method for Determination of Individual Components in Spark Ignition Engine Fuels by 100–Metre Capillary (with Precolumn) High Resolution Gas, Annu. B. ASTM Stand. [Google Scholar]
 Aye M.M.S., Zhang N. (2005) A Novel Methodology in Transforming Bulk Properties of Refining Streams into Molecular Information, Chem. Eng. Sci. 60, 6702–6717. [CrossRef] [Google Scholar]
 Baltanas M.A., Froment G.F. (1985) Computer Generation of Reaction Networks and Calculation of Product Distributions in the Hydroisomerization and Hydrocracking of Paraffins on Ptcontaining Bifunctional Catalysts, Comp. Chem. Eng. 9, 1, 71–81. [Google Scholar]
 Baltanas M.A., van Raemdonck K.K., Froment G.F., Mohedas S.R. (1989) Fundamental Kinetic Modeling of Hydroisomerization and Hydrocracking on Noble Metalloaded Faujasites. 1. Rate Parameters for Hydroisomerization, Ind. Eng. Chem. Res. 28, 7, 899–910. [CrossRef] [Google Scholar]
 BattinLeclerc F., Glaude P.A., Warth V., Fournet R., Scacchi G., Côme G.M. (2000) Computer Tools for Modelling the Chemical Phenomena Related to Combustion, Chem. Eng. Sci. 55, 15, 2883–2893. [CrossRef] [Google Scholar]
 Beirnaert H.C., Alleman J.R., Marin G.B. (2001) A Fundamental Kinetic Model for the Catalytic Cracking of Alkanes on a USY Zeolite in the Presence of Coke Formation, Ind. Eng. Chem. Res. 40, 5, 1337–1347 [CrossRef] [Google Scholar]
 Benson S.W. (1976) Thermochemical Kinetics: Methods for the Estimation of Thermochemical Data and Rate Parameters, 2nd ed., Wiley Interscience, ISBN 9780471067818. [Google Scholar]
 Benson S.W., Cohen N. (1993) Estimation of Heats of Formation of Organic Compounds by Additivity Methods, Chemical Reviews 93, 2419–2438. [Google Scholar]
 Bera T., Thybaut J.W., Marin G.B. (2011) SingleEvent MicroKinetics of Aromatics Hydrogenation on Pt/HZSM22, Ind. Eng. Chem. Res. 50, 23, 12933–12945. [CrossRef] [Google Scholar]
 Bera T., Thybaut J.W., Marin G.B. (2012) Extension of the SingleEvent Microkinetic Model to Alkyl Substituted Monoaromatics Hydrogenation on a Pt Catalyst, ACS Catalysis 2, 7, 1305–1318. [CrossRef] [Google Scholar]
 Bertoncini F., Courtiade M., Thiébaut D. (2013) Gas Chromatography and 2DGas Chromatography for Petroleum Industry, Editions Technip, Paris (France). [Google Scholar]
 Bhan A., Hsu S.H., Blau G., Caruthers J.M., Venkatasubramanian V., Delgass W.N. (2005) Microkinetic Modeling of Propane Aromatization over HZSM5, Journal of Catalysis 235, 35–51. [CrossRef] [Google Scholar]
 Billaud F., Elyahyaoui K., Baronnet F., Kressmann S. (1991) Chemical kinetic modeling of nhexane pyrolysis of ACUCHEM, CHEMKIN and MORSE software packages, Chem. Eng. Sci. 46, 11, 2941–2946. [CrossRef] [Google Scholar]
 Blinov M.L., Yang J., Faeder J.R., Hlavacek W.S. (2005a) Graph Theory for Rulebased Modeling of Biochemical Networks, Transactions on Computational Systems Biology VII: Lecture Notes in Computer Science 4230, 89–106. [CrossRef] [Google Scholar]
 Blinov M.L., Faeder J.R., Yang J., Goldstein B., Hlavacek W.S. (2005b) ‘Onthefly’ or ‘generatefirst’ modeling?, Nature Biotechnology 23, 11, 1344–1345 [CrossRef] [Google Scholar]
 Blurock E.S. (1995) Reaction: System for Modeling Chemical Reactions, J. Chem. Inf. Comput. Sci. 35, 3, 607–616. [CrossRef] [Google Scholar]
 Boduszynski M.M. (1987) Composition of Heavy Petroleums. 1. Molecular Weight, Hydrogen Deficiency, and Heteroatom Concentration as a Function of Atmospheric Equivalent Boiling Point up to 1400.degree.F (760.degree.C), Energy & Fuels 1, 2–11. [CrossRef] [Google Scholar]
 Boduszynski M.M. (1988) Composition of Heavy Petroleums. 2. Molecular Characterization, Energy & Fuels 2, 597–613. [CrossRef] [Google Scholar]
 Boek E.S., Yakovlev D.S., Headen T.F. (2009) Quantitative Molecular Representation of Asphaltenes and Molecular Dynamics Simulation of Their Aggregation, Energy & Fuels 23, 3, 1209–1219. [CrossRef] [Google Scholar]
 Bonnardot J. (1998) Modèlisation cinétique des réactions d’hydrotraitement par regroupement en familles chimiques, Ph.D. Thesis, Université Claude Bernard – Lyon (France). [Google Scholar]
 Bounaceur R., Warth V., Marquaire P.M., Scacchi G., Domine F., Dessort D., Pradier B., Brevart O. (2002) Modeling of Hydrocarbons Pyrolysis at Low Temperature. Automatic Generation of Free Radicals Mechanisms, J. Anal. Appl. Pyrol. 64, 1, 103–122. [CrossRef] [Google Scholar]
 Bouyssiere B., Leonhard P., Profrock D., Baco F., López García C., Wilbur S., Prange A. (2004) Investigation of the sulfur speciation in petroleum products by capillary gas chromatography with ICPcollision cellMS detection, Journal of Analytical Atomic Spectroscopy 19, 5, 700–702. [CrossRef] [Google Scholar]
 Broadbelt L.J., Stark S.M., Klein M.T. (1994) ComputerGenerated Pyrolysis Modeling – OntheFly Generation of Species, Reactions, and Rates, Ind. Eng. Chem. Res. 33, 4, 790–799. [CrossRef] [Google Scholar]
 Broadbelt L.J., Klein M.T., Dean B.D., Andrews S.M. (1995) Thermal Degradation of AliphaticAromatic Polyamides: Kinetics of N,N′Dihexylisophthalamide Neat and in Presence of Copper Iodide, Journal of Applied Polymer Science 56, 7, 803–815. [CrossRef] [Google Scholar]
 Brown J.K., Ladner W.R. (1960) A Study on the Hydrogen Distribution in Coallike Materials by HighResolution Nuclear Magnetic Resonance Spectroscopy II: A Comparison with InfraRed Measurement and the Conversion to Carbon Structure, Fuel 39, 87–96. [Google Scholar]
 Brønsted J.N. (1928) Acid and Basic Catalysis, Chem. Rev. 5, 3, 231–338. [CrossRef] [Google Scholar]
 Brønsted J.N., Pedersen K.J. (1924) Die katalytische Zersetzung des Nitramids und ihre physikalischchemische Bedeutung, Zeitschrift für Phys. Chemie 108, 185–235. [Google Scholar]
 Bruk L.G., Gorodskii S.N., Zeigarnik A.V., ValdésPérez R.E., Temkin O.N. (1998) Oxidative Carbonylation of Phenylacetylene Catalyzed by Pd(II) and Cu(I): Experimental Tests of FortyOne ComputerGenerated Mechanistic Hypotheses, J. Mol. Catal. A: Chem. 130, 1–2, 29–40. [CrossRef] [Google Scholar]
 Buda F., Heyberger B., Fournet R., Glaude P.A., Warth V., BattinLeclerc F. (2006) Modeling of the GasPhase Oxidation of Cyclohexane, Energy & Fuels 20, 4, 1450–1459. [CrossRef] [Google Scholar]
 Bytautas L., Klein D.J. (1998) Chemical Combinatorics for AlkaneIsomer Enumeration and more, J. Chem. Inf. Comput. Sci. 38, 1063–1078. [CrossRef] [Google Scholar]
 Campbell D.M., Klein M.T. (1997) Construction of a Molecular Representation of a Complex Feedstock by Monte Carlo and Quadrature Methods, Appl. Catal. A Gen. 160, 41–54. [CrossRef] [Google Scholar]
 Cantor D.M. (1978) Nuclear magnetic resonance spectrometric determination of average molecular structure parameters for coalderived liquids, Analytical Chemistry 50, 8, 1185–1187. [CrossRef] [Google Scholar]
 Cayley A. (1875) Über die analytischen Figuren, welche in der Mathematik Baeume genannt warden, Chem. Ber. 8, 1056–1059. [CrossRef] [Google Scholar]
 Chapman N.B., Shorter J. (eds) (1978) Correlation analyses in chemistry: Recent advances, Plenum, New York. [CrossRef] [Google Scholar]
 CharonRevellin N., Dulot H., López García C., Jose J. (2011) Kinetic Modeling of Vacuum Gas Oil Hydrotreatment using a Molecular Reconstruction Approach, Oil Gas Sci. Technol. – Rev. d’IFP Energies Nouvelles 66, 3, 479–490. [CrossRef] [EDP Sciences] [Google Scholar]
 ChavarríaHernández J.C., Ramírez J., Gonzalez H., Baltanas M.A. (2004) Modelling of nHexadecane Hydroisomerization and Hydrocracking Reactions on a Mo/H BetaAlumina BiFunctional Catalyst Using the Single Event Concept, Catal. Today 98, 1–2, 235–242. [CrossRef] [Google Scholar]
 ChavarríaHernández J.C., Ramírez J., Baltanas M.A. (2008) SingleEvent LumpedParameter Hybrid (SELPH) Model for NonIdeal Hydrocracking of nOctane, Catal. Today 130, 2–4, 455–461. [CrossRef] [Google Scholar]
 ChavarríaHernández J.C., Ramírez J. (2009) Modeling Ideal and Nonideal Hydrocracking of Paraffins Using the SingleEvent LumpedParameter Hybrid (SELPH) Model, Ind. Eng. Chem. Res. 48, 3, 1203–1207. [CrossRef] [Google Scholar]
 Chevalier C., Warnatz J., Melenk H. (1990) Automatic Generation of Reaction Mechanisms for Description of Oxidation of Higher Hydrocarbons, Ber. Bunsenges. Phys. Chem. 94, 1362–1367. [Google Scholar]
 Chevalier C., Pitz W.J., Warnatz J., Westbrook C.K., Melenk H. (1992) Hydrocarbon Ignition: Automatic Generation of Reaction Mechanisms and Applications to Modeling of Engine Shock, Proc. Combust. Inst. 24, 1362–1367. [Google Scholar]
 Chinnick S.J., Baulch D.L., Ayscough P.B. (1988) An Expert System for Hydrocarbon Pyrolysis Reactions, Chemometrics and Intelligent Laboratory Systems 5, 1, 39–52. [CrossRef] [Google Scholar]
 Choudhury I.R., Thybaut J.W., Balasubramanian P., Denayer J.F.M., Martens J.A., Marin G.B. (2010) Synergy between Shape Selective and NonShape Selective Bifunctional Zeolites Modelled via the SingleEvent MicroKinetic (SEMK) Methodology, Chem. Eng. Sci. 65, 1, 174–178. [CrossRef] [Google Scholar]
 Christensen G., Apelian M.R., Hickey K.J., Jaffe S.B. (1999) Future Directions in Modeling the FCC Process: An Emphasis on Product Quality, Chem. Eng. Sci. 54, 2753–2764. [CrossRef] [Google Scholar]
 Clymans P.J., Froment G.F. (1984) Computer Generation of Reaction Paths and Rate Equations in the Thermal Cracking of Normal and Branched Paraffins, Comp. Chem. Eng. 8, 2, 137–142. [CrossRef] [Google Scholar]
 Cochegrue H., Gauthier P., Verstraete J.J., Surla K., Guillaume D., Galtier P., Barbier J. (2011) Reduction of Single Event Kinetic Models by Rigorous Relumping: Application to Catalytic Reforming, Oil Gas Sci. Technol. – Rev. d’IFP Energies Nouvelles 66, 3, 367–397. [CrossRef] [EDP Sciences] [Google Scholar]
 Côme G.M., Warth V., Glaude P.A., Fournet R., BattinLeclerc F., Scacchi G. (1996) ComputerAided Design of GasPhase Oxidation Mechanisms. Application to Modeling of nHeptane and IsoOctane Oxidation, 26^{th} Int. Symposium on Combustion, Symposium (International) on Combustion 26, 1, 755–762. [CrossRef] [Google Scholar]
 Connors K.A. (1990) Chemical Kinetics: The Study of Reaction Rates in Solution, VCH Verlagsgesellschaft Weinheim, New York (NY). [Google Scholar]
 Constantinou L., Gani R. (1994) New Group Contribution Method for Estimating Properties of Pure Compounds, AIChE J. 40, 1697–1710. [Google Scholar]
 DanialFortain P. (2010) Étude de la réactivité des résidus pétroliers en hydroconversion, Ph.D. Thesis, Université Bordeaux 1, Bordeaux (France). [Google Scholar]
 DanialFortain P., Gauthier T., Merdrignac I., Budzinski H. (2010) Reactivity Study of Athabasca Vacuum Residue in Hydroconversion Conditions, Catalysis Today 150, 255–263. [CrossRef] [Google Scholar]
 Daubert T.E., Danner R.P. (1989) Physical and Thermodynamic Properties of Pure Compounds: Data Compilation, Hemisphere, New York. [Google Scholar]
 De Bruycker R., Pyl S.P., Reyniers M.F., Van Geem K.M., Marin G.B. (2015) Microkinetic model for the pyrolysis of methyl esters: From model compound to industrial biodiesel, AIChE Journal 61, 12, 4309–4322. [CrossRef] [Google Scholar]
 de Oliveira L.P. (2013) Développement d’une méthodologie de modélisation cinétique de procédés de raffinage traitant les charges lourdes, Ph.D. Thesis, École Normale Supérieure de Lyon (France). [Google Scholar]
 de Oliveira L.P., Verstraete J.J., Kolb M. (2012) A Monte Carlo Modeling Methodology for the Simulation of Hydrotreating Processes, Chem. Eng. J. 207–208, 94–102. [CrossRef] [Google Scholar]
 de Oliveira L.P., Trujillo Vazquez A., Verstraete J.J., Kolb M. (2013a) Molecular Reconstruction of Petroleum Fractions: Application to Various Vacuum Residues, Energy & Fuels 27, 3622–3641. [CrossRef] [Google Scholar]
 de Oliveira L.P., Verstraete J.J., Kolb M. (2013b) Moleculebased Kinetic Modeling by Monte Carlo Methods for Heavy Petroleum Conversion, Science China Chem. 16, 11, 1608–1622. [CrossRef] [Google Scholar]
 de Oliveira L.P., Verstraete J.J., Kolb M. (2013c) Development of a General Modelling Methodology for Vacuum Residue Hydroconversion, Oil & Gas Science and Technology – Revue IFP Energies nouvelles 68, 6, 1027–1038. [CrossRef] [EDP Sciences] [Google Scholar]
 de Oliveira L.P., Verstraete J.J., Kolb M. (2014) Simulating Vacuum Residue Hydroconversion by means of MonteCarlo Techniques, Catal. Today 220–222, 208–220. [Google Scholar]
 Dente M., Pierucci S., Ranzi E., Bussani G. (1992) New Improvements in Modeling Kinetic Schemes for Hydrocarbons Pyrolysis Reactors, Chem. Eng. Sci. 51, 2629–2634. [CrossRef] [Google Scholar]
 Depauw G.A., Froment G.F. (1997) Molecular Analysis of the Sulphur Components in a Light Cycle Oil of a Catalytic Cracking Unit by Gas Chromatography with Mass Spectrometric and Atomic Emission Detection, J. Chromatogr. A 761, 231–247. [CrossRef] [Google Scholar]
 Dewachtere N.V., Santaella F., Froment G.F. (1999) Application of a SingleEvent Kinetic Model in the Simulation of an Industrial Riser Reactor for the Catalytic Cracking of Vacuum Gas Oil, Chem. Eng. Sci. 54, 15–16, 3653–3660. [CrossRef] [Google Scholar]
 DeWitt M.J., Dooling D.J., Broadbelt L.J. (2000) Computer Generation of Reaction Mechanisms Using Quantitative Rate Information: Application to LongChain Hydrocarbon Pyrolysis, Ind. Eng. Chem. Res. 39, 7, 2228–2237. [CrossRef] [Google Scholar]
 Dickinson E.M. (1980) Structural comparison of petroleum fractions using proton and ^{13}C NMR spectroscopy, Fuel 59, 290–294. [CrossRef] [Google Scholar]
 Di Maio F.P., Lignola P.G. (1992) KING, a Kinetic Network Generator, Chem. Eng. Sci. 47, 9–11, 2713–2718. [CrossRef] [Google Scholar]
 Dutriez T., Courtiade M., Thiebaut D., Dulot H., Bertoncini F., Vial J., Hennion M.C. (2009) HighTemperature TwoDimensional Gas Chromatography of Hydrocarbons up to nC(60) for Analysis of Vacuum Gas Oils, Journal of Chromatography A 1216, 14, 2905–2912. [CrossRef] [PubMed] [Google Scholar]
 Dutriez T., Courtiade M., Thiebaut D., Dulot H., Hennion M.C. (2010) Improved Hydrocarbons Analysis of Heavy Petroleum Fractions by High Temperature Comprehensive TwoDimensional Gas Chromatography, Fuel 89, 9, 2338–2345. [CrossRef] [Google Scholar]
 Dutriez T., Borras J., Courtiade M., Thiebaut D., Dulot H., Bertoncini F., Hennion M.C. (2011) Challenge in the speciation of nitrogencontaining compounds in heavy petroleum fractions by high temperature comprehensive twodimensional gas chromatography, Journal of Chromatography A 1218, 21, 3190–3199. [CrossRef] [PubMed] [Google Scholar]
 Dzidic I., Balicki M.D., Rhodes I.A.L., Hart H.V. (1988) Identification and Quantification of Nitrogen and Sulfur Compounds in Catalytically Cracked Heavy Gas Oils by Isobutane/CI GC/MS and GC Using Selective Detectors, J. Chromatogr. Sci. 26, 236–240. [CrossRef] [Google Scholar]
 Evans M.G., Polanyi M. (1935) Some applications of the transition state method to the calculation of reaction velocities, especially in solution, Trans. Faraday Soc. 31, 875–894. [CrossRef] [Google Scholar]
 Evans M.G., Polanyi M. (1936) Further considerations on the thermodynamics of chemical equilibria and reaction rates, Trans. Faraday Soc. 32, 1333–1360. [Google Scholar]
 Evans M.G., Polanyi M. (1938) Inertia and driving force of chemical reactions, Trans. Faraday Soc. 34, 11–24. [Google Scholar]
 Faeder J.R., Hlavecek W.S., Reischl I., Blinov M.L., Metzger H., Redondo A., Wofsy C., Goldstein B. (2003) Investigation of early events in Fc epsilon RImediated signaling using a detailed mathematical model, Journal of Immunology 170, 7, 3769–3781. [CrossRef] [Google Scholar]
 Faeder J.R., Blinov M.L., Goldstein B., Hlavacek W.S. (2005) Rulebased Modeling of Biochemical Networks, Complexity 10, 22–41. [CrossRef] [Google Scholar]
 Fafet A., Bonnard J., Prigent F. (1999a) New Developments in Mass Spectrometry for GroupType Analysis of Petroleum Cuts (First Part), Oil Gas Sci. Technol.  Rev. l’IFP Energies Nouvelles 54, 439–452. [Google Scholar]
 Fafet A., Bonnard J., Prigent F. (1999b) New Developments in Mass Spectrometry for GroupType Analysis of Petroleum Cuts (Second Part), Oil Gas Sci. Technol.  Rev. l’IFP Energies Nouvelles 54, 453–462. [CrossRef] [EDP Sciences] [Google Scholar]
 Faulon J.L. (1994) Stochastic Generator of Chemical Structure. 1. Application to the Structure Elucidation of Large Molecules, J. Chem. Inf. Model. 34, 1204–1218. [CrossRef] [Google Scholar]
 Faulon J.L., Vandenbroucke M., Drappier J.M., Behar F., Romero M. (1990) Modélisation des structures chimiques des macromolécules sédimentaires: le logiciel XMOL, Oil & Gas Science and Technology – Rev. IFP Energies nouvelles 45, 2, 161–180. [Google Scholar]
 Faulon J.L., Sault A.G. (2001) Stochastic Generator of Chemical Structure. 3. Reaction Network Generation, J. Chem. Inf. Model. 41, 894–908. [Google Scholar]
 Faulon J.L., Visco D.P., Roe D. (2005) Enumerating Molecules, in Reviews in Computational Chemistry, Vol. 21, K.B. Lipkowitz, R. Larter, T.R. Cundari, D.B. Boyd (eds), J. Wiley & Sons, pp. 209–275. [Google Scholar]
 Feng W., Vynckier E., Froment G.F. (1993) SingleEvent Kinetics of Catalytic Cracking, Ind. Eng. Chem. Res. 32, 12, 2997–3005. [CrossRef] [Google Scholar]
 Ferreira C. (2009) Modélisation de l’hydrotraitement de résidus pétroliers en lit fixe. Étude de la réactivité de charges, Ph.D. Thesis, École Normale Supérieure de Lyon (France). [Google Scholar]
 Ferreira C., Marques J., TayakoutFayolle M., Guibard I., Lemos F., Toulhoat H., Ramôa Ribeiro F. (2010) Modeling Residue Hydrotreating, Chem. Eng. Sci. 65, 322–329. [CrossRef] [Google Scholar]
 Ferreira C., TayakoutFayolle M., Guibard I., Lemos F. (2014) Hydrodesulfurization and Hydrodemetallization of Different Origin Vacuum Residues: New Modeling Approach, Fuel 129, 267–277. [CrossRef] [Google Scholar]
 Fierro V., Duplan J.L., Verstraete J.J., Schuurman Y., Mirodatos C. (2001) A nonstationary kinetics approach for the determination of the kinetic parameters of the protolytic cracking of methylcyclohexane, in Reaction kinetics and the development of catalytic processes, G.F. Froment, K.C. Waugh (eds), Studies in Surface Science and Catalysis 133, 341–348. [Google Scholar]
 Fierro V., Schuurman Y., Mirodatos C., Duplan J.L., Verstraete J.J. (2002) Study of the cracking reaction of linear and branched hexanes under protolytic conditions by nonstationary kinetics, Chemical Engineering Journal 90, 1–2, 139–147. [CrossRef] [Google Scholar]
 Finley S.D., Broadbelt L.J., Hatzimanikatis V. (2009) Computational Framework for Predictive Biodegradation, Biotechnology and Bioengineering 104, 6, 1086–1097. [CrossRef] [PubMed] [Google Scholar]
 Fischer I., Fischer P. (1974) Analysis of HighBoiling Petroleum Streams by HighResolution Mass Spectrometry, Talanta 21, 867–875. [CrossRef] [PubMed] [Google Scholar]
 Flory P.J. (1936) Molecular Size Distribution in Linear Condensation Polymers, J. Am. Chem. Soc. 58, 1877–1885. [CrossRef] [Google Scholar]
 Fontain E., Reitsam K. (1991) The Generation of Reaction Networks with RAIN. 1. The Reaction Generator, J. Chem. Inf. Comput. Sci. 31, 1, 96–101. [CrossRef] [MathSciNet] [Google Scholar]
 Frenkel M., Hong X., Wilhoit R. (eds) (2000) TRC Thermodynamic Tables – Hydrocarbons, Washington DC. [Google Scholar]
 Frenklach M. (1987) Modeling of Large Reaction Systems, Complex Chemical Reaction Systems, Mathematical Modelling and Simulation, Springer Series in Chemical Physics, 2–16. [Google Scholar]
 Froment G.F. (1991) Fundamental Kinetic Modeling of Complex Processes, in: Chemical Reactions in Complex Systems: the Mobil Workshop, A.V. Sapre, F.J. Krambeck (eds), Van Nostrand Reinhold, New York, pp. 77–100. [CrossRef] [Google Scholar]
 Froment G.F. (1999) Kinetic Modeling of AcidCatalyzed Oil Refining Processes, Catalysis Today 52, 153–163. [CrossRef] [Google Scholar]
 Froment G.F. (2005) Single Event Kinetic Modeling of Complex Catalytic Processes, Catalysis Reviews 47, 83–124. [CrossRef] [Google Scholar]
 Froment G.F., Depauw G.A., Vanrysselberghe V. (1994) Kinetic Modeling and Reactor Simulation in Hydrodesulfurization of Oil Fractions, Industrial & Engineering Chemistry Research 33, 12, 2975–2988. [CrossRef] [Google Scholar]
 Gaffuri P., Faravelli T., Ranzi E., Cernansky N.P., Miller D., d’Anna A., Ciajolo A. (1997) Comprehensive Kinetic Model for the LowTemperature Oxidation of Hydrocarbons, AIChE J. 43, 5, 1278–1286. [CrossRef] [Google Scholar]
 Gauthier T., Heraud J.P., Kressmann S., Verstraete J.J. (2007) Impact of vaporization in a residue hydroconversion process, Chemical Engineering Science 62, 18–20, 5409–5417. [CrossRef] [Google Scholar]
 Gauthier T., DanialFortain P., Merdrignac I., Guibard I., Quoineaud A.A. (2008) Studies on the evolution of asphaltene structure during hydroconversion of petroleum residues, Catalysis Today 130, 429–438. [CrossRef] [Google Scholar]
 Giddings J.C. (1987) Concepts and Comparisons in Multidimensional Separation, J. High Resolut. Chromatogr. 10, 319–323. [Google Scholar]
 Gillespie D.T. (1976) A general method for numerically simulating the stochastic time evolution of coupled chemical reactions, J. Comput. Phys. 22, 403–434. [Google Scholar]
 Glaude P.A., Warth V., Fournet R., BattinLeclerc F., Côme G.M., Scacchi G. (1997) Modelling of nHeptane and IsoOctane Gas Phase Oxidation at Low Temperature by using ComputerAided Designed Mechanisms, Bull. Soc. Chim. Belg. 106, 6, 343–348. [Google Scholar]
 Glaude P.A., Warth V., Fournet R., BattinLeclerc F., Côme G.M., Scacchi G. (1998) Modelling of the Oxidation of nOctane and nDecane using an Automatic Generation of Mechanisms, International Journal of Chemical Kinetics 30, 949–959. [CrossRef] [Google Scholar]
 Green W.H., Barton P.I., Bhattacharjee B., Matheu D.M., Schwer D.A., Song J., Sumathi R., Carstensen H.H., Dean A.M., Grenda J.M. (2001) Computer Construction of Detailed Chemical Kinetic Models for GasPhase Reactors, Ind. Eng. Chem. Res. 40, 23, 5362–5370. [CrossRef] [Google Scholar]
 Grenda J.M., Androulakis I.P., Dean A.M., Green W.H. (2003) Application of Computational Kinetic Mechanism Generation to Model the Autocatalytic Pyrolysis of Methane, Ind. Eng. Chem. Res. 42, 1000–1010. [CrossRef] [Google Scholar]
 Guillaume D. (2006) Network Generation of Oligomerization Reactions: Principles, Industrial & Engineering Chemistry Research 45, 13, 4554–4557. [CrossRef] [Google Scholar]
 Guillaume D., Surla K., Galtier P. (2003a) From Single Events Theory to Molecular Kinetics  Application to Industrial Process Modelling, Chemical Engineering Science 58, 21, 4861–4869. [CrossRef] [Google Scholar]
 Guillaume D., Surla K., Galtier P. (2003b) Single Events Modelling. Part I: Review, Communication at ECCE4  4th European Congress of Chemical Engineering, Granada (Spain), Sept. 21–25. [Google Scholar]
 Guillaume D., Valéry E., Surla K., Galtier P., Verstraete J., Schweich D. (2003c) Single Events Modelling. Part II: Extension to Large Networks, Communication at ECCE4  4th European Congress of Chemical Engineering, Granada (Spain), Sept. 21–25. [Google Scholar]
 Guillaume D., Valéry E., Verstraete J.J., Surla K., Galtier P., Schweich D. (2011) Single Events Kinetic Modeling without Explicit Generation of Large Networks: Application to Hydrocracking of Long Paraffins, Oil & Gas Science and Technology – Rev. IFP Energies nouvelles 66, 3, 399–422. [CrossRef] [EDP Sciences] [Google Scholar]
 Hammett L.P. (1937) The Effect of Structure upon the Reactions of Organic Compounds. Benzene Derivatives, J. Am. Chem. Soc. 59, 1, 96–103. [CrossRef] [Google Scholar]
 Hatzimanikatis V., Li C.H., Ionita J.A., Broadbelt L.J. (2004) Metabolic Networks: Enzyme Function and Metabolite Structure, Current Opinion in Structural Biology 14, 300–306. [CrossRef] [PubMed] [Google Scholar]
 Hatzimanikatis V., Li C.H., Ionita J.A., Henry C.S., Jankowski M.D., Broadbelt L.J. (2005) Exploring the Diversity of Complex Metabolic Networks, Bioinformatics 21, 1603–1609. [CrossRef] [PubMed] [Google Scholar]
 Haulle F.X. (2002) Modélisation cinétique de l’hydrotraitement en lit fixe des résidus pétrolièrs : Étude de la réactivité des composés soufrés, Ph.D. Thesis, Université Paris VI, Paris (France). [Google Scholar]
 Hendrickson J.B., Parks C.A. (1991) Generation and Enumeration of Carbon skeletons, J. Chem. Inf. Comput. Sci. 31, 101–107. [CrossRef] [Google Scholar]
 Henry C.S., Broadbelt L.J., Hatzimanikatis V. (2007) Thermodynamicsbased Metabolic Flux Analysis, Biophysical Journal 92, 5, 1792–1805. [CrossRef] [PubMed] [Google Scholar]
 Henry C.S., Broadbelt L.J., Hatzimanikatis V. (2010) Discovery and Analysis of Novel Metabolic Pathways for the Biosynthesis of Industrial Chemicals: 3Hydroxypropanoate, Biotechnology and Bioengineering 106, 3, 462–473. [PubMed] [Google Scholar]
 Heyberger B., BattinLeclerc F., Warth V., Fournet R., Come G.M., Scacchi G. (2001) Comprehensive Mechanism for the GasPhase Oxidation of Propene, Combustion and Flame 126, 1780–1802. [CrossRef] [Google Scholar]
 Hill A., Tomshine J., Wedding E., Sotirpolous V., Kaznessis Y. (2008) SynBioSS: The Synthetic Biology Modeling Suite, Bioinformatics 24, 51, 2551–2553. [CrossRef] [PubMed] [Google Scholar]
 Hillewaert L.P., Dierickx J.L., Froment G.F. (1988) Computer Generation of Reaction Schemes and Rate Equations for Thermal Cracking, AIChE J. 34, 1, 17–24. [CrossRef] [Google Scholar]
 Hindmarsch A.C. (1980) LSODE and LSODI, Two New Initial Value Ordinary Differential Equation Solvers, A.C.M. Signum Newsletter 15, 4, 19–21. [Google Scholar]
 Hindmarsch AC. (1983) ODEPACK, a Systematized Collection of ODE Solvers, in Scientific Computing, R.S. Stepleman et al. (eds), IMACS, NorthHolland, Amsterdam, pp. 55–64. [Google Scholar]
 Hirsch E., Altgelt K.H. (1970) Integrated Structural Analysis. Method for the Determination of Average Structural Parameters of Petroleum Heavy Ends, Anal. Chem. 42, 1330–1339. [CrossRef] [Google Scholar]
 Hognon C., Simon Y., Marquaire P.M. (2012) Hydrogen Production by Homogeneous Partial Oxidation of Propane, Energy & Fuels 26, 3, 1496–1508. [CrossRef] [Google Scholar]
 Hou G., Klein M.T. (1999) Molecular Modeling of Gas Oil Hydrodesulfurization, Abstracts of Papers of the American Chemical Society 218, U610–U611. [Google Scholar]
 Hudebine D. (2003) Reconstruction moléculaire de coupes pétrolières, Ph.D. Thesis, École Normale Supérieure de Lyon (France). [Google Scholar]
 Hudebine D., Verstraete J.J. (2004) Molecular Reconstruction of LCO Gas Oils from Overall Petroleum Analyses, Chem. Eng. Sci. 59, 22–23, 4755–4763. [CrossRef] [Google Scholar]
 Hudebine D., Verstraete J.J. (2011) Reconstruction of Petroleum Feedstocks by Entropy Maximization. Application to FCC Gasolines, Oil Gas Sci. Technol. – Rev. d’IFP Energies nouvelles 66, 3, 437–460. [CrossRef] [EDP Sciences] [Google Scholar]
 Hudebine D., Verstraete J.J., Chapus T. (2011) Statistical Reconstruction of Gas Oil Cuts, Oil Gas Sci. Technol. – Rev. d’IFP Energies nouvelles 66, 3, 461–477. [CrossRef] [EDP Sciences] [Google Scholar]
 Ihlenfeldt W.D., Gasteiger J. (1996) ComputerAssisted Planning of Organic Syntheses: The Second Generation of Programs, Angew. Chem. (International Edition in English) 34, 23–24, 2613–2633. [CrossRef] [Google Scholar]
 Iyer S.D., Joshi P.V., Klein M.T. (1998) Automated Model Building and Modeling of Alcohol Oxidation in High Temperature Water, Environmental Progress 17, 4, 221–233. [CrossRef] [Google Scholar]
 Jacob S.M., Gross B., Voltz S.E., Weekman V.W. (1976) A Lumping and Reaction Scheme for Catalytic Cracking, AIChE J. 22, 701–713. [CrossRef] [Google Scholar]
 Jaffe S.B., Freund H., Olmstead W.N. (2005) Extension of StructureOriented Lumping to Vacuum Residua, Ind. Eng. Chem. Res. 44, 9840–9852. [CrossRef] [Google Scholar]
 Jin F., Froment G.F. (2013) Automatic Generation of Reaction Networks for Complex Processes, Computers and Applied Chemistry 30, 1, 1–7. [Google Scholar]
 Joback K.G., Reid R.C. (1987) Estimation of PureComponent Properties from GroupContributions, Chem. Eng. Commun. 57, 233–243. [CrossRef] [Google Scholar]
 Johansen N.G., Ettre L.S., Miller R.L. (1983) Quantitative Analysis of Hydrocarbons by Structural Group Type in Gasolines and Distillates. I Gas Chromatography, J. Chromatogr. A 256, 393–417. [CrossRef] [Google Scholar]
 Joshi P.V. (1998) Molecular and Mechanistic Modeling of Complex Process Chemistries, Ph.D. Thesis, University of Delaware, Delaware (USA). [Google Scholar]
 Joshi P.V., Freund H., Klein M.T. (1999) Directed Kinetic Model Building: Seeding as a Model Reduction Tool, Energy & Fuels 13, 877–880. [CrossRef] [Google Scholar]
 Joshi P.V., Iyer S.D., Klein M.T. (1998) Computer assisted modeling of gas oil Fluid Catalytic Cracking (FCC), Preprints  American Chemical Society. Division of Petroleum Chemistry 42, 3, 537–680, 654–656. [Google Scholar]
 Karaba A., Zamostny P., Lederer J., Belohlav Z. (2013) Generalized Model of Hydrocarbons Pyrolysis Using Automated Reactions Network Generation, Ind. Eng. Chem. Res. 52, 44, 15407–15416. [CrossRef] [Google Scholar]
 Katare S., Caruthers J.M., Delgass W.N., Venkatasubramanian V. (2004) An Intelligent System for Reaction Kinetic Modeling and Catalyst Design, Ind. Eng. Chem. Res. 43, 3484–3512. [CrossRef] [Google Scholar]
 Klein M.T., Hou G., Quann R.J., Wei W., Liao K.H., Yang R.S.H., Campain J.A., Mazurek M.A., Broadbelt L.J. (2002) BioMol: ComputerAssisted Biological Modeling of Complex Chemical Mixtures and Biological Processes at the Molecular Level, Environ. Health Persp. 110, 6, 1025–1029. [CrossRef] [Google Scholar]
 Klein M.T., Hou G., Bertolacini R.J., Broadbelt L.J., Kumar A. (2006) Molecular Modeling in Heavy Hydrocarbon Conversions, CRC Press, Taylor & Francis Group, Boca Raton, FL (USA). [Google Scholar]
 Klinke D.J., Broadbelt L.J. (1999) Construction of a Mechanistic Model of FischerTropsch Synthesis on Ni(111) and Co(0001) Surfaces, Chem. Eng. Sci. 54, 15–16, 3379–3389. [CrossRef] [Google Scholar]
 Knight S.A. (1967) Analysis of Aromatic Petroleum Fractions by Means of Absorption Mode Carbon13 NMR Spectroscopy, Chem. Ind. (November), 1920–1923. [Google Scholar]
 Korre S.C. (1995) Quantitative structure/Reactivity correlations as a reaction engineering tool: Applications to hydrocracking of polynuclear aromatics, Ph.D. Thesis, University of Delaware, 1994. [Google Scholar]
 Kowalewski I., Vandenbroucke M., Huc A.Y., Taylor M.J., Faulon J.L. (1996) Preliminary results on molecular modeling of asphaltenes using structure elucidation programs in conjunction with molecular simulation programs, Energy & Fuels 24, 10, 97–107. [CrossRef] [Google Scholar]
 Kruse T.M., Woo O.S., Broadbelt L.J. (2001) Detailed Mechanistic Modeling of Polymer Degradation: Application to Polystyrene, Chemical Engineering Science 56, 3, 971–979. [CrossRef] [Google Scholar]
 Kruse T.M., Woo O.S., Wong H.W., Kahn S.S., Broadbelt L.J. (2002) Mechanistic Modeling of Polymer Degradation: A Comprehensive Study of Polystyrene, Macromolecules 35, 20, 7830–7844. [CrossRef] [Google Scholar]
 Kruse T.M., Wong H.W., Broadbelt L.J. (2003) Mechanistic Modeling of Polymer Pyrolysis: Polypropylene, Macromolecules 36, 25, 9594–9607. [CrossRef] [Google Scholar]
 Kumar H., Froment G.F. (2007a) A Generalized Mechanistic Kinetic Model for the Hydroisomerization and Hydrocracking of LongChain Paraffins, Ind. Eng. Chem. Res. 46, 12, 4075–4090. [CrossRef] [Google Scholar]
 Kumar H., Froment G.F. (2007b) Mechanistic Kinetic Modeling of the Hydrocracking of Complex Feedstocks, such as Vacuum Gas Oils, Ind. Eng. Chem. Res. 46, 18, 5881–5897. [CrossRef] [Google Scholar]
 Kumar P., Thybaut J.W., Svelle S., Olsbye U., Marin G.B. (2013a) SingleEvent Microkinetics for Methanol to Olefins on HZSM5, Ind. Eng. Chem. Res. 52, 4, 1491–1507. [CrossRef] [Google Scholar]
 Kumar P., Thybaut J.W., Teketel S., Svelle S., Beato P., Olsbye U., Marin G.B. (2013b) SingleEvent MicroKinetics (SEMK) for Methanol to Hydrocarbons (MTH) on HZSM23, Catal. Today 215, 224–232. [CrossRef] [Google Scholar]
 Laxmi Narasimhan C.S., Thybaut J.W., Marin G.B., Martens J.A., Denayer J.F., Baron G.V. (2003a) PoreMouth Physisorption of Alkanes on ZSM22: Estimation of Physisorption Enthalpies and Entropies by Additivity Method, J. Catal. 218, 1, 135–147. [CrossRef] [Google Scholar]
 Laxmi Narasimhan C.S., Thybaut J.W., Marin G.B., Jacobs P.A., Martens J.A., Denayer J.F., Baron G.V. (2003b) Kinetic Modeling of PoreMouth Catalysis in the Hydroconversion of nOctane on PtHZSM22, J. Catal. 220, 2, 399–413. [CrossRef] [Google Scholar]
 Laxmi Narasimhan C.S., Thybaut J.W., Marin G.B., Denayer J.F., Baron G.V., Martens J.A., Jacobs P.A. (2004) Relumped SingleEvent Microkinetic Model for Alkane Hydrocracking on ShapeSelective Catalysts: Catalysis on ZSM22 Pore Mouths, Bridge Acid Sites and Micropores, Chemical Engineering Science 59, 22–23, 4765–4772. [CrossRef] [Google Scholar]
 Laxmi Narasimhan C.S., Thybaut J.W., Martens J.A., Jacobs P.A., Denayer J.F., Marin G.B. (2006) A Unified SingleEvent Microkinetic Model for Alkane Hydroconversion in Different Aggregation States on Pt/HUSY Zeolites, J. Phys. Chem. B 110, 13, 6750–6758. [CrossRef] [PubMed] [Google Scholar]
 Laxmi Narasimhan C.S., Thybaut J.W., Denayer J.F., Baron G.V., Jacobs P.A., Martens J.A., Marin G.B. (2007) Aggregation State Effects in ShapeSelective Hydroconversion, Ind. Eng. Chem. Res. 46, 25, 8710–8721. [CrossRef] [Google Scholar]
 Le Lannic K. (2006) Désulfuration profonde de résidus pétroliers. Élaboration d’un modèle cinétique, Ph.D. Thesis, École Normale Supérieure de Lyon (France). [Google Scholar]
 Levine S.E., Broadbelt L.J. (2009) Detailed Mechanistic Modeling of HighDensity Polyethylene Pyrolysis: Low Molecular Weight Product Evolution, Polymer Degradation and Stability 94, 5, 810–822. [CrossRef] [Google Scholar]
 Li G., Rabitz H. (1989) A General Analysis of Exact Lumping in Chemical Kinetics, Chem. Eng. Sci. 44, 1413–1430. [CrossRef] [Google Scholar]
 Li C., Henry C.S., Jankowski M.D., Ionita J.A., Hatzimanikatis V., Broadbelt L.J. (2004) Computational discovery of biochemical routes to specialty chemicals, Chem. Eng. Sci. 59, 5051–5060. [CrossRef] [Google Scholar]
 Libanati C. (1992) Monte Carlo Simulation of Complex Reactive Macromolecular Systems, Ph.D. Thesis, University of Delaware, Delaware (USA). [Google Scholar]
 Liguras D.K., Allen D.T. (1989) Structural models for catalytic cracking. 1. Model compound reactions, Industrial & Engineering Chemistry Research 28, 6, 665–673; Structural models for catalytic cracking. 2. Reactions of simulated oil mixtures, Industrial & Engineering Chemistry Research 28, 6, 674–683. [CrossRef] [Google Scholar]
 Linstrom P.J., Mallard W.G. (2015) NIST Chemistry WebBook, NIST Standard Reference Database Number 69, National Institute of Standards and Technology, Gaithersburg, MD (USA). [Google Scholar]
 López Abelairas M., de Oliveira L.P., Verstraete J.J. (2016) Application of Monte Carlo Techniques to LCO Gas Oil Hydrotreating: Molecular Reconstruction and Kinetic Modelling, Catal. Today 271, 188–198. [CrossRef] [Google Scholar]
 López García C. (2000) Analyse de la réactivité des composés soufrés dans les coupes pétrolières : Cinétique et modélisation de l’hydrotraitement, Ph.D. Thesis, Université Claude Bernard – Lyon (France). [Google Scholar]
 López García C., Becchi M., GrenierLoustalot M.F., Païsse O., Szymanski R. (2002) Analysis of Aromatic Sulfur Compounds in Gas Oils Using GC with Sulfur Chemiluminescence Detection and HighResolution MS, Anal. Chem. 74, 3849–3857. [CrossRef] [PubMed] [Google Scholar]
 López García C., RoyAuberger M., Chapus T., Baco F. (2003) Analysis and kinetic modeling in ULSD hydrotreating, Abstr. Pap. Am. Chem. Soc. 226, U558. [Google Scholar]
 López García C., Hudebine D., Schweitzer J.M., Verstraete J.J., Ferré D. (2010) Indepth Modeling of Gas Oil Hydrotreating: From Feedstock Reconstruction to Reactor Stability Analysis, Catal. Today 150, 279–299. [CrossRef] [Google Scholar]
 Lozano Blanco G., Thybaut J.W., Surla K., Galtier P., Marin G.B. (2006) FischerTropsch Synthesis: Development of a Microkinetic Model for Metal Catalysis, Oil & Gas Science and Technology – Rev. IFP 61, 4, 489–496. [CrossRef] [EDP Sciences] [Google Scholar]
 Lozano Blanco G., Thybaut J.W., Surla K., Galtier P., Marin G.B. (2008) SingleEvent Microkinetic Model for FischerTropsch Synthesis on IronBased Catalysts, Ind. Eng. Chem. Res. 47, 16, 5879–5891. [CrossRef] [Google Scholar]
 Lozano Blanco G., Surla K., Thybaut J.W., Marin G.B. (2011) Extension of the SingleEvent Methodology to Metal Catalysis: Application to FischerTropsch Synthesis, Oil & Gas Science and Technology – Rev. IFP 66, 3, 423–435. [CrossRef] [EDP Sciences] [Google Scholar]
 MagnéDrisch J. (1995) Cinétique des réactions d’hydrotraitement de distillats par décomposition en familles et par coupes étroites, Ph.D. Thesis, Université Pierre et Marie Curie, Paris (France). [Google Scholar]
 Mahe L., Dutriez T., Courtiade M., Thiebaut D., Dulot H., Bertoncini F. (2011) Global approach for the selection of high temperature comprehensive twodimensional gas chromatography experimental conditions and quantitative analysis in regards to sulfurcontaining compounds in heavy petroleum cuts, Journal of Chromatography A 1218, 3, 534–544. [CrossRef] [PubMed] [Google Scholar]
 Marcus R.A. (1968) Theoretical relations among rate constants, barriers, and Brönsted slopes of chemical reactions, J. Phys. Chem. 72, 891–899. [CrossRef] [Google Scholar]
 MarreroMorejón J., PardilloFontdevila E. (1999) Estimation of Pure Compound Properties using GroupInteraction Contributions, AIChE J. 45, 615–621. [CrossRef] [Google Scholar]
 Martens G.G., Froment G.F. (1999) Kinetic Modeling of Paraffins Hydrocracking based upon Elementary Steps and the Single Event Concept, in Reaction Kinetics and the Development of Catalytic Processes, G.F. Froment, KC. Waugh (eds), Elsevier Science BV, Studies in Surface Science and Catalysis 122, 333–340. [CrossRef] [Google Scholar]
 Martens G.G., Marin G.B., Martens J.A., Jacobs P.A., Baron G.V. (2000) A Fundamental Kinetic Model for Hydrocracking of C_{8} to C_{12} Alkanes on Pt/USY Zeolites, Journal of Catalysis 195, 2, 253–267. [CrossRef] [Google Scholar]
 Martens G.G., Marin G.B. (2001) Kinetics for Hydrocracking based on Structural Classes: Model Development and Application, AIChE J. 47, 7, 1607–1622. [CrossRef] [Google Scholar]
 Martens G.G., Thybaut J.W., Marin G.B. (2001) Single Event Rate Parameters for Hydrocracking of Cycloalkanes on Pt/USY Zeolites, Ind. Eng. Chem. Res. 40, 8, 1832–1844. [CrossRef] [Google Scholar]
 Martinis J.M., Froment G.F. (2006) Alkylation on Solid Acids. Part 2. SingleEvent Kinetic Modeling, Ind. Eng. Chem. Res. 45, 954–967. [CrossRef] [Google Scholar]
 Matheu D.M., Lada T.A., Green W.H, Dean A.M., Grenda J.M. (2001) Ratebased Screening of Pressuredependent Reaction Networks, Computer Physics Communications 138, 3, 237–249. [CrossRef] [Google Scholar]
 Matheu D.M., Green W.H., Grenda J.M. (2003a) Capturing pressuredependence in automated mechanism generation: Reactions through cycloalkyl intermediates, International Journal of Chemical Kinetics 35, 3, 95–119. [CrossRef] [Google Scholar]
 Matheu D.M., Dean A.M., Grenda J.M., Green W.H. (2003b) Mechanism Generation with Integrated Pressure Dependence: A New Model for Methane Pyrolysis, J. Phys. Chem. A 107, 8552–8565. [CrossRef] [Google Scholar]
 Mayeno A.N., Yang R.S.H., Reisfeld B. (2005) Biochemical Reaction Network Modeling: Predicting Metabolism of Organic Chemical Mixtures, Environmental Science and Technology 39, 5363–5371. [CrossRef] [Google Scholar]
 Miller R.L., Ettre L.S., Johansen N.G. (1983) Quantitative Analysis of Hydrocarbons by Structural Group Type in Gasolines and Distillates : II. Liquid Chromatography, J. Chromatogr. A 259, 393–412. [CrossRef] [Google Scholar]
 Mitsios M., Guillaume D., Galtier P., Schweich D. (2009) SingleEvent Microkinetic Model for LongChain Paraffin Hydrocracking and Hydroisomerization on an Amorphous Pt/SiO_{2}Al_{2}O_{3} Catalyst, Ind. Eng. Chem. Res. 48, 3284–3292. [Google Scholar]
 Mizan T.I., Hou G., Klein M.T. (1998) Mechanistic Modeling of the Hydroisomerization of High Carbon Number Waxes, AIChE Meeting New Orleans. [Google Scholar]
 Mizan T.I., Klein M.T. (1999) Computerassisted Mechanistic Modeling of nHexadecane Hydroisomerization over Various Bifunctional Catalysts, Catal. Today 50, 1, 159–172. [CrossRef] [Google Scholar]
 Mochida I., Yoneda Y. (1967a) Linear Free Energy Relationships in Heterogeneous Catalysis. I. Dealkylation of Alkylbenzenes on Cracking Catalysts, J. Catal. 7, 386–392. [CrossRef] [Google Scholar]
 Mochida I., Yoneda Y. (1967b) Linear Free Energy Relationships in Heterogeneous Catalysis. II. Dealkylation and Isomerization Reactions on Various Solid Catalysts, J. Catal. 7, 393–396. [CrossRef] [Google Scholar]
 Mochida I., Yoneda Y. (1967c) Linear Free Energy Relationships in Heterogeneous Catalysis. III. Temperature Effects in Dealkylation of Alkylbenzenes on the Cracking Catalysts, J. Catal. 7, 223–230. [CrossRef] [Google Scholar]
 Mondello L., Lewis A.C., Bartle K.D. (2002) Multidimensional Chromatography, John Wiley & Sons, Inc. [Google Scholar]
 Montgomery D., Boyd M. (1959) New Method of Hydrocarbon Structural Group Analysis, Analytical Chemistry 31, 8, 1290–1298. [CrossRef] [Google Scholar]
 Moustafa T., Froment G.F. (2003) Kinetic Modeling of Coke Formation and Deactivation in the Catalytic Cracking of Vacuum Gas Oil, Ind. Eng. Chem. Res. 42, 1, 14–25. [CrossRef] [Google Scholar]
 Németh A., Vidóczy T., Héberger K., Kúti Z., Wágner J. (2002) MECHGEN: Computer Aided Generation and Reduction of Reaction Mechanisms, J. Chem. Inf. Comput. Sci. 42, 2, 208–214. [CrossRef] [PubMed] [Google Scholar]
 Neurock M. (1992) A Computational Chemical Reaction Engineering Analysis of Complex Heavy Hydrocarbon Reaction Systems, Ph.D. Thesis, University of Delaware, Delaware (USA). [Google Scholar]
 Neurock M., Nigam A., Trauth D.M., Klein M.T. (1994) Molecular Representation of Complex Hydrocarbon Feedstocks through Efficient Characterization and Stochastic Algorithms, Chem. Eng. Sci. 49, 4153–4177. [CrossRef] [Google Scholar]
 Park T.Y., Froment G.F. (2001a) Kinetic Modeling of the MTO Process – I. Model Formulation, Ind. Eng. Chem. Res. 40, 4172–4186. [CrossRef] [Google Scholar]
 Park T.Y., Froment G.F. (2001b) Kinetic Modeling of the MTO Process – II. Experimental Results, Model Discrimination and Parameter Estimation, Ind. Eng. Chem. Res. 40, 4187–4196. [CrossRef] [Google Scholar]
 Park T.Y., Froment G.F. (2004) Analysis of Fundamental Reaction Rates in the MethanolToOlefins Process on ZSM5 as a Basis for Reactor Design and Operation, Ind. Eng. Chem. Res. 43, 3, 682–689. [CrossRef] [Google Scholar]
 Peng B. (1999) Molecular Modelling of Petroleum Processes, Ph.D. Thesis, University of Manchester  Institute of Science and Technology, Manchester (UK). [Google Scholar]
 Petrakis L., Allen D. (1987) NMR for Liquid Fossil Fuels, Analytical Spectroscopy Library – Volume 1, Elsevier. [Google Scholar]
 Petzold L. (1983) Automatic Selection of Method for Stiff and NonStiff Systems of ODEs, SIAM Journal of Scientific Computing 4, 1, 136–148. [Google Scholar]
 Pfaendtner J., Broadbelt L.J. (2008a) Mechanistic Modeling of Lubricant Degradation. 1. Structure−Reactivity Relationships for FreeRadical Oxidation, Ind. Eng. Chem. Res. 47, 9, 2886–2896. [CrossRef] [Google Scholar]
 Pfaendtner J., Broadbelt L.J. (2008b) Mechanistic Modeling of Lubricant Degradation. 2. The Autoxidation of Decane and Octane, Ind. Eng. Chem. Res. 47, 9, 2897–2904. [CrossRef] [Google Scholar]
 Pitault I., Nevicato D., Forissier M., Bernard J.R. (1994) Kinetic Model Based on a Molecular Description for Catalytic Cracking of Vacuum Gas Oil, Chem. Eng. Sci. 49, 4249–4262. [CrossRef] [Google Scholar]
 Prickett S.E., Mavrovouniotis M.L. (1997a) Construction of Complex Reaction Systems – I. Reaction Description Language, Comput. Chem. Eng. 21, 11, 1219–1235. [CrossRef] [Google Scholar]
 Prickett S.E., Mavrovouniotis M.L. (1997b) Construction of Complex Reaction Systems – II. Molecule Manipulation and Reaction Application Algorithms, Comput. Chem. Eng. 21, 11, 1237–1254. [CrossRef] [Google Scholar]
 Prickett S.E., Mavrovouniotis M.L. (1997c) Construction of Complex Reaction Systems – III. An Example: Alkylation of Olefins, Comput. Chem. Eng. 21, 12, 1325–1337. [CrossRef] [Google Scholar]
 Pyl S.P., Hou Z., Van Geem K.M., Reyniers M.F., Marin G.B., Klein M.T. (2011) Modeling the Composition of Crude Oil Fractions Using Constrained Homologous Series, Ind. Eng. Chem. Res. 50, 10850–10858. [CrossRef] [Google Scholar]
 Qian K., Rodgers R.P., Hendrickson C.L., Emmett M.R., Marshall A.G. (2001) Reading Chemical Fine Print: Resolution and Identification of 3000 NitrogenContaining Aromatic Compounds from a Single Electrospray Ionization Fourier Transform Ion Cyclotron Resonance Mass Spectrum of Heavy Petroleum Crude Oil, Energy & Fuels 15, 492–498. [CrossRef] [Google Scholar]
 Qian S.A., Li C.F., Zhang P.Z. (1984) Study of structural parameters on some petroleum aromatic fractions by ^{1}H NMR/IR and ^{13}C, ^{1}H NMR spectroscopy, Fuel 63, 268–273. [CrossRef] [Google Scholar]
 Quann R.J., Jaffe S.B. (1992) StructureOriented Lumping: Describing the Chemistry of Complex Hydrocarbon Mixtures, Ind. Eng. Chem. Res. 31, 11, 2483–2497. [CrossRef] [Google Scholar]
 Quann R.J., Jaffe S.B. (1996) Building Useful Models of Complex Reaction Systems in Petroleum Refining, Chem. Eng. Sci. 51, 10, 1615–1635. [CrossRef] [Google Scholar]
 Quann R.J. (1998) Modeling the Chemistry of Complex Petroleum Mixtures, Environmental Health Perspectives 106, 6, 1441–1448. [CrossRef] [PubMed] [Google Scholar]
 QuintanaSolórzano R., Thybaut J.W., Marin G.B., Lødeng R., Holmen A. (2005) SingleEvent MicroKinetics for Coke Formation in Catalytic Cracking, Catal. Today 107, 8, 619–629. [CrossRef] [Google Scholar]
 QuintanaSolórzano R., Thybaut J.W., Marin G.B. (2007a) A SingleEvent Microkinetic Analysis of the Catalytic Cracking of (Cyclo)Alkanes on an Equilibrium Catalyst in the Absence of Coke Formation, Chem. Eng. Sci. 62, 18–20, 5033–5038. [CrossRef] [Google Scholar]
 QuintanaSolórzano R., Thybaut J.W., Galtier P., Marin G.B. (2007b) SingleEvent MicroKinetics for Coke Formation during the Catalytic Cracking of (Cyclo)Alkane/1Octene Mixtures, Catalysis Today 127, 1, 17–30. [CrossRef] [Google Scholar]
 QuintanaSolórzano R., Thybaut J.W., Galtier P., Marin G.B. (2010) Simulation of an Industrial Riser for Catalytic Cracking in the Presence of Coking using SingleEvent MicroKinetics, Catalysis Today 150, 3–4, 319–331. [CrossRef] [Google Scholar]
 Rangarajan S., Bhan A., Daoutidis P. (2010) RuleBased Generation of Thermochemical Routes to Biomass Conversion, Ind. Eng. Chem. Res. 49, 21, 10459–10470. [CrossRef] [Google Scholar]
 Ranzi E., Faravelli T., Gaffuri P., Sogaro A. (1995) Low Temperature Combustion – AutomaticGeneration of Oxidation Reactions and Lumping Procedures, Combust. Flame 102, 179–192. [CrossRef] [Google Scholar]
 Ratkiewicz A., Truong T.N. (2003) Application of Chemical Graph Theory for Automated Mechanism Generation, J. Chem. Inf. Comput. Sci. 43, 36–44. [CrossRef] [PubMed] [Google Scholar]
 Ratkiewicz A., Truong T.N. (2006) Automated Mechanism Generation: From Symbolic Calculation to Complex Chemistry, International Journal of Quantum Chemistry 106, 1, 244–255. [CrossRef] [Google Scholar]
 Read R.C. (1976) The Enumeration of Acyclic Chemical Compounds, in Balaban A.T. (ed.), Chemical Applications of Graph Theory, Academic Press, New York. [Google Scholar]
 Revellin N., Dulot H., López García C., Baco F., Jose J. (2005) Specific nitrogen boiling point profiles of vacuum gasoils, Energy & Fuels 19, 6, 2438–2444. [CrossRef] [Google Scholar]
 Sato S. (1997) The Development of a Support Program for the Analysis of Average Molecular Structures by Personal Computer, Sekiyu Gakkaishi (in Japanese) 40, 46–51. [CrossRef] [Google Scholar]
 Sato S., Matsumura A., Urushigawa Y., Metwally M., AlMuzaini S. (1998) Structural Analysis of Weathered Oil from Kuwait’s Environment, Environment International 24, 1–2, 77–87. [CrossRef] [Google Scholar]
 Schweitzer J.M., Galtier P., Schweich D. (1999) A Single Events Kinetics Model for Hydrocracking of Paraffins in a ThreePhase Reactor, Chemical Engineering Science 54, 2441–2452. [CrossRef] [Google Scholar]
 Semenov N.N. (1954) O Nekotorykh Problemakh Khimicheskoi Kinetiki i Reaktsionnoi Sposobnosti (On Some Problems in Chemical Kinetics and Reactivity), Acad. Sci. USSR, Moscow. [Google Scholar]
 Semenov N.N. (1958) Some Problems in Chemical Kinetics and Reactivity, Princeton University Press, Princeton (NJ), Vol. I & Vol. II. [Google Scholar]
 Shahrouzi J.R. (2010) Simulation et réduction de schéma cinétique de systèmes réctionnels complexes par des méthodes stochastiques (application à l’oligomérisation), Ph.D. Thesis, Université Pierre et Marie Curie, Paris (France). [Google Scholar]
 Shahrouzi J.R., Guillaume D., Rouchon P., Da Costa P. (2008) Stochastic Simulation and Single Events Kinetic Modeling: Application to Olefin Oligomerization, Ind. Eng. Chem. Res. 47, 4308–4316. [CrossRef] [Google Scholar]
 Shannon C.E. (1948) A Mathematical Theory of Communication, Bell Syst. Tech. J. 27, 379–423, 623–656. [Google Scholar]
 Sheremata J.M., Gray M.R., Dettman H.D., McCaffrey W.C. (2004) Quantitative Molecular Representation and Sequential Optimization of Athabasca Asphaltenes, Energy & Fuels 18, 1377–1384. [CrossRef] [Google Scholar]
 Simon Y., Baronne F., Marquaire P.M. (2007) Kinetic Modeling of the Oxidative Coupling of Methane, Ind. Eng. Chem. Res. 46, 7, 1914–1922. [CrossRef] [Google Scholar]
 Song J., Raman S., Yu J., Wijaya C.D., Stephanopoulos G., Green W.H. (2003) Development of Automatic Chemical Reaction Mechanism Generation Software Using ObjectOriented Technology, Prepr. Pap.Am. Chem. Soc., Div. Fuel Chem. 48, 2, 516–517. [Google Scholar]
 SoteloBoyás R., Froment G.F. (2009) Fundamental Kinetic Modeling of Catalytic Reforming, Ind. Eng. Chem. Res. 48, 3, 1107–1119. [CrossRef] [Google Scholar]
 Speight J.G. (1970) A Structural Investigation of the Constituents of Athabasca Bitumen by Proton Magnetic Resonance Spectroscopy, Fuel 49, 76–90. [CrossRef] [Google Scholar]
 Speight J.G. (1991) The Chemistry and Technology of Petroleum, 2nd ed., CRC Press. [Google Scholar]
 Stangeland B.E. (1974) A Kinetic Model for the Prediction of Hydrocracker Yields, Industrial & Engineering Chemistry Process Design and Development 13, 1, 71–76. [CrossRef] [Google Scholar]
 Surla K., Vleeming H., Guillaume D., Galtier P. (2004) A single events kinetic model: nbutane isomerization, Chemical Engineering Science 59, 22–23, 4773–4779. [CrossRef] [Google Scholar]
 Surla K., Guillaume D., Verstraete J.J., Galtier P. (2011) Kinetic Modeling using the SingleEvent Methodology: Application to the Isomerization of Light Paraffins, Oil & Gas Science and Technology – Revue d’IFP Energies nouvelles 66, 3, 343–365. [CrossRef] [EDP Sciences] [Google Scholar]
 Susnow R.G., Dean A.M., Green W.H., Peczak P., Broadbelt L.J. (1997) RateBased Construction of Kinetic Models for Complex Systems, J. Phys. Chem. A 101, 20, 3731–3740. [CrossRef] [Google Scholar]
 Suzuki T., Itoh M., Takegami Y., Watanabe Y. (1982) Chemical Structure of TarSand Bitumens by ^{13}C and ^{1}H NMR Spectroscopic Methods, Fuel 61, 402–410. [CrossRef] [Google Scholar]
 Svoboda G.D., Vynckier E., Debrabandere B., Froment G.F. (1995) Single Event Rate Parameters for Paraffin Hydrocracking on a Pt/USY Zeolite, Industrial & Engineering Chemistry Research 34, 3793–3800. [CrossRef] [Google Scholar]
 Takegami Y., Watanabe Y., Suzuki T., Mitsudo T., Itoh M. (1980) Structural Investigation on ColumnChromatographed Vacuum Residues of Various Petroleum Crudes by ^{13}C Nuclear Magnetic Resonance Spectroscopy, Fuel 59, 253–259. [CrossRef] [Google Scholar]
 Temkin O.N., Zeigarnik A.V., Kuz’min A.E., Bruk L.G., Slivinskii E.V. (2002) Construction of the Reaction Networks for Heterogeneous Catalytic Reactions: FischerTropsch Synthesis and Related Reactions, Russian Chemical Bulletin, International Edition 51, 1, 1–36. [Google Scholar]
 Teng S.T., Williams A.D. (1994) Detailed Hydrocarbon Analysis of Gasoline by GCMS (SIPIONA), J. High Resolut. Chromatogr. 17, 469–475. [Google Scholar]
 Thybaut J.W., Marin G.B., Baron G.V., Jacobs P.A., Martens J.A. (2001) Alkene Protonation Enthalpy Determination from Fundamental Kinetic Modeling of Alkane Hydroconversion on Pt/H(US)YZeolite, J. Cat. 202, 324–339. [CrossRef] [Google Scholar]
 Thybaut J.W., Marin G.B. (2003) Kinetic Modeling of the Conversion of Complex Hydrocarbon Feedstocks by Acid Catalysis, Chem. Eng. Technol. 26, 4, 509–514. [CrossRef] [Google Scholar]
 Thybaut J.W., Choudhury I.R., Denayer J.F., Baron G.V., Jacobs P.A., Martens J.A., Marin G.B. (2009) Design of Optimum Zeolite Pore System for Central Hydrocracking of LongChain nAlkanes Based on a SingleEvent MicroKinetic Model, Topics in Catalysis 52 9, 1251–1260. [CrossRef] [Google Scholar]
 Toch K., Thybaut J.W., Marin G.B. (2015) Ethene oligomerization on NiSiO_{2}Al_{2}O_{3}: Experimental investigation and SingleEvent MicroKinetic modelling, Appl. Cat. A: General 489, 292–304. [CrossRef] [Google Scholar]
 Tomlin A.S., Turányi T., Pilling M.J. (1997) Mathematical Tools for the Construction, Investigation and Reduction of Combustion Mechanisms, in: M.J. Pilling (ed.), LowTemperature Combustion and Autoignition, Elsevier, Amsterdam. [Google Scholar]
 Toulhoat H., Hudebine D., Raybaud P., Guillaume D., Kressmann S. (2005) THERMIDOR: A New Model for Combined Simulation of Operations and Optimization of Catalysts in Residues Hydroprocessing Units, Catal. Today 109, 135–153. [CrossRef] [Google Scholar]
 Trauth D.M., Stark S.M., Petti T.F., Neurock M., Klein M.T. (1994) Representation of the Molecular Structure of Petroleum Resid through Characterization and Monte Carlo Modeling, Energy & Fuels 8, 576–580. [CrossRef] [Google Scholar]
 Tuan H.P., Janssen H.M., Cramers C.A., Kuipervan Loo E.M., Vlap H. (1995) Evaluation of the Performance of Various Universal and Selective Detectors for Sulfur Determination in Natural Gas, J. High Resolut. Chromatogr. 18, 333–342. [CrossRef] [Google Scholar]
 Ugi I., Bauer J., Bley K., Dengler A., Dietz A., Fontain E., Gruber B., Herges R., Knauer M., Reitsam K., Stein N. (1993) Computer Assisted Solution of Chemical Problems – The Historical Development and the Present State of the Art of a New Discipline of Chemistry, Angew. Chem. Int. Ed. Engl. 32, 201–227. [CrossRef] [Google Scholar]
 ValdesPerez R.E., (1992a) Algorithm to generate reaction pathways for computerassisted elucidation, Journal of Comparative Chemistry 13, 9, 1079–1088. [CrossRef] [Google Scholar]
 ValdesPerez R.E., (1992b) A necessary condition for catalysis in reaction pathways, Journal of Physical Chemistry 96, 5, 2394–2396. [CrossRef] [Google Scholar]
 ValdesPerez R.E. (1994) Human/computer interactive elucidation of reaction mechanisms: application to catalyzed hydrogenolysis of ethane, Catalysis Letters 28, 1, 79–87. [CrossRef] [Google Scholar]
 ValdesPerez R.E., Zeigarnik A.V. (1997) Interactive elucidation without programming of reaction mechanisms in heterogeneous catalysis, Journal of Molecular Catalysis A: Chemical 119, 1–3, 405–414. [CrossRef] [Google Scholar]
 Valéry E. (2002) Application de la théorie des événements constitutifs à l’hydrocraquage de paraffines Lourdes, Ph.D. Thesis, Université Claude Bernard – Lyon (France). [Google Scholar]
 Valéry E., Guillaume D., Surla K., Galtier P., Verstraete J.J., Schweich D. (2007) Kinetic Modeling of Acid Catalyzed Hydrocracking of Heavy Molecules: Application to Squalane, Ind. Eng. Chem. Res. 46, 14, 4755–4763. [CrossRef] [Google Scholar]
 Van de Vijver R., Vandewiele N.M., Vandeputte A.G., Van Geem K.M., Reyniers M.F., Green W.H., Marin G.B. (2015a) Rulebased ab initio kinetic model for alkyl sulfide pyrolysis, Chem. Eng. Sci. 278, 385–393. [CrossRef] [Google Scholar]
 Van de Vijver R., Vandewiele N.M., Bhoorasingh P.L., Slakman B.L., Khanshan F.S., Carstensen H.H., Reyniers M.F., Marin G.B., West R.H., Van Geem K.M. (2015b) Automatic Mechanism and Kinetic Model Generation for Gas and SolutionPhase Processes: A Perspective on Best Practices, Recent Advances, and Future Challenges, International Journal of Chemical Kinetics 47, 4, 199–231. [CrossRef] [Google Scholar]
 Van Geem K.M., Reyniers M.F., Marin G.B., Song J., Mattheu D.M., Green W.H. (2006) Automatic Network Generation using RMG for Steam Cracking of nHexane, AIChE J. 52, 2, 718–730. [CrossRef] [Google Scholar]
 Van Geem K.M., Hudebine D., Reyniers M.F., Wahl F., Verstraete J.J., Marin G.B. (2007) Molecular Reconstruction of Naphtha Steam Cracking Feedstocks Based on Commercial Indices, Comput. Chem. Eng. 31, 1020–1034. [CrossRef] [Google Scholar]
 Van Geem K.M., Reyniers M.F., Marin G.B. (2008) Challenges of Modeling Steam Cracking of Heavy Feedstocks, Oil & Gas Science and Technology – Rev. IFP 63, 1, 79–94. [CrossRef] [EDP Sciences] [Google Scholar]
 van Krevelen D.W. (1952) Einiger neuere Einsichten, die chemische Struktur in Steinkohlen betreffend, BrennstoffChemie 33, 260–268. [Google Scholar]
 Vandegehuchte B.D., Thybaut J.W., Martinez A., Arribas M.A., Marin G.B. (2012) nHexadecane hydrocracking SingleEvent MicroKinetics on Pt/Hbeta, Applied Catalysis A: General 441, 10–20. [CrossRef] [Google Scholar]
 Vandegehuchte B.D., Thybaut J.W., Martens J.A., Marin G.B. (2015) Maximizing nalkane hydroisomerization: the interplay of phase, feed complexity and zeolite catalyst mixing, Catalysis Science & Technology 5, 4, 2053–2058. [CrossRef] [Google Scholar]
 Vandewiele N.M., Van Geem K.M., Reyniers M.F., Marin G.B. (2012) Genesys: Kinetic model construction using chemoinformatics, Chem. Eng. Journal 207–208, 526–538. [CrossRef] [Google Scholar]
 Vendeuvre C., RuizGuerrero R., Bertoncini F., Duval L., Thiébaut D., Hennion M.C. (2005) Characterisation of MiddleDistillates by Comprehensive TwoDimensional Gas Chromatography (GC×GC): A Powerful Alternative for Performing Various Standard Analyses of MiddleDistillates, J. Chromatogr. A 1086, 21–28. [CrossRef] [PubMed] [Google Scholar]
 Verstraete J. (1997) Kinetische Studie van de Katalytische Reforming van Nafta over een PtSn/Al_{2}O_{3} Katalysator, Ph.D. Thesis, Universiteit Gent (Belgium). [Google Scholar]
 Verstraete J.J., Revellin N., Dulot H., Hudebine D. (2004) Molecular Reconstruction of Vacuum Gas Oils, Abstr. Pap. Am. Chem. Soc. 227, U1070. [Google Scholar]
 Verstraete J.J., Le Lannic K., Guibard I. (2007) Modeling FixedBed Residue Hydrotreating Processes, Chem. Eng. Sci. 62, 18–20, 5402–5408. [CrossRef] [Google Scholar]
 Verstraete J.J., Schnongs P., Dulot H., Hudebine D. (2010) Molecular Reconstruction of Heavy Petroleum Residue Fractions, Chem. Eng. Sci. 65, 304–312. [CrossRef] [Google Scholar]
 Vleduts G.E. (1963) Concerning One System of Classification and Codification of Organic Reactions, Inf. Storage Retr. 1, 117–146. [CrossRef] [Google Scholar]
 Vynckier E., Froment G.F. (1991) Modeling of the Kinetics of Complex Processes based upon Elementary Steps, in: Astarita G., Sandler S.I. (eds), Kinetic and Thermodynamic Lumping of Multicomponent Mixtures, Elsevier B.V., Amsterdam, pp. 131–161. [CrossRef] [Google Scholar]
 Warth V., BattinLeclerc F., Fournet R., Glaude P.A., Come G.M., Scacchi G. (2000) ComputerBased Generation of Reaction Mechanisms for GasPhase Oxidation, Comput. Chem. 24, 541–560. [CrossRef] [PubMed] [Google Scholar]
 Watson B.A., Klein M.T., Harding R.H. (1996) Mechanistic Modeling of nHeptane Cracking on HZSM5, Ind. Eng. Chem. Res. 35, 1506–1516. [CrossRef] [Google Scholar]
 Wauquier J.P. (1994) Crude Oil  Petroleum Products  Process Flowsheets, in Petroleum Refining, Editions Technip, Paris. [Google Scholar]
 Weekman V.W., Nace D.M. (1970) Kinetics of Catalytic Cracking Selectivity in Fixed, Moving, and Fluid Bed Reactors, AIChE J. 16, 397–404. [CrossRef] [Google Scholar]
 Wei J., Kuo J.C.W. (1969) Lumping Analysis in Monomolecular Reaction Systems. Analysis of the Exactly Lumpable System. Ind. Eng. Chem. Fundam. 8, 114–123. [CrossRef] [Google Scholar]
 Wei W., Bennett C.A., Tanaka R., Hou G., Klein M.T. (2008) Detailed Kinetic Models for Catalytic Reforming, Fuel Processing Technology 89, 4, 344–349. [CrossRef] [Google Scholar]
 Willems P.A., Froment G.F. (1988a) Kinetic Modeling of the Thermal Cracking of Hydrocarbons. 1. Calculation of Frequency Factors, Ind. Eng. Chem. Res. 27, 11, 1959–1966. [CrossRef] [Google Scholar]
 Willems P.A., Froment G.F. (1988b) Kinetic Modeling of the Thermal Cracking of Hydrocarbons. 2. Calculation of Activation Energies, Ind. Eng. Chem. Res. 27, 11, 1966–1971. [CrossRef] [Google Scholar]
 Williams R.B. (1958) Characterization of Hydrocarbons in Petroleum by Nuclear Magnetic Resonance Spectrometry, in: Symposium on composition of petroleum oils, determination and evaluation, ASTM Spec. Technol. Publ. 224, 168–194. [Google Scholar]
 Wold S., Sjöström M. (1978) Linear Free energy Relationships as Tools for Investigating Chemical Similarity. Theory and Practice, in: Chapman N.B., Shorter J. (eds), Correlation analysis in Chemistry, Plenum Press, New York (NY), pp. 1–54. [CrossRef] [Google Scholar]
 Wong H.W., Li X., Swihart M.T., Broadbelt L.J. (2004) Detailed Kinetic Modeling of Silicon Nanoparticle Formation Chemistry via Automated Mechanism Generation, J. Phys. Chem. A 108, 46, 10122–10132. [CrossRef] [Google Scholar]
 Wu Y., Zhang N. (2010) Molecular Characterization of Gasoline and Diesel Streams, Ind. Eng. Chem. Res. 49, 12773–12782. [CrossRef] [Google Scholar]
 Xue G.P., Weng H.X., Thybaut J.W., Marin G.B. (2014) Catalytic Cracking of Cycloparaffins Admixed with Olefins: 2. SingleEvent Microkinetic (SEMK) Assessment, China Petroleum Processing and Petrochemical Technology 16, 2, 84–90. [Google Scholar]
 Zeigarnik A.V., ValdésPérez R.E., Temkin O.N., Bruk L.G., Shalgunov S.I. (1997) ComputerAided Mechanism Elucidation of Acetylene Hydrocarboxylation to Acrylic Acid Based on a Novel Union of Empirical and Formal Methods, Organometallics 16, 14, 3114–3127. [CrossRef] [Google Scholar]
 Zhang Y. (1999) A Molecular Approach for Characterisation and Property Predictions of Petroleum Mixtures with Applications to Refinery Modelling, Ph.D. Thesis, University of Manchester  Institute of Science and Technology, Manchester (UK). [Google Scholar]
Cite this article as: L.P. de Oliveira, D. Hudebine, D. Guillaume and J.J. Verstraete (2016). A Review of Kinetic Modeling Methodologies for Complex Processes, Oil Gas Sci. Technol 71, 45.
All Tables
Number of possible isomers of Alkanes, Alkenes, Alkynes (Cayley, 1875; Read, 1976; Hendrickson and Parks, 1991; Bytautas and Klein 1998; Faulon et al., 2005)
All Figures
Figure 1. Illustration of an elementary reaction. 

In the text 
Figure 2. Illustration of a complex reaction. 

In the text 
Figure 3. Illustration of the relations between the thermodynamic parameters of the individual elementary steps and the thermodynamic parameters of a complex reaction. 

In the text 
Figure 4. Illustration of the evolution of the lumped kinetic models for the catalytic cracking process: (a) Weekman and Nace (1970); (b) Jacob et al. (1976); (c) Pitault et al. (1994). 

In the text 
Figure 5. Evolution of the IFPEN lumped kinetic models for atmospheric gasoil hydrotreating: (a) Bonnardot (1998); (b) López García (2000); (c) López García et al. (2010). 

In the text 
Figure 6. Evolution of the IFPEN lumped kinetic models for fixedbed residue hydroprocessing: (a) Haulle (2002); (b) Le Lannic (2006); (c) Ferreira (2009). 

In the text 
Figure 7. Illustration of the SR method for creating an asphaltene molecule (Neurock et al., 1994). 

In the text 
Figure 8. Flow diagram of the SR method with the optimization loop. (Trauth et al., 1994). 

In the text 
Figure 9. Flow diagram of the direct SRREM algorithm (a) and indirect SRREM algorithm (b) for reconstruction of petroleum VR. (de Oliveira et al., 2013a). 

In the text 
Figure 10. Gillespie’s SSA (Gillespie, 1976). 

In the text 
Figure 11. Results of stochastic and deterministic simulations for the reaction with c_{1} = 0.5 s^{−1}, c_{2} = 0.1 s^{−1}, starting from X_{A}(0) = 100, X_{B}(0) = 0. 

In the text 
Figure 12. Reaction network for the catalytic reforming of nhexane (without accounting for hydrogenolysis reactions). 

In the text 
Figure 13. Reaction pathway and activated complex (a) for a β scission step and (b) for a PCP isomerization step. L_{1} to L_{5} are substituents which are part of the activated zone, while A and B are lateral chains, which do not participate in the elementary step (Guillaume et al., 2011). 

In the text 
Figure 14. Representation of the DFRN (with 47 species), the DLRN (with only the first 20 species), and SLRN (with only the 20 species created by the ‘high probability’ pathways). 

In the text 
Figure 15. Mass percent product distributions of (a) olefins and (b) carbenium ions for the oligomerization of C4 olefins up to C16 at 67% of conversion. Feed composition is 50% isobutene, 25% 1butene and 25% 2butene. The initial molecular population is 100,000. In the two limited reaction networks, the PCP steps are not considered and the protonation/deprotonation equilibrium is also considered in the DLRN. 

In the text 
Figure 16. Twostep stochastic kinetic modeling approach. 

In the text 
Figure 17. Comparison between experimental and simulated effluent weight fractions for the hydrocarbon families (SAT – Saturates, MONO – Monoaromatics, DI – Diaromatics, TRI+ – Polyaromatics with at least three aromatic cycles) and sulfur families (BT – Benzothiophenes, DBT – Dibenzothiophenes) during hydrotreating of two different LCO gas oils. 

In the text 
Figure 18. Comparison between experimental data and simulation results for (a) residue conversion, (b) sulfur removal, and (c) product yields during hydroconversion of Ural VR at 395 °C and 410 °C. 

In the text 