Energy & Fuels 1992,6, 793-803
793
Determination of Kerogen Activation Energy Distribution Padmanabhan Sundararaman' Chevron Oil Field Research Company, P.O. Box 446, La Habra, California 90633-0446
Paul H. Merz Chevron Research and Technology Company, P.O.Box 1627, Richmond, California 94802-0627
Roy G. Mann Chevron Information and Technology Company, P.0 Box 42832, Houston, Texas 77242-2832 Received May 7,1992. Revised Manuscript Received July 30, 1992
An experimental procedure for determining the activation energy distribution of isolated organic matter from source rocks is described. The method uses a nonisothermalpyrolysis technique coupled with computer analysis of the data. The computer analysis allows three options for the unknown distribution: a single activation energy for homogeneous organic matter, a log normal distribution, or a discrete distribution. In the majority of cases, the discrete distribution is superior to both the single-energy treatment and log normal distribution. Estimation of the discrete distribution is mathematically an ill-posed problem, and attention is given to an algorithm that combats the numerical difficulties. The method described assumes that there is a single preexponential factor for all possible activation energies associated with the discrete distribution. This assumption can lead to erroneous results, especially for organic matter associated with a very broad activation energy distribution, i.e., type 111kerogens. In order to obtain a robust solution, the pyrolysis experimentsshould be performed at two or three widely differing heating rates. In the case of the discrete distribution, an energy spacing of 1kcal/mol or less should be used to obtain a robust solution.
Introduction An important aspect of modern petroleum exploration is the estimation of the maturity of organic matter present in source rocks of sedimentary basins. The two most important factors for estimating maturity are (1) temperature history of the sedimentary basin; and (2) activation energy distribution of the organic matter, i.e., the energy required to convert the macromolecular organic matter present in the source rocks to simpler molecules present in oil and gas. A variety of experimental procedures have been evaluated for determining the activation energy distribution of kerogen, the macromolecular organic matter of petroleum source rocks. These methods include oil evolution1 and sealed bomb pyrolysis.2 A micropyrolysis technique for estimatingthe activation energy distribution of kerogen was developed by Ungerer and Pelet3 and by Braun and B ~ r n h a m .This ~ ~ ~ method consists of pyrolyzing the kerogen in a Rock-Eval I apparatus (Girdel Corp., France) at widely differing linear heating rates. From the rates of evolution of products at differqnt temperatures, the activation energy distribution can be calculated. In this paper we describe the details of a procedure for obtainingthe activationenergy distribution of kerogens
* Address correspondence to this author. Current address: Chevron Canada Resources, 500 Fifth Avenue 5. W., Calgary, Alberta, Canada T2P OL7. (1) Campbell, J. H.; Gellegos, G.; Gregg, M. Fuel 1980,59,727-732. (2) Lewan, M. D. Philos. Tram.R. SOC.London A 1986,315,123-134. (3) Ungerer, P.; Pelet, R. Nature 1987,327, 52-54. (4) Braun, R. L.; Bumham, A. K. Energy Fuels 1987,1, 153-161. (5) Burnham,A. K.;Braun, R. L. "Kineticsof Polymer Decomposition"; Lawrence Livermore Laboratory Report UCID-21298.
including the computationalalgorithm and the estimation of confidence intervals. Emphasis is placed on the problems associated with this approach, especially in estimating the activation energy distribution of highly heterogeneous organic matter.
Principles The principles of using the activation energy distribution have been described in detail by Anthony and Howard? The use of nonisothermal kinetics for determining activation energy distribution of coals has been described by Juntgen and Van Heek.' The rate of a single first-order reaction under isothermal conditions is given by the expression dwldt = -kw which can be integrated and put in the form w(t) = woe-kt where k is the rate constant, w(t) is the weight at time t , and wo is the initial weight of the reactant. The temperature dependence of the rate constant is given by the Arrhenius expression k = Ae-E/RT where A is the Arrhenius factor and E is the activation energy of the reaction. Consequently, under isothermal conditions, w ( t ) = wo exp(-Ae-E'RTt) Kerogen is an extremely complex macromolecule made up of carbon, hydrogen, oxygen, sulfur, and small amounts (6) Anthony, D. B.; Howard, J. B. AIChE 1976,22(4), 625-656. (7) Juntgen, H.; van Heek, K. H. Fuel 1968,47(103).
0 1992 American Chemical Society
Sundararaman et al.
794 Energy & Fuels, Vol. 6, No. 6,1992
of other elements; it contains several different types of chemical bonds. Because of the heterogeneous nature of the kerogen, its decomposition cannot be described adequately by a single chemical reaction. A mathematical model for describing the decomposition of kerogen was first introduced by Tissots and fully described by Tissot and E~pitali6~ and by Tissot and Welte.Io Briefly, the model treats kerogen as a mixture of several different species, each decomposing to produce certain amounts of products. The decomposition of kerogen is treated as a set of several parallel first-order reactions with different activation energies all occurring simultaneously but independently of each other. Quantities to be estimated from the experimental data include the frequency or density of occurrence of the various activation energies and their associated Arrhenius factors. For kerogens the number of reactions can be very large and the estimation of these parameters is a monumental task. The problem is simplified if we assume that the different reactions have the same A factors and differ only in their activation energies; this assumption is made in our procedure. The total amount of product released is given by summing the contribution from each reaction. If the number of reactions and activation energies is large, then the density of E can be treated as having a continuous density function f(E). The procedure that we have developed allows for three different approaches to determine activation energy: (1)the kerogen is homogeneous and is characterized by a single activation energy; (2) the kerogen is heterogeneous and f(E) is presumed to be a log Gaussian distribution; and (3) the kerogen is heterogeneous and f(E) is an arbitrary density function. The third approach involves a discretization of the continuumof activation e n e r g i e ~and ~ - ~imposesthe fewest assumptions on the outcome of the estimation procedure. In this approach, the activation energies are typically assumed to be between 40 and 70 kcal/mol and are chosen at regularly spaced intervals of 1 or 2 kcal/mol.
min). Since dT/dt = q Equation 1can be rewritten as (3)
Integration of eq 3 over temperature results in (4)
where woi is the initial amount of the i-th species present at the initial temperature TOand 8 is a dummy variable for temperature. The ensemble is properly represented as a continuum. In order to obtain the total weight w, integration is required over the spectrum of activation energies:
In the above equation, El and E,, are lower and upper bounds for the collection of activation energies. The function f ( E )is a weight density for the amount of kerogen corresponding to activation energy E. The total amount of kerogen that cracks is given by the integral
The responsemeasured by the flame ionizationdetector is proportional to the rate at which the kerogen is cracking. The observed quantity y reported to the computer is proportional to -dw/dT; ignoring the proportionality constant, we may write
Differentiation of eq 5 with respect to T results in an equation for y Y=
Model The collection of chemical reactions is modeled as an ensemble of independent, parallel reactions, each one characterized by its activation energy E.a7 Assume that the cracking follows first-order chemical kinetics. Then, the decay rate for the set of species corresponding to the i-th activation energy is dwi/dt = -k(T)wi where, for Arrhenius temperature dependence,
(1)
k ( T ) = A @ , ) exp(-Ei/RT) (2) In the above, wi is the weight for the i-th (pseudo-)species, t is time (min), A is the preexponential for the Arrhenius relationship (min-9, R is the gas constant (kcal/K/mol), T is absolute temperature (K), and Ei is an activation energy (kcal/mol) for the i-th species. A linearly programmed temperature schedule is used in the experiments T=To+qt for some initial temperature TOand a heating rate q (K/ (8) Tmot, B. Reu. Imt. Fr. Pet. 1969,24, 470-501. (9) Tiesot, B.; Ekpitali6, J. Rev. Zmt. Fr. Pet. 1975, 30, 743-777. (10) Tiseot, B.; Welte, D. In Petroleum Formation and Occurrence; Springer-Verlag; New York, 1984; pp 683-692.
:J
W",E)f(E)dE
(7)
where
Equation 7, a Fredholm integral equation of the first kind with kemelK, expressesthe fundamental relationship between the data and the unknowns. The measured data are expressed by the function y, and our task is to invert eq 7 and solve for f = f ( E ) . From this point on, the preexponential A is treated as a single unknown constant, independent of energy E, as previously mentioned. Inversion of eq 7 is an ill-posed problem-the solution is a discontinuous function of the data. The difficulty can be appreciated by recognizing that any high-frequency components of f are highly attenuated by the linear operator in (71, thereby perturbingy verylittle. As aresult, small errors in the data typically alter the solution by very large amounts unless a special computational procedure is used specifically to address the ill-posed nature of the problem. Lerche" and Liu and L e r ~ h e provide ~ ~ J ~ a thorough formulation of the least-squares problem of recovering Geol. 1991,23, 137-166. (12) Liu, J.; Lerche, I. Math. Geol. 1990, 22, 989-1009. (13) Liu, J.; Lerche, I. Math. Geol. 1991,23,325-347. (11) Lerche, I. Math.
Energy & Fuels, Vol. 6, No. 6, 1992 795
Kerogen Actiuation Energy Distribution
f(nas an inverse problem which requires the solution of a first-kind Fredholm integral equation. Unfortunately, the numerical results of Liu and Lerche12are an instance of the oscillations in the solution of a first-kind integral equation when the ill-posednessis not directly confronted. In particular, the numericalprocedureof Liu and Lerche12 is to use a AE spacing of approximately5 kcal/mol and to recover distributions f(E) whose histograms typically are composed of three nonzero spikes separated by zero amounts as if the kerogen were composed of three distinct components. Liu and Lerche merely assigned the amount zero whenever the least-squares treatment indicated a negative (infeasible) value for f(E). Furthermore, the spacing AE = 5 kcal/mol used by Liu and Lerche is inadequate. But, the effect of using a smaller AE with an unmodified least-squares treatment such as that of Liu and Lerche is to exacerbate the swings between positive and negative values of f ( E ) . This is the likely explanation for their recommendation for use of such a coarse AE spacing. It is very important to address the ill-posedness directly. Solving the constrained least-squares problem does properly combat the ill-posedness. OLeary and Rust14and Varah15J6discuss methods for numerically solving first-kind integral equations. Briefly, the techniques include (1)regularization, which involves introducing an additional free parameter and attempting to achieve asmooth functionf; or (2) projection techniques that restrict f to lie in a prescribed subspace; or (3) the imposition of side conditions. In this study, we have chosen the third alternative. The computational procedure of our study combatsthe effects of ill-posedness by imposing nonnegativity constraints on f as side conditions. As an additional option, constraining the integral off can yield further improvement if needed. Our method does not involvethe use of smoothness aasumptions and an ancillary Lagrange parameter, which can bias the estimation. Variants of the Model Three distinct cases have been allowed for and explored during this study: (1)the kerogen is homogeneous and characterized by a single activation energy; in this case, f is a delta function; (2) the kerogen is heterogeneous and the density f follows a lognormal distribution; (3) the kerogen is heterogeneous and no biasing assumptions are made regarding the character of the density f. A single computer program incorporates the three variants as options. Each of the three cases is discussed below. If there is just one value of E, then eqs 7 and 8 are still correct but the interpretationof f(E) is as a delta function. Rewrite eq 7 for just one energy to read
where w o is the amount of kerogen initially. Equation 9 contains three unknowns: WO,A , and E. The second case considered is that the density may follow a log normal (also called log Gaussian) distribution. Physically occurring quantities frequently follow a normal distribution; and physical quantities, which must be positive, frequently follow a log normal distribution. In this case, eq 7 can be recast as ~~~
~
~
~~
~
(14) O'hary, D.P.;Rust,
~
~
~
~
~
~
B. W.SZAM J. Sci. Stat. Comput. 1986,
7(2). (15) Varah, J. M.SIAM Reo. 1979,21(1). (16)Varah, J. M.SIAM J . Sci. Stat. Comput. 1983, 4 ( 2 ) .
where Z=lnE
In this approach, there are four unknowns to be estimated: WO,A, pz, and UZ. The newly introduced quantities, pz and UZ, denote the mean and standard deviation of Z = In E, respectively; they are measures of centrality and spread of the distribution. Each of the preceding two approaches makes biasing assumptions as to the functional form of the unknown density f. The advantage of each of the first two cases is that mathematically the problem is well-posed (the solution is a continuous function of the data), and the computation is intrinsically easier and more robust. The third case requires recovering an approximationof f without presuming anything regarding the functional form of f. A disadvantage of this approach is that the problem is mathematically ill-posed; the solution is a discontinuousfunction of the data. Specialmeasures must be taken in numerically recovering f in order to ensure that the ill-posed nature does not introduce spurious,largemagnitude oscillations in f. The oscillations result from large gain factors that a solution method may apply to the experimental error. Computation of Kernel The kernel K in (8) is fundamental to each of the three methods of investigation; it is important that the kernel be evaluated precisely. Consequently, special attention was given to the evaluation of K in order to minimize any error that might be introduced. This section describes a rearrangement of (8) that exhibits the integral so that an established algorithm for evaluation of the exponential integral can be used, and that avoids underflow on the computer. It is easier to work with the logarithm of the kernel. Define H by H=lnK so that
Under the change of variable u = 1/19,H becomes
Denote the integral on the right-hand side (RHS) by B. Break this integral into two terms according to
Apply the changes of variable u = Tu to the first integral on the RHS and u = TOUto the second integral on the RHS. Then B = T Jlmu-2e-E~IRT du - T
U-2e-Eu/RTQdu
O L
For moderate-sized values of EIRT, the quantity exp(-EIRT) can be quite small. In order to avoid
796 Energy &Fuels, Vol. 6, No. 6,1992
underflow,it is advantageousto work with a scaled version of the exponential integral, namely, ex JIm v-pe-xudv
E&)
The expression for B now reads
B = Te-%,(r)
- Toe-ToEz(ro)
where 7 7,
Sundararaman et al.
E = (E,,&, ...,EpIT and
= ( f ( E J ,f(E,),...,f(EPHT Also, define the n X p matrix K by f
Kij = c,K(Ti,Ej)
= EJRT
We have chosen the trapezoidal rule of integration so that is defined by c = (c1, c2, ...,
= EJRT,
cj = 'J,A, if j = 1, p
In the computerprogram,E2 is evaluated usinga procedure that is due to Amo~.~'-'~ Algorithm for Cases 1 and 2 Cases 1 and 2 are both mathematically well-posed. Standard nonlinear least-squares parameter estimation may be used. The computer program used here was MINPACK.,, Case 1, when the kerogen is homogeneous and there is only one activation energy to estimate, is the simplest. Approximationsto the data y are computed directly using expression 9. In case 2, there is the added complication of evaluating the integral 10 over the logarithm of activation energy. Again, special measures were taken in order to minimize any numerical error that could be introduced during the integration. Use of procedure QAG from QUADPACKZ1 allowed for evaluation of the integral 10 subject to very tight tolerances.
Algorithm for Case 3 Denote the measured data vector by T
Y = (Yl, Y,, .*.,Y,) where the yi have been measured at enough temperatures Ti to resolve the profile of the curves y = ~(5'7). More than one heating rate q is highly recommended, i.e., three or four different heating rates over as broad a range as is practical. In order to approximate the integral 7, the range of energies [E&,] is discretized and a quadrature is used. A number p > 1 of energies is chosen to resolve the profile off, and the energies are discretized with uniform spacing A according to E, = E , < E, < ... < E , = E , where
Introduce the vectors (17) Amos,D. E. ACM Trans. Math. Software 1980,6(3), 365-377. (18) Amos,D. E. ACM Trans. Math. Software 1980, 6(3), 420-448. (19) Amos, D. E. ACM Tram. Muth. Software 1983, 9(4), 525. (20) Mod, J.; Garbow, B. S.; Hillstrom, K. E. 'User Guide for MINPACK-l", Applied Mathematics Division, Argonne National Laboratory, August 1980, ANL-80-74. (21) Piessens, R.; deDoncker-Kapenga,E.;Oberhuber,C. W.; Kahaner, D. K. QUADPACK, A Subroutine Package for Automatic Integration; Springer-Verlag: New York, 1983.
cj = A, otherwise
The measured data may be at different heating rates; the value of q in eq 8 is assigned appropriately. Now the system of linear equations Kf = y (12) approximates the linear integral equation (7). For n 1 p, eq 12 could be solved for f i n the least-squares sense by solving
min IJKf- y1I2 f
where 11.( denotes Euclidean norm. The result will be a solution vector f whose entries oscillate wildly between large positive and negative values. This is due to the very poor condition number of K, which in turn is a reflection of the fact that inversion of a first-kind integral equation is ill-posed. As a class, ill-posedproblems are very difficult to solve numerically; see, for example, Varah.15916 We exploit the fact that no coordinate off can physically be negative. This side condition mitigates spurious oscillationsin f . Specifically,the computer program solves the least-squares problem min llKf - y1I2 with nonnegativity constraints
(13)
f10 Recall that K cannot be calculated unless a value for A is specified. So actually the constrained optimization (13) is embedded in a loop that is designed to choose A such that the root mean squared residual is minimized. In other words, for a given A, let r(A) be given by r(A) = Kf - y, where f solves (13) for that value of A. Then the outer loop chooses A such that Ilr(A)II is minimized. Algorithm NNLS (nonnegativeleast squares) is used to solve (13) as described by Lawson and Hanson.22 The outer loop which varies A is based on parabolic interpolation as described, for example, by Gill, Murray, and Wright.,3
Integral Constraint In the computer program, one of the first computations performed is to scale the data so the integral under the observed weight-evolved curve is unity. The measured data vector y is integrated using the trapezoidal rule according to (22) Lawson, C. L., Hanson, R. J. Solving Least Squares Problem; Prentice-Hall: New York, 1974; pp 158-184. (23) Gill, P . E.; Murray, W.; Wright, M. H. Practical Optimization; Academic Press: New York, 1981; pp 91-92.
Kerogen Activation Energy Distribution
Energy & Fuels, Vol. 6, No. 6,1992 797
“-1
n
H = 2JTJ+ 2zriV2ri
(17)
1=1
-
and the data vector is subsequently scaled according to
Yll This massaged data vector is always used regardless of any option that is chosen by the user. Scaling facilitates diagnosing possible inadequacy of the fit to the data. The magnitude of any departure of w o from unity can be a useful gauge of model inadequacy. The computer program allows an option: the user may impose the constraint wo = 1on the solution. This is simple for either of case 1or case 2 becausethe quantity wo appears explicitly in the right-hand sides of both eq 9 and 11.For case 3, the constraint is added to the least-squaresproblem as follows: Y
where J is the n X p Jacobian matrix of derivatives of residuals with respect to parameters and where V2riis the Hessian matrix of ri. Consistent with standard practice (for instance, Draper and Smith25),we approximate H using H =J ~ J (18) which is adequate when the residuals are small. Then the confidence region is approximated by the ellipsoid
9 = (bl(b- b)JTJ(b- 6 ) Ipb2F@,n- p ; a ) )
(19)
min (JKf- y(I2
(14)
The printout provides the standard errors &,b for the parameters and individual confidence intervals for the various parameters becausethese are frequently requested. The estimated standard error for parameter b k is given by
CTf= 1
(15)
and the individual confidence interval is
such that
f10 Use of this constraint (15) can further mitigate spurious oscillationsfor case 3 and numerically combats the effects of the ill-posed problem. Equation 14 is solved by use of algorithm WNNLS, which is due to Haskell and H a n ~ o n .Of ~ ~course, the solution of (14) is still embedded in a loop that iterates on the preexponential A, as described previously.
Confidence Intervals Computation of confidence intervals for the fitted parameters is performed automatically by the computer program for each of the three cases: single energy, log normal distribution, and general density function. The followingdescription covers the first two cases. Let a denote a number such as 0.05 (corresponding to 95% confidence),and let F denote Fisher’sF distribution. The 1-a inference region for the (generic)vector of parameters b can be written as
R = {b(rTrIfTf + pa2F@,n- p;a)J
(16)
where
r=y-$
n = total number of data points p = number of parameters being fitted
and the ‘hat” notation refers to variables evaluated at 6, which minimizes the sum of squared residuals rTr. A quadratic approximation to rTr is used. The p X p Hessian matrix of second derivatives of rTrwith respect to the parameters is given by (24) Haskell, K. H.; Haneon, R. J. Math. Program. 1981,21,98-118.
where t denotes Student’s distribution. The correlations among the parameters, which are also frequently requested, are summarized as a matrix C where
It is important to d e e m p h a s i ~ the e ~ ~individual confidence intervals despite the fact that they are frequently requested. The preexponential is usually strongly correlated with the activation energies. Consequently, the is usually very elongated. Under these ellipsoid circumstances,any interpretation of individualconfidence intervals as composinga ‘box” inp-dimensional parameter space whose sides are of length a b k is very misleading-such a box is much too large. Thus, the researcher is encouraged to use graphical representations of the ellipsoid (19) in lieu of a tabulation of the standard errors (20) because the common misinterpretation of (20) is to underestimate the confidence when there is high correlation among the parameters. The importanceand utility of recognitionof the high correlation between Arrhenius factor and activation energy has also been noted by Essenhigh and Misra.26 Computation of confidence intervals for the third case for a general distribution is substantially complicated by the presence of the inequality constraints that require that values of density be nonnegative. The formulation of the region R is a modified version of (16). Our procedure follows that of O’Leary and Rust.14 Comparing Models There is no simple prescription for choosing among the three models even after the computer analysis. However, some methodology can be described for this choice. One simplemeasure is to prepare a plot. As an example, results for a particular kerogen (Toarcian Shale, Germany) are presented here in order to demonstrate some means of (25) Draper, N.; Smith, H. Applied Regression Analysis; John Wiley & Sons, Inc.: New York, 1981. (26) Essenhigh, R. H.; Miera, M. K. Energy Fuels 1990,4, 171-177.
798 Energy &Fuels, Vol. 6, No.6, 1992
om, .
ffiuREIa
Sundararaman et al.
HEATINORATE q=o.39dcgKhuin
I
FIGURE lb 0.02,
HEATING RATB q = 3.9 dog wmin
1
1
TEMPERATURE (dcgK)
FIGURE IC 0.02,
1
-
HEATING RATE 9 39
FIGURE I f
Wmin I
3 l P
0.35 -'
I
LOG NORMAL AND GENERAL DISTRIBLII?ON pncnl dsmbvtlca
A
lnA=37.89
Figure 1. Comparison of fits to flame ionization detector data (a-e) and comparison of computed distributions (0. Solid line indicates general distribution and dashed line a log normal distribution. deciding between a log Gaussian distribution and a general
distribution. For this sample, pyrolysis was performed a total of five times, once at a heating rate of 0.39 K/min, twice at a heating rate of 3.9K/min, and twice at a heating rate of 39 K/min. Figure la-e exhibits the best fitsto the flame ionization detector data resulting from a single log Gaussian distribution (dashed line) and also a single general distribution (solid line). The distributions obtained from the two models (i.e., general and log normal) are plotted in Figure If. Visually, the general distribution is to be preferred overall. Looking at patterns in the residuals r is appropriate; long runs of positive or negative residuals are undesirable. Also, comparison of the sums of squared
residuals fTf can be helpful. In particular, there is the "extra s u m of squaresn analysis, which can be performed on nested models. (Refer, for example, to Bates and Watts.27) In this analysis, the single-peak model is considered to be nested within the log normal model because the first is a special case of the second-setting the variance of the log normal distribution to zero results in the single-peak model. Likewise, the log normal distribution model is nested within the generaldistribution model. (27) Bates, D. M.; Watts, D. G. Nonlinear Regression Analysis and Its Applications; John Wiley & Sons: New York, 1988. (28) Lakshmanan,C . C.; Bennett, M. L.; White, N. Energy Frtele 1991, 5, 110-117.
Energy & Fuels, Vol. 6, No. 6,1992 799
Kerogen Activation Energy Distribution Table I. Sample Extra Sum of Squama Analyeie (Toarcian Shale, Germany) sum of degrees of mean 8 0 ~ 1 ~ 8 squares freedom square Fratio P 6 6.60X 1O-g 64 >0.99999 extra 3.96X 10-5 genl distribn 2.48 X 1Oa 242 1.02 X lO-' 248 lognormal 6 . 4 4 X l W
As an example, for data set D4742203, Table I summarizesthe extra sum of squares analysis used to compare the log normal fit versus the general distribution fit. The column labeled "sum of squares" tabulates fTf; "degrees of freedom" tabulates v = n - p; "mean square" refers to fTf/v;"Fratio" is the ratio of the two mean squares; and "P"is a probability. The degrees of freedom for each of the models is the number of data points minus the number of parameters (minus one if the unit area constraintis imposed). The degreesof freedom associated with "extra" is obtained as the difference of the other two degrees of freedom. The value P in the table is obtained from that for an F(6,242) random variable where the numbers in parentheses show the degrees of freedom for the numerator and denominator. The above computation stronglypoints toward preferringthe general distribution model over the log normal distribution model for this particular kerogen since P is so large. Data Acquisition Data are acquired using a stand-alone smart data acquisition box capable of 1 p V resolution and capable of multiplexingtwo channels. One channel is assigned to an amplified and cold-junction compensated thermocouple signal, while the other is assigned to the electrometer's chart recorder output for the flame ionization detector.. The acquisition box takes one voltage reading from each channel at two points per second for temperature ramp rates of 40 and 4 K/min and at one point per second for 0.4 K/min ramps. The analog channel that represents the temperature is calibrated periodically in an explicit step. The pyrolysis temperature is ramped at 40 K/min, while the analog temperature and a digital thermocouple thermometer are sampled (sequentially as rapidly as possible) once per second. The digital temperatures are fitted by a single straight line by linear least-squares fit. Originally, we anticipated using a more complicated fit to compensate for thermocouple nonlinearity, but our error analysis indicates that the algorithm is relatively insensitive to errors of the size and type caused by this rough interpolation. Temperature calibration is performed weekly or when lost by power outage or instrument malfunction.
Data Preprocessing Once pyrolysis data for a specific sample have been acquired at several different heating rates, activation energy analysis can begin. To perform the analysis, each data file is read, smoothed (using a simple boxcar average of 15 points), normalized (offset and drift are removed, the "valid" region selected visually by the operator, and the integral of the 'valid" region is normalized to unity by simple trapezoidal rule), and sampled at an operatorspecifiedtemperature interval. In this manner, reduction of each data set to approximately 60data points enhances the speed of the analysis without sacrificingany accuracy.
50 60 70 Activation Energy (kcallmole)
40
40 50 60 70 Activation Energy (kcallmole)
Figure 2. Activation energy distributionsof type I (lacustrine) kerogens. Type I kerogensgenerallyshow singleactivationenergy illustrating the homogeneous nature of the organic matter.
Results and Discussion Activation energy distributionsfor kerogensvary widely in both shape and magnitude. Type I kerogens, i.e., kerogensderived from lacustrineorganicmatter, generally have a singleactivationenergy (Figure 2). These kerogens are made up of very uniform organic matter. Type I kerogens which have undergone degradation during deposition typically show a broader activation energy distribution. The activation energy distribution of type I1 kerogens, i.e., kerogens derived from marine organic matter, vary both in magnitude and shape (Figure3). The shape of the distribution varies from very narrow to broad and the magnitudecan vary from 48 to 55 kd/mol. Within a singlecore hole there can be wide variation in both shape and magnitude (Figure 4). Type I11 kerogens and coals generally have broad activation energy distributions (Figure 5). Braun et al.29 recently published the activation energy distributions of several type I and type I1 kerogens. Accordingto these authors, type I kerogensgenerallyshow a single activation energy at 54-56 kcal/mol,whereas type I1 kerogens have broader activation energy distributions with a node around 51-56 kcal/mol. These results are in excellent agreement with our data (Figures 2 and 3).
Energy Spacing and Multiple Minima
As described earlier,the continuumof activationenergy E is uniformly divided. The Arrhenius factor A is chosen to be the same for all values of E, and A is used as an optimizing variable. The value of In A is varied until the root mean squared discrepancy between observed and estimated rates of organic evolved is at a minimum. Ungerer and Pelet3divided the continuum of activation energy E into regular intervals of 2 kcal/mol with upper and lower limitsof 40and 80kcal/mol, respectively. Braun et al.,29on the other hand, chose 1 kcal/mol spacing. Our own computational experience on a wide variety of kerogens demonstrates that the choice of the energy spacing is important. Improper choice of entrgy spacing can lead to multiple, local minima. Plots of a(rms) vs In A for 2 and 1 kcal/mol spacings show that for 2 kcal/mol (29)Braun, R. L.;Burnham, A. K.;Reynolds J. G.;Clarbon, J. Energy Fuels 1991,5,192-204.
E.
Sundataraman et al.
800 Energy & Fuels, Vol. 6, No. 6, 1992
Activation Energy (kcallmole)
40
50
Activation Energy (kcallmole)
60
40
50
Activation Energy (kcallmole)
60
Figure 3. Activation energy distributionsof type I1 (marine) kerogens. The shape and the magnitude of type I1 kerogens are highly variable.
180
190
5 9.
205
Q,
n
210
230
50
55
Figure 4. Activation energy distributions of several kerogens from a single core hole illustrating the variability both in shape and magnitude. The variability is probably due to the variation in the organic matter during diagenesis.
spacing there can be multiple minima (Figure 6), whereas narrowingthe energyspacing to 1kcal/mol leads to a single minimum. Ungerer and Pelet3 overcame the problem of converging to a local minimum rather than the global minimum by repeating the analysis with various random values of A and f ( E )to overcome it. Burnham and Braun5 chose to search the parameter space systematically by shifting the initial convergent value of A to a predictable value for which another local minimum in the residual exists and again convergingto an optimum value of A. In our experience, 1kcal/mol spacing can usually be relied
on to produce a single minimum. Exceptions are kerogens that have homogeneous organic matter and thus have a single activation energy. In these cases (Figure 7), even use of 1kcal/mol spacing leads to multiple minima. This is anticipated mathematically because the model of a continuum of energies is not in accord with a single energy. In the case of single activation energies, it is appropriate touse a different approach in order to obtain the activation energy and frequency factor of the kerogen. Our computer procedure specifically includes an option for a homogeneous kerogen with a single value of E.
Energy & Fuels, Vol. 6, No. 6,1992 801
Kerogen Activation Energy Distribution Llgnltw
North Dakota
5
0
8
0
7
0
8
0
5
0
8
0
7
0
8
0
55
5
0
8
0
7
0
8
0
Figure 5. Activation energy distributions of three type I11 kerogens. Type I11 kerogens generally show a broad activation energy distribution. The node usually occurs at high values illustrating the refractory nature of the organic matter. 0.4 7
1
0.3
ga
0.2
0.1
1I
4
.
I
1 Ty e I1 Kerogen 1 Rcai Spacing
0.3
0.0
a 30 40 50
20
Ln (A)
Figure 6. Plots of root mean square (rms)residual against In A showing (a) that use of 2 kcal of energy spacing can lead to multiple minima each having different solution and (b) in these circumstances use of 1 kcal/mol spacing will lead to a single minimum and hence to a unique solution. Multiple Heating Rates In order to obtain satisfactory reproducibility and reliability, Ungerer and Pelet3and Braun and Burnham4v5 recommended that the analysis be performed at various heating rates with the higher and the lower heating rates differing by 2 orders of magnitude. In our experience, at least two heating rates and in some instances three are necessary to obtain a satisfactory solution. Plots of a and In A show (Figure 8a) that use of one heating rate will not lead to a unique solution. Thus, an activation energy distribution obtained from a sing!e heating rate is not reliable. Use of two widely different heating rates results in a substantial improvement in the robustness of the solution (Figure 8b). An additional analysis using a third heating rate can additionally improve the solution. Such improvement in the quality of estimation of the activation energy distribution parameters is mathematically anticipated. As an example, consider the case in which the distribution is assumed to be log normal. The confidenceregion is the multidimensional ellipsoiddefined
Figure 7. Plot of root mean square (rms) residual against In A for kerogens with very narrow activation energy distribution showingthat use of 1kcal/mol spacingleads to multiple minima, each having a different solution. by eq 19. The volume of the ellipsoid is proportional to )JTJI-1/2. Because one wants this volume to be minimum, one maximizes the value of the determinant IJTJI.This is the well-known experimental design criterion of D-optimality.27 Table 11 illustrates the benefit of performing the experimentsat widely different heating rates. The values of the determinant lJTJlare tabulated for different possible experimental designs. The left-hand column identifies the design using a three-digit code for the three heating rates 0.39, 3.9, and 39 K/min. The values of the three digits indicate the number of experiments at that heating rate. Thus, 121 indicates one run a t 0.39 K/min, two runs at 3.9 K/min, and one run at 39 K/min. It is important to recognizethat the listed entries are independent of any data on y. They correspond to particular temperature values T a t which FID measurements could be taken. In fact, in Table 11, they were based on the T values for Toarcian Shale, Germany. A glance at the table reveals the benefit conferred by multiple experiments at very different heating rates. Each of the entries for two heating rates (rows 011, 110, and 101) has values of JJTJ) that are larger than those for a single heating rate by at least a factor of 17 (compare row 110 with row 001). Further improvement by a factor of at least 2.4 can be achieved using three heating rates (comparerow 111 with 101). Note also that there is added value in replicated runs, but the advantage is not as great as that from widely different heating rates. The matrix JTJalso determines the correlation among the parameters. The table displays the correlation between the Arrhenius factor In A and the mean log energy pln E. In all cases,the correlation is very strong, larger than 0.998. These two parameters are very strongly coupled. No significant reduction in the correlation between these two parameters is obtained by either adding heating rates or repeated experiments. Decreasing (or increasing) the estimate of both of these parameters by the same small amount can explain the observed data virtually as well as the optimal parameter estimates. Lakshmanan, Bennett, and Whitez9refer to the strong correlation between mean log energy pin E and preexponential factor In A as the compensation effect. These researchers find the relationship empirically and display insightful graphs of the sum of squared residuals surface
802 Energy & Fuels, Vol. 6, No. 6,1992 0.08 1 0.07
-
Sundararaman et al. 1
0.4
0.4CImln 4.0C/min 40Cimin
0.02
0.3
-
0.2
-
0.1
-
0.0
28
30
32
34
36
38
40
Ln (A)
42
Two Heatlng Rates (40,4 C/mln)
I 30 32 34 36 38 40 42
28
Ln (A)
Figure 8. Plots of root mean square (rms) residual against In A showing (a) that use of single heating rate does not lead to a unique solution and (b) that use of two heating rates leads to a unique solution. Table 11. lJTJlValues for Toarcian Shale, Germany design det(JTJ) corr(ln A. R) 1.1x 10-8 0.9998 001 8.8 X lo+ 010 0.9998 100 0.9998 6.6 X lo4 2.2 x 10-7 011 0.9995 1.9x 10-7 110 0.9995 0.9985 6.1x 10-7 101 1.5 X lo4 0.9989 111 0.9991 121 2.7 X lo4 0.9990 221 5.7 x 104 ~
as a function of both p h E and In A. The graphs visually elucidate the correlation. These graphs show the valley in the surface, which is oriented along the approximate linear relationship between the two parameters. Lakshmanan,Bennet, and White29alsosuggestthat this diminishesthe practical value of these artificialmaturation experiments (high temperatures for short duration) for evaluation of sedimentary basins over geologic periods of time (comparativelylow temperatures for long duration). However, the disparity that they cite in predictions of geologic yield is substantially inflated by the statistically unnecessarily severe spread among the parameter estimates that is chosen as an example. Nielsen and DahP also report a high positive correlation between the preexponential factor and the activation energy. It is important to recognize that experimental error is not the source of this correlation. Rather, the c a w is the experimentaldesign as reflected by eq 22 being independent of 8. The correlation is not attributable to experimental error as Nielsen and Dahl suggest. If the correlationbetween these two parameters is to be reduced, then an entirely different experiment must be designed-a mere reduction of laboratory error will not diminish the correlation between the preexponential factor and the activation energy. The same conclusion as ours was reached previously by Essenhigh and Misra.26 A descriptionof the geometry of the ellipsoidalinference region 9 of (19) and the parameter confidence intervals (21) is helpful. When parameters are highly correlated, the ellipsoid 9 is very eccentric-long and thin. When the parameters have small correlation, the ellipsoid 9 is approximately spherical. Most of the points which fall inside the box whose vertices are determined by the confidence intervals are unreasonable and are outside (30) Nieleen, S.
B.;Dahl, B.Mar. Pet. Geol. 1991, 8, 483-492.
?( when the parameters are highly correlated. Thus, the confidenceinterval is useful when a range for an individual parameter is estimated. But, when strongly correlated, multiple parameters are consideredjointly, the confidence intervals are misleading and the ellipsoid should be used. When multiple parameters are considered, confidence intervals are deemphasized. If an experimental design is kept the same but laboratory procedures are improved so that error is reduced, then the effect is to diminish the size (volume) of the confidence ellipsoid 3; however, the expected correlations among the parameters and the eccentricity of the ellipsoid remain the same.
Kerogens with Broad Activation Energies In the calculation of.activation energy using a discrete distribution model, the problem is simplified if it is assumed that the Arrhenius factor A is the same for all the reactions. This has been justified because of the limited variation in the A factor3 and its small influence on the rate constant, compared to the influence of temperature and activation energy. However, our experiments have shown that the simplifyingassumption is not alwaysvalid, particularlyin the case of kerogens with broad activation energy distributions such as coals and vitriniterich kerogens. Figure 9 shows the activation energy distributions of pure liptinite and pure vitrinite as well as those of various mixtures of these two macerals. Pure liptinite has a very narrow activation energy distribution centered around 52 kcal/mol, whereas pure vitrinite has a broad activation energy distribution centered around 69 kcal/mol. Kerogen is a physical mixture of different macerals. Each of these macerals decomposes at different rates depending on its activation energy distribution and frequency factor. The rates are independent of other macerals present in the kerogens. Thus, in a 50/50 mixture of liptinite and vitrinite, one would expect liptinite to decompose first to produce products characteristic of liptinite, and the vitrinite to decompose at a later stage due to ita high activation energy distribution and frequency factor. The activation energy distribution of the mixtures of these two kerogens was determined in two different ways. Figure 9a shows the activation energy distributions of three different mixtures of liptinite and vitrinite calculated from that of pure macerals by multiplying the activation energy distribution of the macerals by the corresponding
Kerogen Actiuation Energy Distribution Llptlnlte:Vltrlnlte
Energy & Fuels, Vol. 6, No. 6, 1992 803
Llpt1nlte:Vltrlnlte
"1
11
50
50
25
75
40
5 0 5 0
40
Actlvatlon Energy
1,
., .
25
., . .
75
~~
20 0
so
SO
70
80
Actlvatlon Energy
b)
(b)
Figure 9. (a) Calculated activation energy distributions of mixtures of liptinite and vitrinite. The activation energy distributionswere calculatedby multiplyingthe activationenergy distributionsof pure liptiniteand pure vitrinite by the appropriate proportion. (b)Measured activationenergy distributionsof pure liptinite,pure vitrinite, and mixtures of vitrinite and liptinite in different proportions. proportion. For example, the activation energy distribution of a 5050 mixture of liptinite and vitrinite in Figure 9a was calculated by multiplying the activation energy distribution of liptinite by 0.5 and that of vitrinite by 0.5 and displaying the net result as the activation energy distribution of the mixture. This leads to two distinct distributions centered around 52 and 69 kcal/mol. Figure 9b shows activation energy distributions of mixtures obtained by physically mixing liptinite and vitrinite in various proportions; pyrolyzing the mixtures at three widely differing heating rates, 40,4, and 0.4 K/min; and processing the resulting pyrolysis curves assuming severalparallel reactions with different activation energies and a common Arrhenius factor. The activation energy distribution obtained by physically mixing the macerals and pyrolyzing the mixture shown in Figure 9b does not exhibit the bimodal distribution of Figure 9a. This results from the assumption that the Arrhenius factor is the same for all activation energies. This confirms the hypothesis that the use of a single Arrhenius factor for all the activation energies leads to a false activation energy distribution. This in turn can lead to a false conclusion in maturetion models of sedimentary basins.
Future Research Previous studies have addressed the determination of the kerogen activation energy distribution for a single energy or for a normal distribution. A principal contribution of this study is the unbiased determination of the distribution subject to no assumptions as to the nature of the distribution. Two simplifyingassumptionswere made in order to render the mathematical determination of a general activation energy distribution feasible. Kerogen consists of a multitude of organic species. It is to be anticipated that the chemical pathways for the pyrolysis of kerogen are coupled, but the analysis in this paper assumes that the reactions proceed in parallel and that the reactions are decoupled. Until an extraordinarily compact chemical pathway scheme is postulated for kerogen pyrolysis, it seems necessary to preserve this simplifying assumption. The second critical assumption is that the Arrhenius preexponential factor A is independentof activation energy E. Early in the evolution of this study, we attempted to allow for some dependence, such as In A = u + bE; that is, the logarithm of preexponential is linear in activation energy. It quickly became apparent that this level of generalization had to be abandonedbecause of nonsensical answers. The effects of the problem ill-posedness and high correlation between preexponential and activation energy precluded determination of a general distribution when In A = a + bE. The temperature schedule employed here is strictly linear, ramping up. A strategically designed temperature schedule can better elucidate the chemistry of a set of parallel, first-order, irreversible kinetics. Superior temperature schedules would combine alternate heating, cooling, and intervals of constant temperature. This result is quantifiable via D-optimality by computing and comparing values of 1JTJI-1/2fordifferent experimental designs. This was done in this paper, for example, in the section Multiple Heating Rates. As mentioned earlier, deriving activation energies from a pyrolysis experiment is an illposed, inverse problem. The degree of ill-posedeness and extent of correlations among the parameters can be reduced by an improved experimental design such as mixed temperature schedules. The authors recognize that designing and implementing a laboratory procedure and apparatus that will accomplish mixed periods of heating and cooling interspersed with periods of constant temperature is a challenging prospect.
Acknowledgment. We acknowledge the work of Mr. M. Stauffer for technical assistance with the kerogen activation energy distribution determination. We also thank Dr. Bert Rust who generously supplied us with computer code for the determination of confidence intervals for the case of least-squaressubjectto nonnegativity constraints. We are also grateful to the reviewers for their advice and also for alerting us to a recently published paper by Braun, Burnham, Reynolds, and Clarkson (ref 29) on this topic as well as an overview series of three papers by Lerche (ref 11)and Liu and Lerche (refs 12 and 13) and to the Editor who directed our attention to the work of Essenhigh and Misra (ref 26).
