110
Energy & Fuels 1991,5, 110-117
Implications of Multiplicity in Kinetic Parameters to Petroleum Exploration: Distributed Activation Energy Models C. C. Lakshmanan, M. L. Bennett, and N. White* BHP Melbourne Research Laboratories, 245-273 Wellington Road, Mulgrave 31 70, Victoria, Australia Received April 9, 1990. Revised Manuscript Received July 27, 1990
A grid search procedure has been developed to estimate the kinetic parameters of kerogen decomposition represented by the distributed activation energy models (DAEM). This procedure uses approximate analytical solutions to the time-temperature integrals of DAEM. These approximations reduce the computational efforts involved in parameter estimation substantially and allow the rapid estimation and graphical representation of the three-dimensionalerror surface of the objective function over a wide range of kinetic parameters while retaining the true nonlinearity of the model. The existence of a number of sets of kinetic parameters which describe the observed data, generally within the limits of experimental accuracy, is also shown. These multiple sets of mean activation energy of the distribution and the corresponding preexponential factor correlated to a form referred to as the "compensation effect". This has been demonstrated by using published data on kerogen decomposition and experimental data. While these multiple sets of kinetic parameters can predict the laboratory data to similar levels of accuracy, they behave differently when extrapolated to geological heating rates. This then results in significantly different predictions for timing of onset of hydrocation generation. It is concluded that kinetic models used for the extrapolation from laboratory to geological heating rates are unlikely to provide any more guidance than a range of estimates for kerogen decomposition in the best case of known thermal history. This is shown to be due to lack of resolution arising from the analysis of the laboratory data which results in nonunique kinetic parameters. A kinetic model based on a distribution of activation energies initially developed for kerogen decomposition has been shown to predict length reduction of fission tracks in apatite (a mineral used as a thermal history indicator) due to thermal annealing in the laboratory. Compensation effect is observed in this case also.
Introduction Kinetic models have been used together with appropriate thermal histories over geological time scales to estimate the timing of petroleum generation and the quantity of hydrocarbons generated from sedimentary organic matter Under geological conditions, the petroleum-generating reactions proceed for millions of years at low temperatures (60-150 "C). In order to simulate these reactions in the laboratory, a known amount of kerogen is subjected to isothermal and/or nonisothermal conditions in a suitable equipment which allows the measurement of the progressive thermal decompositionof kerogen. These experiments, known as artificial maturation experiments, are performed at high temperatures (200-600O C ) and for a short duration (seconds to hours) to compensate the slow heating over millions of years to provide a basis for the development of a kinetic model. Kinetic models based on single first-order reactions have previously been reported to be inadequate.' More recent developments have been based on finite and infinite numbers of first order reaction^."^ These models use (1) Price, L. C. J. Pet. Geol. 1983, 6(1), 5-38. (2) Tissot, B. P.; Espitalie, J. Rev. Inst. Fr. Pet. 1976, 30, 743-777. (3) Tiesot, B. P.; Pelet, R.; Ungerer, P. AAPG Bull. 1987, 71(12), 1445-1466. (4) Ungerer, P.; Pelet, R. Nature 1987,327(7), 52-54. (5) Braun, R. L.; Burnham, A. K. Energy Fuels 1987, I, 153-161. (6) Mackenzie, A. S.; Quigley, T. M. AAPG Bull. 1988, 72(4), 399-415. (7) Lakshmanan, C.; White, N. Proc. 15th. Aust. Chem. Engg. Conf., CHEMECA-87 1987, 44.1-44.1.
continuous and/or discrete distributions of activation energies to account for the effect of time and temperature on the progress of petroleum-generating reactions. One of the requirements of these models is that the parameters involved explain the extent of the reactions for both laboratory and geological heating rates. Although a large amount of work on distributed activation energy models (DAEM) has been published as cited above, a systematic analysis of the existence of nonuniqueness of kinetic parameters (compensationeffect), their significance, and limitations for the situations to be discussed has not been reported. One of the objectives of this paper is to give a detailed account of these aspects of DAEM. The other objective of this paper is to investigate if DAEM can be used to explain the time and temperature on the annealing of fission tracks. Quantitative modeling of thermal annealing of fission tracks has been used to constrain thermal history of sedimentary basins.B-'O These publications describe special mathematical procedures for the prediction of thermal annealing of fission tracks as models based on a simple first-order kinetic expression do not simulate the laboratory data. It must be emphasized that the annealing of fission tracks is a complex process yet to be fully understood and the approach adopted in (8)Gleadow, A. J. W.; Duddy, I. R.; Lovering, J. F. APEA J . (Aust. Pet. Eng. Assoc.) 1983,23, 93-102. (9) Green, P. F.; Duddy, I. R.; Gleadow, A. J. W.; Tingate, P. R.; Laslett, G. M. Chem. Ceol. (Isot. Geosci. Sec) 1986, 59, 237-253. (10) Duddy, I. R.; Green, P. F.; Laslett, C. M. Chem. Geol. 1988, 73, 25-38.
0887-0624/91/2505-0110$02.50/00 1991 American Chemical Society
Kinetic Parameters in Petroleum Exploration
Energy &Fuels, Vol. 5, No. 1, 1991 111
this work is to find a mathematical description of experimental data rather than an accurate description of the physics of annealing process. Distributed Activation Energy Model (DAEM) The distributed activation energy model (DAEM) assumes that the reaction involved form a set of infinite number of parallel first-order reactions and the activation energies of these reactions can be represented by a continuous probability density function. Though it is possible to obtain this function from experimental measurements, Gaussian distribution will be assumed throughout this paper. According to DAEM, the instantaneous residual weight of kerogen, expressed as a fraction of the initial reactive kerogen, during its decomposition (isothermal and/or non-isothermal) is given by I
\
r
Similarly the rate of decomposition of kerogen as a function of time and temperature is given by
Figure 1. F(E)versus E for laboratory and geological heating
rates.
consideration to define multiple sets. The time-temperature integral in DAEM can be analytically evaluated only for isothermal conditions. But it is possible to obtain analytical values for the approximate forms of integral for nonisothermal conditions. Using Approximate Forms of t h e Integrand There are two approximations published in literature related to combustion and flames. They are based on the initial work of Pitt'l which approximates the time-temperature integral using step functions. The difference between them is basically in the choice of certain numerical constants. Consider eq 2 describing the rate of decomposition of reactive kerogen. By use of the approximation sexp( -
6)
7
d T = R F ex,( -
")
RT
to evaluate the Arrhenius integral, the rate of decomposition can be written as Equations 1 and 2 consider all reactions with positive values of activation energy. Rt in eqs 1and 2 represents the fraction of the initial reactive kerogen present at any time t and f ( E ) in eqs 1and 2 represents the probability density function. For a Gaussian distribution
Grid Search Technique-Outline of the Methodology Numerical values of the kinetic parameters (mean activation energy, standard deviation, and the preexponential factor) are required to use DAEM for a given thermal history (laboratory or geological). As there are only three parameters to be estimated, a three-dimensional grid search procedure has been adopted. This approach retains the true nonlinearity of the model in the estimation of its parameters. A region for grid search is defined by specifying the lower the upper limits of the values of the kinetic parameters (mean activation energy, standard deviation, and the preexponential factor). The number of grid points in each direction is specified and the (sum of the squares of residual between predicted and experimental values) SSR obtained by using the parameters which characterize the grid point. This process is repeated until SSR is obtained for each grid point. Then the error surface (a three-dimensionalrepresentation of the sum of the squares of deviation) is drawn and is analyzed for its shape and the associated valley. The best set of parameters is then determined as the coordinates of that grid point which produced the lowest SSR. In this process a number of sets of parameters may be found for which the SSR is either equal or approximately equal to the SSR corresponding to the best set of kinetic parameters. It must be noted that a detailed knowledge of the sensitivities of the parameters and of the errors associated with the experimental measurements are usually required to make this subjective
mt
--dt = k o x
-
1 ex.[
-
A further approximation is made where the term To2 exp(-(E/RTo)) is neglected after comparing with the corresponding term F exp(-(E/RZ")). Note that this approximation uses the highest temperature during the time interval (0 to t). On the basis of these two approximations, the rate of decomposition can be given as a+
where F(E) is defined as
Note that F(E) is a function involving a doubleexponential term. If F(E) is now plotted against E on a linear scale, it can be observed that F(E) is zero (nearly) until a critical value of activation energy (E,) is attained and after that F(E) steeply rises to a value close to unity in an almost step like manner (see Figure 1). Effectively this will mean that the integration in eq 4 need not be performed for values lower than this critical value. The following approximations are developed by using this step like behavior of F(E). These approximations basically differ in their approach to obtain the critical activation energy (E,). (11) Pitt, G.J. Fuel 1962,42, 267-274.
Lakshmanan et al.
112 Energy & Fuels, Vol. 5 , No. 1, 1991
Approximation Due to Suuberg12 Suuberg12reported that the critical activation energy (E,) can be chosen as the point at which the double exponential has a maximum slope on a F(E)versus E plot or the point at which the double exponential has a value of 0.5 or the point at which the argument of the double exponential has a value of 0.5. He reported that choosing the point at which the argument of the double exponential has a value of 0.5 gives the best approximation of the integral and developed equations to predict the instantaneous residual weight. The following equations use this choice of critical activation energy (Ec). The details of deriving these equations will not be provided, as refs 12 and 13 can be consulted. With the following notation A = - 1- koT (5) 2 P In In A (1 + In A)
]
R, = 1 - y2 erfc ( y , )
InA
dt
- ko ex.[ q2 - ,Eo ]i
EO?
kcal/mol 34.1 60.6 35.1 61.6
k,, l / h 0.105E+13b 0.212E+21 0.1053+13 0.2123+21
if yc > 0 if y c < 0
1 erfc (u,)
SSR 0.241342 0.321342 0.547342 0.238342
Data: Digitized Rock-Eva1 data from Mackenzie and Quigley! all tables E+n or E-n represents power of 10.
50 75
i , 250
350
450
o\ 0,
o'a.ao
~
r ^
_.
550
C
(6)
Figure 2. Model predictions and experimental data: comparison
(7)
choice given in eq 5. The instantaneous residual weight at any time t and the corresponding rate of reaction can be calculated by using this choice of A and eqs 6-12.
(9)
(10)
The corresponding equations for the rate of decomposition of reactive kerogen are given by --Mt
0,
kcal/mol SezenI3 1.16 2.06 Suuberg12 1.16 2.06
TEMPERATURE,deg
the instantaneous fractional residual weight of reactive kerogen at any time can be given by
R, = '/z erfc (y,)
Table I. Comparison of Parameters Estimated by Using Approximate Integrands"
if u, > 0
2R2P
(12) The term erfc (argument) represents complementary error function corresponding to the argument.
Approximation Due to SezenI3 Sezen13approached the problem in a way similar to that of Suuberg12but his choice of critical activation energy was different. He chose the critical point as the turning point on F(E) versus E plot during the steep rise. Sezen13also obtained explicit expressions for critical activation energy as a function of preexponential factor, rate of change of temperature with time, and the highest temperature during the time interval considered. He derived the equations to predict the variation of residual weight and the rate of reaction as a function of time and temperature for the case of linear heating rate only. However, these can be modified to other situations (isothermal and nonisothermal cooling conditions with linear cooling rates). According to Sezen13 A = koT/P (13) Notice the difference between this choice of A and the (12) Suuberg, E. M. Combust. Flame 1983,50, 243-246. (13) Sezen, Y. Int. J. Heat Mass Transfer 1988, 32(6), 1319-1322.
(grid search 50 X 50 X 50).
Validity of the Approximation for Laboratory and Geological Heating Rates The approximationsgiven above were developed for coal combustion. Since the temperatures and heating rates encountered in coal combustion are considerably higher than the corresponding values for kerogen decomposition (laboratory and geological), the validity of these approximations needs to be tested within the context of kerogen decomposition studied in this work. Accordingly, plots of F(E) versus activation energy for a range of temperatures of interest to kerogen decomposition are produced (Figure 1)for both laboratory and geological heating rates (typical) which suggests that F(E) can be approximated by a step function. Notice the shift in the set of curves to the lower side of activation energy as heating rate is increased. Since the approximate form of the integrand leads to closed form analytical expressions for instantaneous residual weights and rates of decomposition of kerogen, the computational times (see Table 111) are reduced drastically compared to a procedure involving numerical integration." This allows the use of a very fine grid thus providing an opportunity to look at the error surface in detail. The methodology outlined in this paper will be demonstrated by using single and multiple heating rate experimental data. These include instantaneous and rate data. Simulations Using Instantaneous Residual Weights In order to validate the grid search procedure developed in this work, Rock-Eva1 data reported by Mackenzie and Quiglef (immature samples of Brown Limestone Formation, Gulf of Suez, heated at 25 "C/min) are used. Grid search analyses using eqs 6-10 together with Sezen's13 or Suuberg's12 approximation are performed at every grid point. The estimated kinetic parameters and processing time taken are reported in Tables I and 111, respectively. Table I shows that the two approximations result in comparable estimates of the kinetic parameters. The parameters estimated by using the approximate forms of the integrands are comparable with those estimated by using numerical integration. Further work reported in this paper uses the approximation proposed by Sezen13since there are no obvious advantages of using Suuberg's approximation.12
Kinetic Parameters i n Petroleum Exploration
I
Energy & Fuels, Vol. 5, No. 1, 1991 113
---------
Figure 5. Effect of SSR on predicted residual weight (grid search 50 X 50 X 50). 100
TOP OF 3l.nlUDOir I*SS"*EC, .........................................................
-.
P.,oi,ml
0 ~ 1 . 1 6 " C C :ma e
Figure 3. Error surface. Data from Mackenzie and Quiglef (grid search 50 X 50 X 50). u = 1.16 kcal/mol. 25
M
0
100 150 TEMPERATURE, deg C
200
250
Figure 6. Effect of SSR on predicted residual weight (see Table I1 for parameters described by various sets). 100
- 3 F CF 0 , . N ' l l C C h I A I S L Y E C , ......................................................
Q
gg 35
75-
so -
P ; 5"r,rrr
s. b.6 sEr#
i
2
3
1
5
i r ? ;wale
6
~
LO
45
50
60
55
65
Ec, kcol/mole
Figure 4. Compensation effect (grid search 50
50
X
50
X
50).
Figure 2 compares the predicted values of residual weight with the experimental values over a range of temperatures. The predicted values in this figure are calculated by using kinetic parameters corresponding to the minimum SSR. However, this grid search indicated that there are a number of sets of parameters which predict the observed residual weight equally well (SSR between and lo-*). These are tabulated along with their respective SSR in Table 11. It can be seen that there are parameters close to the corresponding values reported by Mackenzie and Quigley6 (also indicated as a footnote to Table 11) for labile kerogen but with a SSR different from the best SSR estimated in this work. The error surface (corresponding to the best estimate of standard deviation) is shown in Figure 3, using logarithmic scales (base 10) for SSR and the preexponential factor. This figure illustrates clearly that the error surface is a reasonably well-behaved surface with a distinct valley and the SSR varies over four orders of magnitude. The error surface shown in Figure 3 has 125000 (50 X 50 X 50) grid points.
Multiplicity in Kinetic Parameters-Existence and Implications A semilogarithmic plot of preexponential factor versus the mean activation energy resulting in a straight line is known as the compensation effect. Mialhe et al.14 have suggested that the compensation effect is a consequence of experimental approximations. Figure 4 represents such a plot drawn using values reported in Table 11. Strictly speaking, the compensation effect shown in Figure 4 has to be associated with a value of standard deviation equal to 1.16 kcal/mol. Similar figures can be produced for other values of standard deviation. Table IV provides kinetic parameters from a grid search using data from Mackenzie and QuigleyG which resulted in sum of squares of residuals between and This tabulation shows that the mean activation energy decreases as the ratio of standard (14) Mialhe, P.;Charles, J. P.; Khoury, A. J. Phys. D: Appl. Phys. 1988,21, 383-384.
100 150 TEMPERATURE, deg C
200
250
Figure 7. Effect of SSR on predicted residual weight (grid search 50 X 50 X 50) (see Table I1 for parameters described by various
sets).
deviation to mean activation energy increases. The straight line shown in Figure 4 is a view along the valley of the error surface presented in Figure 3 looking from the corner corresponding to low activation energy and low preexponential factor. Note that the set of kinetic parameters listed in Table I1 have unequal SSR. This indicates an uneven surface along the valley of Figure 3, but for practical purposes it is not significant. Some of the estimates from Table I1 are now chosen to illustrate the significance of multiple estimates for laboratory and geological heating rates. Seven of these sets are chosen including the two extreme points on the straight line representing the compensation effect (corresponding to the lowest and highest mean activation energy), a set corresponding to the lowest SSR (the best estimate) and a set corresponding to the highest SSR. Three additional sets are chosen. These seven sets are marked (*) on their respective rows in Table 11. Also marked (**) is a set of kinetic parameters close to that of Mackenzie and &uigley.'j Figure 5 contains the predictions based on these sets of kinetic parameters and the observed values of residual weight (wt % ) from Mackenzie and Quigley.6 The difference between the predictions are, generally, within experimental resolutions. However, when these sets of kinetic parameters are used to predict the fractional residual weight for typical geological heating rates, they separate out distinctly. Figures 6 and 7 contain predicted values of residual weight (wt %) for geological heating rates indicated in the legends. From these figures, it is clear that multiple estimates, although resulting in similar predictions for laboratory heating rates, behave differently when extrapolated to geological heating rates thus demanding caution be exercised during extrapolation. In order to illustrate the significance of multiple estimates, the onset of oil generation (top of oil window) is assumed to correspond to a 10% conversion of reactive kerogen. The dotted horizontal line in Figures 6 and 7 (top of oil window) intersects the predicted curves a t different
Lakshmanan et al.
114 Energy & Fuels, Vol. 5, No. 1, 1991 Table 11. Sets of Kinetic Parameters Providing Acceptable Sum of the Sauares of Residualsa no. E,, kcal/mol KO,l / h SSR 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34
30.0 31.0 32.0 33.1 34.1 35.1 35.1 36.1 37.1 38.2 39.2 40.2 40.2 41.2 42.2 42.2 43.3 44.3 45.3 46.3 47.3 47.3 48.4 49.4 50.4 51.4 52.4 53.5 54.5 55.5 56.5 57.6 58.6 60.6
0.471E+ll 0.1333+12 0.2223+12 0.6253+12 0.105E+13 0.176E+13 0.2953+13 0.4943+13 0.8293+13 0.2333+14 0.3913+14 0.6553+14 O.llOE+15 0.1843+15 0.3093+15 0.5183+15 0.8693+15 0.1463+16 0.4093+16 0.6873+16 0.115E+17 0.1933+17 0.3243+17 0.5433+17 0.1533+18 0.256E+ 18 0.7203+18 0.121E+19 0.2023+19 0.5693+19 0.9543+19 0.2683+20 0.4503+20 0.2123+21
*
0.339E-01* 0.3853-01 0.1063-01 0.614E-01 0.241E-02* 0.6123-01 0.9423-01 0.770342 0.2983-01 0.2483-01 0.143341 0.9573-01* 0.521E-0 1 0.127E-01 0.593E-01 0.881E-01 0.232E-01 0.3873-01 0.441E-01* 0.319E-01 0.8923-01 0.7393-01 0.371E-01 0.659E-01** 0.527E-01* 0.559341 0.771E-01 0.5743-01 0.922E-01 0.689E-01* 0.805E-01 0.8893-0 1 0.7963-01 0.882E-01*
6o
I
set 1
set 2
0-0
A0--0
set 3
20 I 30
I
LO
GLEN DAVIS KEROGEN [TYPE I! K I M M E R I D G E KEROGEN [ T Y P E I11 YA?LD?IRN C O A L I T Y P E I!i! I
1
50 60 Eo,lkcal/molei
I
70
80
Figure 8. Compensation effect. set 4
set 5
0 TEMPERATURE, deg C
Figure 9. Mahakam Delta. Rate data (grid search 50 X 50 X 50). set 6
set 7
Data: Digitized Rock-Eva1 data from Mackenzie and Quigley.6 Mackenzie and Quigley6 reported the following values for the decomposition of labile kerogen. Mean activation energy (Eo) = 49.68 kcal/mol; Standard deviation (a) = 1.194 kcal/mol; preexponential factor (KO) = 0.5693+17, h-l. These are comparable to the values indicated (**) in this table. u = 1.16 kcal/mol. *Values chosen to illustrate the significance of the multiple estimates are shown with an asterisk.
temperatures and consequently a wide (significant) range for the timing of onset of oil generation is predicted. Figures 5-7 illustrate the need for maximum accuracy in determining the laboratory kinetic data and the sensitivity of the associated errors on the extrapolation to geological heating rates. Some of these estimates can be discarded on the basis of their goodness of fit to the laboratory data when experiments are performed with minimum errors. However, it is possible to see that some of these multiple estimates reported in Figure 5 cannot be resolved any further as the difference between their predictions are well within the limit of experimental accuracy and beyond the resolution that laboratory kinetic data can provide. Based on these arguments it can be seen that the timing for the onset of oil generation, as estimated from single heating rate experimental data, can vary over a significant interval as a consequence of the multiplicity in kinetic parameters which cannot be tightly constrained and/or distinguished by the laboratory experimental data. Several nonisothermal, thermogravimetric experiments are also performed15on three samples, namely, Glen Davis kerogen (Type I), Kimmeridge kerogen (Type 11), and Yallourn coal (Type 111) using linear heating rates. These include experiments conducted at multiple heating rates for every one of the samples. Figure 8 depicts the com(15)Unpublished data obtained at BPH Melbourne Research Laboratories available on request.
Figure 10. Error surface. Rate data4(grid search 51 X 51 X 51), u = 2.76 kcal/mol.
pensation effect for these data. During the initial stages of the development of this methodology, several synthetic residual weight data, generated by using known kinetic parameters, are also used to recover input kinetic parameters. These illustrated the ability of the method to recover parameters to the required tolerance. Compensation effect is also observed in these cases.
Grid Search Using Rates of Kerogen Decomposition In some laboratory maturation experiments (e.g., Rock-Eva1analysis, a widely used experimental procedure to study kerogen decomposition under laboratory heating rates) a sample is subjected to temperatures increasing linearly with time and the rate of decomposition of the sample is measured instead of instantaneous residual weights. Therefore, it becomes necessary to estimate the kinetic parameters using rate data as well. A grid search is performed by using Sezen's appr~ximation'~ on rate data
Kinetic Parameters in Petroleum Exploration
Energy & Fuels, Vol. 5, No. 1, 1991 115
Table 111. Comparison of CPU Times for Unapproximated and Approximated Forms of the Integrands" method arid size k , ranee E, ranee u ranae CPU time, min Sezen12 50 X 50 X 50 lE10-1E21 30-80 1-5 25.40 31 X 31 X 31 1E14-1E17 35-50 1-5 5.58 1E18-1E21 55-80 1-5 0.42 5 x 20 x 20 Suuberg" 50 X 50 X 50 lE10-1E21 30-80 1-5 25.40 31 X 31 X 31 1E14-1E17 35-50 1-5 5.58 55-80 1-5 129.Wb NAG" 5 x 20 x 20 lE18-1E21
" Data: Digitized Rock-Eva1 data from Mackenzie and Quig1ey.B *Numerical integration using subroutine WlAHF of NAG.17 Data General MV/lOOOO was used in all simulations reported in this work. CPU times shown include generation of a one line summary for every combination of kinetic parameters. This summary consists of run no. u, k,, Eo,and SSR. This information can be used draw three-dimensional error surfaces.
B2
- 50
oi
50
60
I
I
I
70
80
90
100
EXPERIMENTAL RESIDUAL LENGTH, %
Figure 1 1. Model predictions and experimental data: comparison isothermal experiments grid search (50 X 50
X
50).
published by Ungerer and Pelet.4 The predicted values are plotted along with the observed data in Figure 9 and the corresponding error surface (with logarithmic scale for SSR and the preexponential factor) is plotted in Figure 10. The agreement shown in Figure 9 is satisfactory considering the additional errors introduced through digitizing. The parameters (Eo = 59.0 kcal/mol; u = 2.76 kcal/mol; k, = 0.3313-20 h-l), estimated by using the three heating rate data together, compared very well with those r e p ~ r t e d . It ~ is possible to show that multiple sets of kinetic parameters exist to describe the rate of decomposition of kerogen in a similar manner as shown for the instantaneous residual weight. A detailed analysis of the results from the grid search on rate data4 suggests that the multiplicity in kinetic parameters is better constrained but not eliminated when results from multiple heating rate experiments are used for parameter estimation.
Application of DAEM to Fission Track Annealing Apatite fission track anslysis is one of the methods used for quantitative modeling of thermal histories of sedimentary basins. This technique depends on the fact that annealing of fission tracks in minerals, like hydrocarbon generation and maturation, is a function of time and temperature and the temperature interval over which it occurs in the mineral apatite is similar to that required for the hydrocarbon generation. Qualitative and quantitative analyses including mathematical expressions to predict thermal annealing of fission tracks have been previously reported.g10J6 In the present study, DAEM is used to explain the isothermalgand nonisothermal datal0 of fission track annealing. According to DAEM, the instantaneous residual length of a track, expressed as a fraction of the initial length, during its annealing (isothermal and/or nonisothermal) is given by
(16) Laslett, G. M.; Green, P. F.; Duddy, I. R.; Gleadow, A. J. W. Chem. Geol. (Isot. Geosci. Sec) 1987, 65, 1-13. (17) NAG. Numerical Algorithm Group, Fortran Library Manual, Mark 11, Oxford, U.K., 1984.
Figure 12. Error surface. Durango apatite isothermal experiments (grid search 50 X 50 X 50), u = 6.14 kcal/mol.
I
50
60
70
,
,
I
BO
90
100
EXPERIMENTAL RESIDUAL LENGTH, %
Figure 13. Model predictions and experimental data: comparison isothermal and nonisothermal experiments (50 X 50 X 50). 100,
20
I
I
I
I
40
60
80
100
TEMPERATURE, deg C
0
Figure 14. Durango apatite. Geological heating rates: DAEM vs Laslett et a1.16 Parameters from grid search (50 X 50 X 50) isothermal data.
where f ( E ) represents the probability density function. Equation 14 is used with the grid search technique presented earlier for kerogen decomposition to estimate mean activation energy, standard deviation, and preexponential factor. Figure 11 presents the comparison between predicted residual length (using the set of parameters resulting in the lowest SSR) and the experimental values for isothermal condition^.^ The corresponding error surface is shown in Figure 12. Similar simulations are performed with isothermal experimental data combined with the nonisothermal and complex heating rate data reported by Duddy et al.l0and the results are shown in Figure 13. The
Lakshmanan et al.
116 Energy & Fuels, Vol. 5, No. 1, 1991 Table IV. Results from Grid Search Parameters Resulting in SSR run no. 0 k0 EO SSR 455 1.00 1.053+12 34.1 9.263-3 2.223+11 32.0 8.173-3 2 803 1.08 34.1 3.733-3 2 955 1.08 1.053+12 34.1 2.413-3 5 455 1.16 1.053+12 4.943+12 36.1 7.703-3 5 607 1.16 34.1 4.863-3 7 955 1.24 1.05E+ 12 4.943+12 36.1 6.953-3 8 107 1.24 1.24 3.913+13 39.2 9.903-3 8 310 41.2 6.303-3 8 462 1.24 1.843+14 4.943+12 36.1 9.523-3 10607 1.33 39.2 10 810 1.33 3.913+13 8.673-3 1.843+14 41.2 3.003-3 10962 1.33 9.293-3 8.693+14 43.3 11114 1.33 41.2 2.523-3 13462 1.41 1.843+14 1.41 8.693+14 43.3 6.613-3 13614 46.3 8.323-3 13817 1.41 6.873+15 48.4 9.233-3 13969 1.41 3.243+16 1.49 1.84E+14 41.2 4.623-3 15962 43.3 6.473-3 16 114 1.49 8.69E+ 14 1.49 6.873+15 46.3 5.603-3 16 317 48.4 4.953-3 16469 1.49 3.243+16 18 462 1.843+14 41.2 9.093-3 1.57 8.643-3 18 614 1.57 8.693+14 43.3 3.243+16 48.4 2.873-3 18969 1.57 1.213+18 19 324 1.57 53.5 8.053-3 46.3 6.723-3 21 317 6.873+15 1.65 21 317 6.873+15 46.3 6.723-3 1.65 48.4 21 469 3.243+16 1.65 2.813-3 21 621 1.533+ 17 1.65 50.4 8.833-3 51.4 9.343-3 21 672 2.563+17 1.65 21 824 1.213+18 53.5 4.533-3 1.65 55.5 9.153-3 21 976 5.693+18 1.65 23 969 3.243+16 1.73 48.4 4.633-3 24 121 1.533+17 50.4 8.923-3 1.73 51.4 8.833-3 24 172 2.56E+ 17 1.73 2.793-3 24 324 1.213+18 1.73 53.5 24 476 5.693+18 1.73 55.5 6.113-3 8.203-3 24 679 4.503+19 1.73 58.6 26 469 3.243+16 1.82 48.4 8.183-3 26 824 1.213+18 1.82 53.5 2.723-3 1.82 5.693+18 55.5 4.703-3 26 976 1.82 27 179 4.503+19 58.6 5.223-3 1.82 27 331 2.123+20 60.6 6.813-3 1.213+18 4.213-3 1.90 29 324 53.5 1.90 29 476 5.693+18 55.5 4.833-3 58.6 3.733-3 1.90 29 679 4.503+19 2.123+20 60.6 4.273-3 1.90 29 831 1.213+18 53.5 7.133-3 1.98 31 824 55.5 1.98 31 976 5.693+18 6.383-3 1.98 32 179 4.503+19 58.6 3.633-3 1.98 32 331 2.123+20 60.6 3.103-3 9.273-3 2.06 34 476 5.693+18 55.5 2.06 4.503+19 58.6 4.833-3 34 679 2.123+20 60.6 3.213-3 2.06 34 831 1.003+21 62.7 9.113-3 2.06 34 983 2.14 37 179 4.50E+ 19 58.6 7.253-3 2.14 37 331 2.123+20 60.6 4.533-3 2.14 37 483 1.003+21 62.7 9.273-3 2.22 2.123+20 39 831 60.6 6.983-3 Data: Digitized Rock-Eva1 data from Mackenzie and Quigley.6
results of the DAEM are compared with corresponding values obtained by using the quantitative analysis published by Laslett et al.l8 for a geological heating rate of 3 "C/My in Figure 14. Similar comparisons are obtained (not shown) for 0.7 and 1.0 "C/My. In all these cases the DAEM predicted a greater reduction in track length compared to the results from a previous model.16 The difference between the predicted residual lengths from these two methods corresponds to a maximum difference of approximately 8-10 "C when thermal histories are predicted from measured track lengths. This suggests that DAEM and the grid search technique developed in this
I
13
35
45
LO
50
E,, kcoilmole
Figure 15. Compensation effect. Durango apatite (grid search 50 x 50 x 50) isothermal and nonisothermal data. work may provide an alternate methodology to predict the results of thermal annealing of fission tracks. Compensation effect for fission track annealing is displayed in Figure 15 based on the results from the grid search. This supports the existence of multiplicity in parameters of the mathematical model, an inference previously reported for kerogen decomposition in this paper. It must be emphasized that the approach indicated in this work using DAEM should be regarded as a mathematical description of the observed reduction in track length and not as an accurate description of the physics of the process.
Summary and Conclusions A method using grid search has been developed for predicting the instantaneous fractional residual weight of reactive kerogen and the rate of decomposition of reactive kerogen and represents a considerable advance in the aspect of parameter estimation. This concerns the use of an explicit analytical expression (approximate)to calculate the instantaneous residual weight and the associated rate of reaction involving reactive kerogen. These analytical approximations allow the user to do a grid search rapidly on three of the kinetic parameters (mean activation energy, preexponential factor, and the standard deviation). A complete picture of the error surface can thus be obtained and the significance of the kinetic parameters and their interaction with the errors in the experimental observation can be studied in detail. The analytical approximations developed in this report are the first of their kind in the context of petroleum exploration. The methodology reported in this paper has been tested with published and unpublished data (both instantaneous and rate data) and the agreement between the predicted and the experimental data is good in both cases. Comparable sets of kinetic parameters are obtained. The existence of multiple sets of kinetic parameters and the compensation effect are established. The significance of this multiplicity on the onset of oil generation is also illustrated for geological heating rates. The analysis presented in this paper suggests that the nonuniqueness in kinetic parameters cannot be resolved fully due to lack of resolution in experimental data and, for this reason, extrapolation from laboratory to geological heating rates is unlikely to provide any more guidance than a range of estimates for kerogen decomposition in the best case of known thermal history. The distributed activation energy model has also been used to predict the reduction in track lengths during thermal annealing of fission tracks in apatite observed in the laboratory and the results suggest that DAEM may provide an alternate methodology to predict the results of thermal annealing of fission tracks in apatite. It can be concluded that the analysis presented in this paper provides means to quickly examine the error surface, existence, and the significance of multiple sets of kinetic parameters which can satisfy both laboratory and geological heating rates for two classes of problems (kerogen
Energy & Fuels 1991,5, 117-122 decomposition and fission track annealing).
A E
EO E, f(E)
function involving double exponential defined in text preexponential factor, h-' residual length of fission track at time t (expressed as a fraction) heating rate, deg K/h or deg K/My gas constant, kcal/(mol K) fraction of initial reactive kerogen remaining at time t dummy variable defined in text time. h dummy variable for time occurring in integrals, h temperature at time t , K temperature at t = 0, K dimensionless quantity defined in eq 8 dimensionless quantity defined in eq 7 dummy argument defined in text standard deviation, kcal/mol
F(E)
Acknowledgment. We acknowledge Drs. Evan Evans and Kerry Gallagher for discussing experimental and theoretical aspects of this work at various stages and Robyn Klepetko for her efforts in obtaining experimental data. The programming work done by Parthasarathy Seshadri (Monash University) is also appreciated. We acknowledge BHP Corporate Research Requirements Committee support in fully funding this work. Nomenclature dimensionless constant defined in eqs 5 and 13 activation energy, kcal/mol mean activation energy, kcal/mol critical activation energy, kcal/mol probability density function defined in text
117
k0
1,
P R Rt S
t t* T
To UC
Yc z U
Hydrocracking Boscan Heavy Oil with a Co-Mo/A1203 Catalyst Containing an H-Mordenite Zeolite Component Rwaichi J. A. Minjat Department of Chemical Engineering, University of Ottawa, Ottawa, Ontario, K l N 9B4, Canada
Marten Ternan* Energy Research LaboratorieslCANMET, Energy, Mines and Resources, Canada, Ottawa, Ontario, K l A OG1, Canada Received May 22, 1990. Revised Manuscript Received August 13, 1990 Co-Mo/A1203 catalysts for hydrocracking heavy oil and residue were modified by adding up to 20 w t % hydrogen-mordenite zeolite. The acidic sites on the external surface of the mordenite crystals were expected to increase cracking reactions. In fact, there was a slight decrease in the +525 "C resid conversion, although vanadium and nickel hydrodemetallization increased as the mordenite content of the catalyst increased. On the other hand, the pseudo turnover frequency for metals removal, i.e, the number of reactions per second per reaction site (or in this case per (nm)2),was greater for catalysts containing greater amounts of mordenite. The catalyst performance was attributed to a combination of two factors. First, both the catalyst bulk density (grams of catalyst per milliliter of reactor volume) and the catalyst specific surface area (m2/g) in pores larger than 3 nm decreased as the mordenite content increased. Hence, a smaller quantity of catalyst could be placed into the reactor and the catalyst that was in the reactor had less surface area per unit mass. Clearly the mordenite changed the structure of the alumina support, which resulted in a net decrease in the effective catalyst surface area. Second, the catalyst acidity, as measured by temperature-programmed desorption of benzofuran, increased as the mordenite component of the catalyst increased. It was concluded that the improved overall hydrodemetallization was caused by both the increased number of acidic sites of the exterior surfaces of the mordenite and the changes in catalyst pore geometry, which improved the rate of diffusion to the catalyst surface. Introduction The primary objective of hydrocracking is to decrease the size of the large molecules in the +525 "C fraction of heavy oils, bitumens, and petroleum vacuum residua. The removal of heteroatoms, sulfur, nitrogen, oxygen, and metals is also of considerable importance. The remaining unconverted +525 OC residuum, after hydrocracking, is frequently used as a feedstock for coking processes. Often
* T o whom correspondence should be addressed. Permanent address: Department of Chemical and Process Engineering, University of Dar es Salaam, Box 35131,Dar es Salaam, Tanzania.
Table I. Catalyst Chemical Composition A1203, SiOz, Moo3, COO,
H-mordenite, wt%
0 1
3 5 10 20
wt% 81.5 81.2 80.0 78.9 76.5 72.4
wt%
wt%
wt%
0.0 0.3 1.7 3.2 6.0 11.2
15.4
3.1 3.1 3.1 3.1 3.0 2.9
15.3
15.1 14.9 14.4 13.5
the concentrations of sulfur and nitrogen in the coke prohibit its use for combustion. In many cases, the metals concentration in the coke is also too large, so it cannot be used for electrode coke. The Boscan heavy oil used in this
0887-0624/91/2505-Ol17$02.50/00 Published 1991 by the American Chemical Society
