1. Introduction World economy is tightly bound to energy, and fossil fuels are still the major source of energy. Conventional crude oil reservoirs are depleting constantly, and heavy oil resources are exceeding the remaining conventional oil (Meyer and Attanasi, 2003; Thomas, 2008). In situ combustion (ISC) is accepted as a thermally efficient method of enhanced oil recovery (EOR) for unconventional heavy oil reservoirs. However, it is difficult to design and implement a successful ISC project because it is a complex EOR process further complicated by heterogeneous chemical reactions. Performing a predictable and successful ISC project directly depends on the understanding of the process and * Corresponding Author: milad.karimian@gmail.com M. Karimian et al. / A Kinetic Investigation into the In Situ Combustion … 19 reaction kinetics. Early studies of crude oil combustion were performed using a reactor or a combustion tube (Fassihi, et al., 1984; Ranjbar and Pusch, 1991; Vossoughi et al., 1983). Thermal analysis methods, including thermogravimetric analysis (TG) and differential scanning calorimetry (DSC), are other common methods to study the chemical reaction kinetics. Tadema (Tadema, 1959) implemented differential thermal analysis (DTA) to study crude oils for the first time. Since then, many investigators (Ambalae et al., 2006; Ciajolo and Barbella, 1984; Khansari et al., 2012) used different methods of thermal analysis to obtain the kinetic parameters of heavy oils. According to the results of thermal analysis (Gundogar and Kok, 2014; Mahinpey et al., 2010; Rezaei et al., 2013), the combustion process can be divided into three reaction regions known as low temperature oxidation (LTO), fuel deposition (FD), and high temperature oxidation (HTO). Some investigators (Kok, 2011a; Turta, 2013) considered FD as a sub region of LTO reactions; therefore, during LTO reaction regions, the distillation of volatile hydrocarbons, the oxidation of the light hydrocarbons, and coke formation occurs consecutively. It is recommended that these reaction regions should be separated for kinetic modeling (Karimian et al., 2016). The rate of crude oil combustion processes can be parameterized using basic kinetic equations in terms of T and α as follows: d𝛼 d𝑡 = 𝑘(𝑇)𝑓(𝛼) (1) where, t is time, and T represents the temperature; α is the extent of conversion, and f(α) stands for the reaction model; finally, k(T) is the rate constant which is represented by the Arrhenius equation: 𝑘(𝑇) = 𝐴exp (− 𝐸𝑎 𝑅𝑇) (2) where, A is the pre-exponential factor, and Ea is the activation energy; R is the gas constant. A set of f(α), Ea, and A is called a kinetic triplet. Once the kinetic triplet is determined, the kinetic model is obtained. Kök et al. (Kök et al., 2004; Kök and Gul, 2013; Kok, 2011b) characterized the combustion properties of various fossil fuels using TG/DTG, DTA, and DSC methods with different additives. They used different approaches of kinetic modeling to evaluate the average activation energy. These models generally assume a reaction model prior to the calculation of Ea (model fitting), or bypass such assumption using isoconversional (model free) kinetic models. A major limitation of model fitting methods with a single heating rate is using average values for both activation energy and preexponential factor and ignoring the variation of Ea and A with conversion. Another problem is raised by the so-called compensation effect; both α and T change simultaneously, so an inaccurate selection of the reaction model in Equation 1 is compensated by the relevant inaccurate rate constant. As a result, there are different sets of k(T) and f(α) which statistically fit the experimental data equally well (Vyazovkin and Linert, 1995; Vyazovkin and Wight, 1999). According to the results of “2000 Kinetic Project” of the International Confederation of Thermal Analysis and Calorimetry (ICTAC) (Brown et al., 2000; Burnham, 2000; Maciejewski, 2000; Roduit, 2000; Vyazovkin, 2000), single-heating-rate kinetic analyses should be avoided and multiple heating rates are recommended. Most of the recent papers concerning crude oil in situ combustion have used multiple-heating-rate experiments followed by isoconversional kinetic models to obtain a variation of activation energy with the extent of conversion (Varfolomeev et al., 2016). Gundogar and Kok 20 Iranian Journal of Oil & Gas Science and Technology, Vol. 6 (2017), No. 4 (Gundogar and Kok, 2014) characterized 6 crude oils using TG/DSC. They obtained the variation of activation energy by model-free methods; the average activation energy and pre-exponential factor were calculated via model-fitting methods. Their results showed that the average activation energy of model-free and model-fitting methods differed considerably, probably due to their basic assumptions. Bai et al. (Bai et al., 2015) analyzed the pyrolysis of oil shale, and they determined Ea by a rigid model-free method and obtained f(α) using the model fitting method of Malek. Having obtained f(α) and Ea, one can calculate the pre-exponential factor and evaluate the kinetic triplet. Kuh-E-Mond reservoir is a heavy oil reservoir located in the southwest of Iran. This reservoir is a potential candidate for in situ combustion process. The kinetic triplet of the oxidation reactions of the Kuh-E-Mond heavy crude oil was investigated to contribute to the better understanding of the reaction kinetics of the ISC. The activation energy and pre-exponential factors of the crude oil sample were calculated as a function of reaction conversion. Moreover, the reaction model was evaluated using an advanced model fitting approach proposed by Shahcheraghi et al. (Shahcheraghi et al., 2015). 2. Material and methods Crude oil was provided from Kuh-E-Mond (KEM) reservoir, located in the southwest of Iran. The chemical composition of the KEM crude oil samples is presented in Figure 1. The reservoir and crude oil characteristics are given in Table 1. The crude oil density was determined in accordance with ASTM D4052, and saturates, aromatics, resins, and asphaltene (SARA) fractions were obtained according to ASTM D-2007. Table 1 Reservoir and sample characteristics of KEM. Property Oil sample Reservoir depth from sea level (m) 1100 (top), 1200 (bottom) Bubble point pressure, Pb (kPa) 6308 Reservoir temperature, TR (K) 332 Oil specific gravity, sp. gr. 0.9806 Oil density (API°) 12.8 Oil viscosity at saturation pressure, µob (cP) 1654 a t 332 K Saturates (wt.%) 12.3 Aromatics (wt.%) 21.7 Resins (wt.%) 51.3 Asphaltene (wt.%) 14.7 Differential scanning calorimetry (DSC) was performed by NETZSCH 409PG thermal analyzer system. The oxidation of crude oil causes self-heating, which induces temperature gradients within the sample (Vyazovkin, 2015); to avoid this effect, the sample thickness was decreased to 100 µm, and high conductivity platinum crucibles were used. The sample masses were maintained nearly the same (± 0.1 mgr) for all the experiments at different heating rates. Temperature calibration was performed using different pure metals with known phase transition temperatures, and the standard error of calibration was ±1.0 K. Two calibration tests (with RbNO3 and Ag2SO4) were repeated before the measurements, and the standard error of reproducibility experiments was ±0.3 K. M. Karimian et al. / A Kinetic Investigation into the In Situ Combustion … 21 The samples were heated up from room temperature (298 K) to 973 K in a dynamic oxidizing atmosphere, consisting 50 ml/min dry air and 10 ml/min N2 as the protective gas. According to ICTAC recommendations (Vyazovkin et al., 2014), a reliable kinetic experiment should consist three to five different heating rates, all less than 20 K/min; therefore, four different linear heating rates (β = 2, 5, 10, and 15 K/min) were selected. Figure 1 Components of KEM heavy oil sample. 3. Theory 3.1. Theoretical background of combustion kinetics Generally, the rate of the oxidation reaction can be described by Equation 1. When the reaction progress is measured by DSC, α is evaluated as follows: 𝛼 = ∫ (d𝐻/d𝑡)d𝑡 𝑡 𝑡0 ∫ (d𝐻/d𝑡)d𝑡 𝑡𝑓 𝑡0 = ∆𝐻 ∆𝐻𝑡𝑜𝑡 (3) where, ∆H is the current heat change, and ∆Htot is the total released heat of the process measured by DSC. Various integral isoconversional methods have been proposed to estimate the activation energy. They can be generally categorized as: (a) rigid and (b) flexible integral methods. Methods of Ozawa and Flynn and Wall (OFW) (Flynn and, Wall 1966; Ozawa, 1965) were the first traditional integral isoconversional methods addressed as rigid integral methods. The integration of Equation 1 for isothermal conditions yields: 𝑔(𝛼) = ∫ d𝛼 𝑓(𝛼) = 𝐴𝑒𝑥𝑝 (− 𝐸 𝑅𝑇) 𝑡 𝛼 0 (4) where, g(α) is the integral form of the reaction model. Solving Equation 4 requires assuming a reaction model, i.e. f(α). In complex heterogeneous reactions like crude oil combustion, many models oversimplify the real reaction scheme, and, due to the compensation effect, they do not show significantly different results, especially in certain ranges of α. Isoconversional methods bypass such assumptions by using multiple heating rates. According to isoconversional principle, the process rate at a constant extent of conversion is only a function of temperature (Vyazovkin, 2015). If the heating 0 0.66 0.23 10.35 2.35 1.95 1.62 4 3.6 2.28 2.74 70.22 0 10 20 30 40 50 60 70 80 H2S N2 CO2 C1 C2 C3 iC4 nC4 iC5 nC5 C6 C7+ mol % Component 22 Iranian Journal of Oil & Gas Science and Technology, Vol. 6 (2017), No. 4 rate (β) remains constant during the experiment, then 𝑇 = 𝑇0 + 𝛽 × 𝑡 or 𝑑𝑇 𝑑𝑡 = 𝛽. Therefore, rearranging Equation 1 and integrating yields: 𝑔(𝛼) = 𝐴𝐸 𝛽𝑅 ∫ exp(−𝑥) 𝑥 2 𝑑𝑥 = 𝐴𝐸 𝛽𝑅 𝑝(𝑥) ∞ 𝑥 (5) The term p(x) in Equation 5 has no analytical solution, and different approximations were proposed to solve it. OFW method is based on Doyle approximation (Doyle, 1962) of p(x); therefore, Equation 5 can be rearranged to obtain Equation 6. 𝑙𝑛(𝛽𝑖 )𝛼 = 𝑙𝑛 ( 𝐴𝐸𝛼 𝑔(𝛼)𝑅 ) − 5.330 − 1.052 𝐸𝛼 𝑅𝑇𝛼,𝑖 (6) where, i indicates heating rates. Eα is estimated from the slope of the linear plot of ln(βi) against 1/Tα,i. The slope at different conversions shows the dependence of Eα on α. The “rigid” integral methods estimate activation energies according to some assumptions (Vyazovkin, 2015). Vyazovkin has proposed a nonlinear isoconversional method (Vyazovkin, 1996; Vyazovkin, 2001) which performs integration over small time segments. The major advantage of Vyazovkin’s method is that there is no need to make any assumptions, and the user has a complete control over the integration; therefore, it is called “flexible”. However, the solution of flexible methods can only be obtained using computer algorithms. Vyazovkin’s method assumes that, at each extent of conversion, g(α) remains unaffected at different heating rates. For n different nonisothermal experiments, Vyazovkin’s method can be rewritten as: 𝛷(𝐸𝛼) = ∑∑ 𝐼(𝐸𝛼, 𝑇𝛼,𝑖)𝛽𝑗 𝐼(𝐸𝛼, 𝑇𝛼,𝑗)𝛽𝑖 𝑛 𝑗≠1 𝑛 𝑖=1 (7) 𝐼(𝐸𝛼, 𝑇𝛼,𝑖) = ∫ exp (− 𝐸𝛼 𝑅𝑇𝛼 ) 𝑑𝑇 𝑇𝛼 𝑇𝛼−∆𝛼 (8) where, βi and βj represent the different heating rates, and Tα-∆α and Tα are the reaction temperatures at extents of conversion of α-∆α and α respectively. Eα can be found as a value that minimizes Eα in Equation 7. Authors provided a code using Visual Basic for Applications (VBA) within Microsoft Excel® to evaluate Eα at different conversion increments. 3.2. Modified distributed activation energy model (DAEM) The distributed activation energy model was originally proposed by Vand (Vand, 1943). The model assumes that many irreversible first-order parallel reactions with different rate parameters occur simultaneously in the determined activation energies (Miura, 1995). The basic equation of DAEM can be written as: 1 − 𝛼 = ∫ exp (− 𝐴 𝛽 ∫ 𝑒 −𝐸/𝑅𝑇𝑑𝑇 𝑇 0 ) 𝑓(𝐸)d𝐸 ∞ 0 (9) where, f(Ea) is the distribution function of the activation energy to account for differences in the activation energies of many first-order irreversible reactions. Miura and Maki (Miura and Maki, 1998) simplified the DAEM to directly evaluate both Eα and Aα at different conversions. The simplified DAEM may be written by the following expression: M. Karimian et al. / A Kinetic Investigation into the In Situ Combustion … 23 ln ( 𝛽 𝑇 2 ) 𝛼 = ln ( 𝐴𝛼𝑅 𝐸𝛼 ) + 0.6075 − 𝐸𝛼 𝑅𝑇 (10) Therefore, the plot of ln(β/T 2 ) versus 1/T at each conversion yields a straight line the slope and intercept of which can be used to determine activation energy values and ln(A) respectively. The activation energy distribution function, f(Eα), can be obtained by numerically differentiating α versus Eα relationship. More details about the simplification procedure can be found elsewhere (Huang et al., 2014; Sonobe and Worasuwannarak, 2008). 3.3. Methodology of determining the advanced reaction model Model fitting methods fall into two basic categories; some methods first guess a reaction model, and they then fit experimental data to that model and finally obtain activation energy and pre-exponential factor (Vyazovkin, 2000). Due to the compensation effect, an inaccurate determination of activation energy is compensated by the erroneous pre-exponential factor. Thus, several reaction models may reasonably fit to the data. Alternatively, some models estimate Eα and A using model-free approaches, and they then estimate reaction model using master plots and comparing the normalized experimental data with theoretical models. The best match is considered as the conversion function of the reaction (Cai et al., 2014; Chen et al., 2015). Shahcheraghi et al. (Shahcheraghi et al., 2015) proposed a new method for the determination of reaction model. They obtained an overall change of the reaction rate with change of conversion or d(dα/dt)/dα. Differentiating Equation 1 with respect to conversion (α) and using chain rules give: ( d 2𝛼 d𝑡 2 ) ( d𝛼 d𝑡 ) 2 = 1 𝐴𝛼 d𝐴𝛼 d𝛼 − 1 𝑅𝑇𝛼 d𝐸𝛼 d𝛼 + 𝐸𝛼 𝑅𝑇𝛼 2 d𝑇𝛼 d𝛼 + 𝑓 ′ (𝛼) 𝑓(𝛼) (11) According to Equation 10, the pre-exponential factor can be rewritten as: 𝐴𝛼 = 𝐸𝛼𝛽 𝑅𝐶𝑇𝛼 2 exp ( 𝐸𝛼 𝑅𝑇𝛼 ) (12) where, C is a constant value. Substituting Aα from Equation 12 into Equation 11 and employing the chain rule of differentiation, Equation 11 can be rearranged as: 𝑓 ′ (𝛼) 𝑓(𝛼) = Sh(𝛼) = ( 1 𝑑𝛼 𝑑𝑡 )( ( 𝑑 2𝛼 𝑑𝑡 2 ) ( 𝑑𝛼 𝑑𝑡) − 𝛽 𝐸𝛼 𝑑𝐸𝛼 𝑑𝑇 + 2𝛽 𝑇𝛼 ) (13) The right-hand side of Equation 13 is evaluated using the experimental data, and the left-hand side can be obtained by the ratio of the derivative of the theoretical reaction model to the reaction model itself. Comparing theoretical Sh(α) with experimental Sh(α) at different extents of conversing reveals the appropriate reaction model. If the variation of activation energy obtained by an isoconversional method could be assumed constant over the course of the reaction, then dEα /dT is equal to zero, and Equation 13 is further simplified as: 𝑓 ′ (𝛼) 𝑓(𝛼) = Sh(𝛼) = ( 1 𝑑𝛼 𝑑𝑡 )( ( 𝑑 2𝛼 𝑑𝑡 2 ) ( 𝑑𝛼 𝑑𝑡) + 2𝛽 𝑇𝛼 ) (14) 24 Iranian Journal of Oil & Gas Science and Technology, Vol. 6 (2017), No. 4 4. Results and discussion 4.1. Results of the DSC tests Figure 2 shows the DSC curves of the combustion of the oil sample at all heating rates. Kok and Gul (Kok and Gul, 2013) mentioned that, by increasing heating rate, the peak amplitude increased, and similar reactions occurred at higher temperatures. It is evident from Figure 2 that the reaction zones shifted to the right at higher heating rates, and the peaks were amplified, whereas their domain decreased. The latter means that more heat is released in less time, and total heat flow is relatively constant under similar peaks at different heating rates. This is in accordance with the isoconversional principle (Vyazovkin, 2015). As mentioned before, the combustion reactions of crude oils are divided into three major stages known as: low temperature oxidation (LTO), fuel deposition (FD), and high temperature oxidation (HTO) regions. According to ICTAC recommendations, it is better to analyze each step of the multistep reactions separately (Vyazovkin et al., 2011). A look into Figure 2 reveals that it is not possible to assign a distinct peak of DSC to LTO and FD regions, whereas corresponding DSC peak of HTO region is quite distinct and clear. As shown in Figure 1, there are numerous fractions in crude oil. The lighter components of the oil mostly contribute to LTO reactions (Dabbous and Fulton, 1974). These components are converted to heavier oxygenated products such as alcohols, aldehydes, ketones, hydroperoxides, and carboxylic acids (Burger and Sahuquet, 1972). Later, these LTO products and heavier components of the crude oil participate in FD to form coke as a fuel for HTO reactions (Khansari, 2014). Consequently, various reactions having different thermal behaviors and occurring simultaneously result in a wavy DSC signal in LTO and FD regions. On the contrary, the composition of the reactants in the HTO region (i.e. high carbon content residue or coke) is fairly uniform; therefore, similar reactions take place in the HTO region, and the corresponding DSC peak is distinctive. It is possible to transform any peaks of the DSC into conversion degree using Equation 3, but isoconversional kinetic models imply that these peaks should be equivalent in the experiments conducted at different heating rates. Therefore, we only focus on the determination of the kinetic triplet of the HTO peaks. Table 2 presents the range, peak temperature, and the total peak area of the HTO reactions. The area under the peaks demonstrates the released heat. Figure 2 DSC curves of the combustion of KEM oil sample at the heating rates of 2, 5, 10, and 15 K min-1 . -2 0 2 4 6 8 10 325 425 525 625 725 825 925 Heat flow ( mW mg -1) Exo→ Temperature (K) 2 K/min 5 K/min 10 K/min 15 K/min M. Karimian et al. / A Kinetic Investigation into the In Situ Combustion … 25 Table 2 Range, peak temperature, and the peak area of the HTO reaction region of the oil sample. Heating rate (K/min) Reaction interval (K) Peak temperature (K) Peak area (J/g) 2 708-790 761 2785 5 741-833 796 2450 10 768-881 821 1841 15 775-916 838 1987 4.2. Determination of activation energy Before the implementation of isoconversional methods, one needs to find the extent of the conversion (α) versus temperature (T). The output signal of DSC was normalized to the degree of conversion using Equation 3. Typically, DSC records so many data points, so the kinetic analysis should be performed with a pre-specified step (∆α) (Vyazovkin et al., 2011). It is unlikely that the experimental α versus T contains points exactly at the selected values of α; therefore, one has to use interpolation or curve fitting to find Tα. This procedure was repeated for all the reaction regions and all the temperature programs (i.e. β values in case of a linear program). Figure 3 illustrates the variation of the conversion as a function of temperature for the oil sample at different heating rates and in different reaction regions with ∆α = 0.05. After rearranging α-T data, OFW and DAEM characteristic curves were plotted according to Equations 6 10 or the HTO region. For example, Figure 4 shows OFW characteristic curve of the crude oil sample. At each extent of conversion, the slope and intercept of the fitted lines for all the heating rates were found, and the apparent activation energy and pre-exponential factor were then calculated accordingly. Activation energy was also calculated according to Vyazovkin’s method, by implementing Equations 7 8 in a code written in Visual Basic for Application (VBA). Some calculated values of the activation energy are summarized in Table 3. Figure 3 α–T curves of the HTO reaction region of the crude oil at the heating rates of 2, 5, 10, and 15 K min-1 . 0.0 0.2 0.4 0.6 0.8 1.0 650 700 750 800 850 900 Conversion ( α) Temperature (K) 2 5 10 15 26 Iranian Journal of Oil & Gas Science and Technology, Vol. 6 (2017), No. 4 Table 3 Eα (kJ/mol) of HTO reaction region of the oil sample at selected conversions. α 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 OFW 124 123 123 123 122 120 118 114 109 DAEM 143 142 142 142 141 140 137 134 128 Vyazovkin 115 114 114 114 104 99 84 77 65 Figure 1 Kinetic characteristic curves of OFW method for HTO reaction region. According to Table 3, Eα values obtained by DAEM were higher than the ones obtained by the methods of OFW and Vyazovkin. To investigate the conversion dependence of the activation energy, mean Eα at each conversion was calculated through Eα values obtained by methods of DAEM, OFW, and Vyazovkin. The trend of mean activation energy is shown in Figure 5; in this figure, the error bars indicate standard error of the mean (SEM) at each conversion (α). It is evident from Figure 5 that, based on methods of OFW, Vyazovkin, and DAEM, Eα is almost invariant with the average of Eα= 118.9± 8.7. Figure 5 Variation of mean activation energy with conversion during HTO reaction region. 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 0.00113 0.00118 0.00123 0.00128 0.00133 0.00138 ln (β) Temperature-1 (K-1 ) 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 50 70 90 110 130 150 0.0 0.2 0.4 0.6 0.8 1.0 Activation energy (kJ.mol -1 ) Conversion (α) M. Karimian et al. / A Kinetic Investigation into the In Situ Combustion … 27 4.3. Evaluation of pre-exponential factor It is possible to calculate the pre-exponential factor (A) using Equation 10 and the kinetic characteristic curve of DAEM method, e.g. Figure 4b. Activation energy is obtained from the slope of the fitted line at each conversion. Replacing Eα in the intercept of the Equation 10 yields ln(A). Values of the pre-exponential factor are presented in Figure 6. The presence of limestone reduced average ln(A) from 39.34 to 38.26. Moreover, pre-exponential factor varies with conversion; according to DAEM model, ln(A) was 39.34±0.53 min-1 , and the dependence of A on α is a bell-shaped distribution which implies using the variation of pre-exponential factor instead of a single average value. Figure 6 Variation of pre-exponential factor versus conversion (α). Figure 7 Experimental and theoretical Sh(α) curves at β= 10 K min-1 for HTO reaction region. 0 5 10 15 20 25 0.0 0.2 0.4 0.6 0.8 1.0 Pre-exponential factor, A (min -1 ) Conversion (α) ×1016 -30 -20 -10 0 10 20 30 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 Sh(α) Conversion (α) Exp. Sh(α) An SB 28 Iranian Journal of Oil & Gas Science and Technology, Vol. 6 (2017), No. 4 4.4. Reaction model estimation Shahcheraghi’s method (Shahcheraghi et al., 2015) was used to find the most probable reaction mechanism. If the difference between the maximum and minimum values of Ea was less than 20 to 30% of the average Eα, the activation energy can be considered independent of conversion and dE/dt = 0 (Sbirrazzuoli, 2007). Therefore, Equation 13 is reduced to Equation 14. The experimental Sh(α) was obtained from the DSC data of the HTO region and was compared to several theoretical Sh(α) curves (Table 4). Figure 7 shows the experimental and best fitted theoretical Sh(α) curves at the heating rate of β= 10 K/min. The Sh(α) curves at different heating rates follow the Avrami reaction model (An) with n= 3, i.e., A3= 3(1 - α)[-ln(1 - α)](2/3). Furthermore, one can use the truncated Sestak–Berggren (SB) model (α m (1 - α)n ), which is flexible enough to treat all types of the conversion dependencies via adjustable parameters (i.e. m and n). These parameters were obtained by an optimization procedure; it was found out that m= 0.6 and n= 0.7 at different heating rates. Theoretical Sh(α) curves from SB model matches the experimental data better than A3 model. Both An and SB parameters remained virtually constant at different heating rates. Furthermore, it can be concluded that the reaction model does not significantly change with heating rate, which is in accordance with isoconversional principle. Table 4 Mathematical expressions of functions Sh(α) for different reaction models. No. Model f(α) Sh(α) 1. Chemical process or mechanism Fn (n=0,1/2,2/3,1) (1 - α)n /|1 - n| -n/(1 - α) 2. Acceleratory rate equations Pn (n = 1) (1/n)α(1-n) (1 - n)/α E1 α 1/α 3. Sigmoidal rate equations An n(1 - α) ×[-ln(1 - α)](1-1/n) [1 - (1/n) + ln(1 - α)] ×[(α - 1)ln(1 - α)]-1 F1 (1 - α) 1/(α - 1) 4. Deceleratory rate equations (phase boundary reaction) R1, F0, P1 1 0 Rn (n = 2,3) (n)(1 - α)1-1/n (1 - n)/[n(1 - α)] Fn (n = 1/2,2/3) (1 - α)n /|1 - n| -n/(1 - α) 5. Truncated Sestak–Berggren SB(m, n) α m (1 - α)n (m/α) - (n/(1 - α)) 5. Conclusions The thermal behavior of Kuh-E-Mond crude oil in a high-temperature combustion reaction region was investigated using DSC technique. The variation of the activation energy and pre-exponential factor as a function of conversion, as well as conversion function, f(α), at different heating rates was obtained. Such variations suggest that the kinetic models which assign a single average Arrhenius parameter to an overall reaction interval should not be used without evidence that confirms the constancy of A and Eα. Moreover, a major limitation of the well-known model fitting methods is the prior assumption of the reaction model, which results in the subsequent inaccurate estimation of Eα and A due to compensation effect. To solve this problem, different isoconversional kinetic models were used to calculate Eα and A independent of the reaction model at several extents of conversion. Subsequently, by using an advanced method based on the variation of Arrhenius parameters, the reaction model was obtained. According to models of OFW, Vyazovkin, and DAEM, the activation energy of the oxidation M. Karimian et al. / A Kinetic Investigation into the In Situ Combustion … 29 reactions of Kuh-E-Mond heavy oil was invariant with conversion, and an average value could be used instead. Despite Eα, the pre-exponential factor varied considerably with conversion; therefore, the variation of pre-exponential factor must be used instead of a single average value. It seemed that the conversion function follows the Avrami–Erofeev (A3) model, but truncated Sestak–Berggren (SB) fits better to the experimental data. However, both models remained virtually the same over the experiments at different heating rates. It is suggested that a ramped-temperature oxidation (RTO) experiment should be conducted within a combustion reactor containing crude oil and reservoir rock and followed by an evolved gas analysis (EGA); then, by applying a proposed kinetic model, more information can be obtained about the effect of reservoir rock on crude oil oxidation reaction. Therefore, ISC can be modeled with greater accuracy and reliability. Nomenclature A Pre-exponential factor An Avrami-Errofeev model DSC: Differential scanning calorimetry DTA Differential thermal analysis dα/dt Reaction rate EGA: Evolved gas analysis EOR Enhanced oil recovery Eα Activation Energy f(Ea) the distribution function of the activation energy f(α) Reaction model FD Fuel deposition g(α) the integral form of the reaction model HTO High temperature oxidation ISC: In situ combustion k(T) Rate constant LTO Low temperature oxidation n Model order OFW Ozawa and Flynn and Wall method p(x) Temperature integral R Gas constant RTO: Ramped-temperature oxidation SB: Sestak–Berggren model Sh(a) Shahcheraghi model TG Thermogravimetric analysis β Heating rate.