Document Type: Research Paper
Authors
^{1} Ph.D. Candidate, Department of Petroleum Engineering, Faculty of Engineering, Shahid Bahonar University of Kerman, Kerman, Iran 2 Iran Society of Young Researchers, Shahid Bahonar University of Kerman, Iran
^{2} Professor, Department of Petroleum Engineering, Shahid Bahonar University of Kerman, Iran
^{3} Professor, Process Department of Chemical Engineering, Faculty of Engineering, Shahid Bahonar University of Kerman, Kerman, Iran
Abstract
Keywords
Main Subjects
. 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 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:
(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:
(2) 
where, A is the preexponential factor, and E_{a} is the activation energy; R is the gas constant. A set of f(α), E_{a}, 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 E_{a} (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 E_{a} and A with conversion. Another problem is raised by the socalled 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), singleheatingrate kinetic analyses should be avoided and multiple heating rates are recommended. Most of the recent papers concerning crude oil in situ combustion have used multipleheatingrate 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 (Gundogar and Kok, 2014) characterized 6 crude oils using TG/DSC. They obtained the variation of activation energy by modelfree methods; the average activation energy and preexponential factor were calculated via modelfitting methods. Their results showed that the average activation energy of modelfree and modelfitting 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 E_{a} by a rigid modelfree method and obtained f(α) using the model fitting method of Malek. Having obtained f(α) and E_{a}, one can calculate the preexponential factor and evaluate the kinetic triplet.
KuhEMond 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 KuhEMond heavy crude oil was investigated to contribute to the better understanding of the reaction kinetics of the ISC. The activation energy and preexponential 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 KuhEMond (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 D2007.
Table 1
Reservoir and sample characteristics of KEM.
Property 
Oil sample 
Reservoir depth from sea level (m) 
1100 (top), 1200 (bottom) 
Bubble point pressure, P_{b }(kPa) 
6308 
Reservoir temperature, T_{R }(K) 
332 
Oil speciﬁc 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 selfheating, 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 RbNO_{3} and Ag_{2}SO_{4}) were repeated before the measurements, and the standard error of reproducibility experiments was ±0.3 K.
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 N_{2} 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:
(3) 
where, ∆H is the current heat change, and ∆H_{tot} 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:
(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 rate (β) remains constant during the experiment, then or . Therefore, rearranging Equation 1 and integrating yields:
(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.
(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:
(7) 

(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 firstorder parallel reactions with different rate parameters occur simultaneously in the determined activation energies (Miura, 1995). The basic equation of DAEM can be written as:
(9) 
where, f(E_{a}) is the distribution function of the activation energy to account for differences in the activation energies of many firstorder 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:
(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 preexponential factor (Vyazovkin, 2000). Due to the compensation effect, an inaccurate determination of activation energy is compensated by the erroneous preexponential factor. Thus, several reaction models may reasonably fit to the data. Alternatively, some models estimate E_{α} and A using modelfree 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:
(11) 
According to Equation 10, the preexponential factor can be rewritten as:
(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:
(13) 
The righthand side of Equation 13 is evaluated using the experimental data, and the lefthand 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:
(14) 
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}.
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 
708790 
761 
2785 
5 
741833 
796 
2450 
10 
768881 
821 
1841 
15 
775916 
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 prespecified 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 preexponential 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}.
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.
4.3. Evaluation of preexponential factor
It is possible to calculate the preexponential 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 preexponential factor are presented in Figure 6. The presence of limestone reduced average ln(A) from 39.34 to 38.26. Moreover, preexponential 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 bellshaped distribution which implies using the variation of preexponential factor instead of a single average value.
Figure 6
Variation of preexponential factor versus conversion (α).
Figure 7
Experimental and theoretical Sh(α) curves at β= 10 K min^{1} for HTO reaction region.
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 E_{a} 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 (A_{n}) with n= 3, i.e., A_{3}= 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 A_{3 }model. Both A_{n} 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(α) 

F_{n}(n=0,1/2,2/3,1) 
(1  α)^{n}/1  n 
n/(1  α) 

P_{n} (n = 1) 
(1/n)α^{(1n)} 
(1  n)/α 
E_{1} 
α 
1/α 


A_{n} 
n(1  α) ×[ln(1  α)]^{(11/n)} 
[1  (1/n) + ln(1  α)] ×[(α  1)ln(1  α)]^{1} 
F_{1} 
(1  α) 
1/(α  1) 


R_{1}, F_{0}, P_{1} 
1 
0 
R_{n} (n = 2,3) 
(n)(1  α)^{11/n} 
(1  n)/[n(1  α)] 

F_{n} (n = 1/2,2/3) 
(1  α)^{n}/1  n 
n/(1  α) 


SB(m, n) 
α^{m}(1  α)^{n} 
(m/α)  (n/(1  α)) 
5. Conclusions
The thermal behavior of KuhEMond crude oil in a hightemperature combustion reaction region was investigated using DSC technique. The variation of the activation energy and preexponential 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 wellknown 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 reactions of KuhEMond heavy oil was invariant with conversion, and an average value could be used instead. Despite E_{α}, the preexponential factor varied considerably with conversion; therefore, the variation of preexponential factor must be used instead of a single average value. It seemed that the conversion function follows the Avrami–Erofeev (A_{3}) 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 rampedtemperature 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 
Preexponential factor 
A_{n} 
AvramiErrofeev 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(E_{a}) 
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: 
Rampedtemperature oxidation 
SB: 
Sestak–Berggren model 
Sh(a) 
Shahcheraghi model 
TG 
Thermogravimetric analysis 
β 
Heating rate 