New Approach To Elucidate Compensation Effect between Kinetic

May 27, 2007 - The study of thermogravimetric curves is commonly used to elucidate the likely processes involved during pyrolysis as well as determine...
1 downloads 0 Views 398KB Size
4382

Ind. Eng. Chem. Res. 2007, 46, 4382-4389

New Approach To Elucidate Compensation Effect between Kinetic Parameters in Thermogravimetric Data Antonio Marcilla, Amparo Gomez, Sergio Menargues, and Juan Carlos Garcia-Quesada* Chemical Engineering Department, UniVersity of Alicante, 03690 South Vicente del Raspeig, Alicante, Spain

The study of thermogravimetric curves is commonly used to elucidate the likely processes involved during pyrolysis as well as determine the corresponding kinetic parameters. One of the main problems reported in the determination of such parameters is their ability to couple each other in such way that different sets of kinetic parameters can describe similar conversion degree curves once a kinetic model has been selected. Least-squared methods can be used to calculate the optimum kinetic parameters that minimize the squared differences between experimental and calculated data. Although the calculation can be successfully done, there are still questions to be answered. In the present work, a procedure to register the compensation effect between kinetic parameters is suggested. Introduction Kinetic analysis is a topic that has been frequently associated to pyrolysis for decades. Data obtained during the thermal decomposition of biomass, polymers, etc., can be used to achieve a semiquantitative analysis or alternatively kinetically treated in order to evaluate the pseudo kinetic parameters associated to the likely processes involved. Kinetic analysis of data obtained during pyrolysis or combustion of materials can be easily found in the literature published during the last few years, such as, for example, refs 1-8. Typically, data obtained from a thermobalance has been used in order to obtain the kinetic parameters associated to the single process reaction rate equation

dR ) Ae(-E/RT)f(R) dt

(1)

where R is the conversion degree, t is the time, A is the preexponential factor, E is the activation energy, T is the absolute temperature, and f(R) is a conversion degree function which depends on the reaction pattern and hence on the kinetic model considered. Although kinetic analysis has been widely utilized, application of eq 1 to calculate the kinetic parameters may be still somewhat controversial, being necessary to make some considerations about different aspects. (1) Is eq 1 really representative of the process under study? (2) How the kinetic analysis procedure employed. (3) The way the parameters are considered. (4) The compensation effect between kinetic parameters. (1) Is eq 1 really representatiVe of the process under study? A very important factor to bear in mind is the type of materials (polymeric materials, biomass, sludge, etc.) studied. Due to their nature, complex reaction patterns, far from a single or elementary step, may occur during their degradation. In this case, conversion degree curves can be used to obtain pseudo kinetic parameters using eq 1 as a simplified hypothesis but not for model validation nor elucidation of the likely processes involved.9 In this case, kinetic parameters obtained should not be asked to have physical meaning since their only role should be the * To whom correspondence should be addressed. Fax: 34 96 590 3826. E-mail: [email protected].

Figure 1. Experimental pyrolysis conversion degree curves obtained at 5, 10, and 35 K/min.

Figure 2. General aspect of the surfaces obtained in a (50% range of activation energy and preexponential factor (a constant reaction order has been employed: n ) 1) around the reference values.

mathematical correlation of experimental data under the conditions studied. (2) The kinetic analysis procedure employed. In the case that eq 1 is accepted to be representative, it is worth mentioning that there is not any standardized procedure to evaluate kinetic parameters. Although the world seems to experience a globalization trend, scientists do not agree, for the moment, about what is the best procedure to obtain kinetic parameters and loads

10.1021/ie061524w CCC: $37.00 © 2007 American Chemical Society Published on Web 05/27/2007

Ind. Eng. Chem. Res., Vol. 46, No. 13, 2007 4383

Figure 3. Calculations obtained by considering a fixed reaction order: (a) aspect of the error surface, (b) projection of the levels of error equal or below 4 × 10-5 on the E/R - ln A plane, and (c) magnification of Figure 3a close to the minimum of the error surface.

Figure 5. Calculations obtained by considering a fixed preexponential factor: (a) aspect of the error surface and (b) projection of the levels of error equal or below 4 × 10-5 on the E/R - n plane.

Figure 4. Calculations obtained by considering a fixed activation energy: (a) aspect of the error surface and (b) projection of the levels of error equal to or below 4 × 10-5 on the n - ln A plane.

of likely procedures can be continuously found in the literature in which each author claims the most suitable method. Never-

theless, the most worrying fact is that each procedure may yield somewhat different kinetic parameters, as stated from the ICTAC kinetic project results.10 The different methodology followed, the inherent error associated to each procedure, together with the compensation effect between kinetic param-

4384

Ind. Eng. Chem. Res., Vol. 46, No. 13, 2007

Figure 6. Aspect of different cutting planes for average and absolute error surfaces, (a and b) at ln A (s-1) ) 6.0 and (c and d) at E ) 60 kJ/mol, by considering a fixed reaction order and different number of heating rates: one single heating rate (0.1 K/min), two heating rates (0.1, 0.5 K/min), four heating rates (0.1, 0.5, 1, 2 K/min), and eight heating rates (0.1, 0.5, 1, 2, 5, 10, 15, 20 K/min).

eters could be the most likely reasons of disagreement between the results obtained by different researchers. (3) The way the kinetic parameters are considered. When speaking about kinetic parameters it is very important to have in mind that the whole kinetic triplet should be considered: activation energy, preexponential factor, and conversion degree function. Nevertheless, it is possible to find in the literature work where authors consider and compare just one of the kinetic parameters (normally the activation energy), although each isolated parameter is useless.11 (4) Compensation effect between kinetic parameters. The compensation effect between kinetic parameters has been treated extensively in the literature.12-18 Most authors have paid attention to the compensation effect between the activation energy and the preexponential factor once the conversion degree function has been selected. These parameters may adopt different values in order to represent the same conversion degree behavior, but it is important to have in mind that not necessarily

all of them have the same grade of accuracy. In general, the dependence reported between the logarithm of the preexponential factor and the activation energy is linear. The compensation effect between these kinetic parameters is so relevant that there are methods for kinetic parameter evaluation that are based on this fact, such as, for example, the invariant kinetic parameters method.19 However, as mentioned before, not only does the activation energy couple with the preexponential factor but there is another factor to have in mind: the conversion degree function considered. For example, if the reaction order model is considered, the characteristic parameter to be considered is the reaction order. Then, if a reaction process can be described by a triplet, the compensation effect of the whole triplet should be analyzed, which is the aim of the present work. The analysis has been applied to generated data as well as to an example of experimental data.

Ind. Eng. Chem. Res., Vol. 46, No. 13, 2007 4385

or deviation versus the activation energy-preexponential factor-reaction order. Due to the impossibility of the construction of 4D graphs, 3D charts can be used alternatively. In this case, three different alternative 3D charts can be used since they can represent the error or deviation versus two of the kinetic parameters (keeping the third as a constant). Those zones of the 3D charts, with an error below that certain value, should correspond to the set of kinetic parameters which compensate for each other in order to represent a certain conversion degree curve. This limit value for a good fit is logically arbitrary and depends on the personal grade of satisfaction. In order to illustrate the shape of the 3D charts, a reference set of kinetic parameters has been considered (ln Aref (s-1) ) 6.0, Eref ) 60.0 kJ/mol, nref ) 1) and a conversion degree curve has been obtained by integration of the reaction rate equation at a given heating rate (for example 0.1 K/min) Figure 7. Experimental and calculated data from optimized kinetic parameters.

Experimental Section A LDPE 780R supplied by Dow Chemical (MFI ) 7 g/10 min and density 923 kg/m3) has been employed. The thermogravimetric behavior of this polymer has been obtained in a Netzsch TG 209 thermobalance. Experiments for determination of the thermal decomposition temperature were carried out with an initial sample mass around 5 mg and at heating rates of 5, 10, and 35 K/min. All samples were in powder form. The samples were placed in open Al2O3 pans. The carrier gases were N2 (99.99% minimum purity) at a flow rate of around 30 mL STP/min. To ensure measurement of the actual sample temperature, a calibration of the temperature was performed using the Curie-point transition of standard metals, and a parabolic correction was applied. Results obtained are shown in Figure 1, where the conversion degree of the decomposition process, R, is represented versus the temperature. Basics of the Procedure As mentioned above, as a consequence of the compensation effect, the kinetic triplet may adopt different values in order to represent similar behaviors but not necessarily identical behaviors. For example, if a certain conversion degree curve is considered as a reference, the conversion degree curves generated using different combinations of the kinetic triplet may reach different grades of similarity with respect to the reference curve. Thus, a measure of such similarity needs to be defined. A possibility might be to evaluate a variation coefficient or alternatively an averaged error which describes the deviation between both curves, taking into account the number of data considered (N)

error )

(Rref - Rgen)2 ∑ all data N

(2)

Thus, all sets of the kinetic triplet that give an error below a given value, considered as representative of a good grade of similarity between the reference and the generated data, should describe the compensation effect of the kinetic triplet. In consequence, if eq 1 is used to represent the process and the reaction order is selected for f(R), the actual compensation effect among the three kinetic parameters (E, A, and n) should be represented in four-dimensional charts, representing the error

dRref ) Ae(-E/RT)(1 - Rref)n f Rref ) dt

∫0t Ae(-E/RT)(1 - Rref)n dt

(3)

Afterward, each kinetic parameter was modified in a relatively wide range

0.5 ln Aref e ln A e 1.5 ln Aref

(4)

0.5Eref e E e 1.5Eref

(5)

0.5nref e n e 1.5nref

(6)

As commented before, the compensation effect among kinetic parameters has to be evaluated by pairs and keeping the third constant. Thus, three different combinations can be obtained

type I:0.5 ln Aref e ln A e1.5 ln Aref, 0.5Eref e E e 1.5Erefnref ) constant type II:0.5 ln Aref e ln A e 1.5 ln Aref, 0.5nref e n e 1.5nref, Eref ) constant type III: 0.5nref e n e 1.5nref, 0.5Eref e E e 1.5Eref, ln Aref ) constant A mesh of around 180 × 180 points has been calculated within the corresponding intervals. For each set of kinetic parameters, the rate equation (eq 1) has been integrated by the fourth-order Runge-Kutta procedure (the number of intervals employed was 1000 in order to ensure that the integration result did not depend on that parameter), obtaining the corresponding generated conversion degree curves

dRgen ) Ae(-E/RT)(1 - Rgen)n f Rgen ) dt

∫0t Ae(-E/RT)(1 - Rgen)n dt

(7)

Then, the error has been calculated by means of the eq 2, allowing construction of Figures 2-5. All calculations have been done by means of a Matlab program. When the parameter interval studied is large, the surfaces look like Figure 2. Nevertheless, when the interval is restricted to the zone near the parameters used as reference, all the figures

4386

Ind. Eng. Chem. Res., Vol. 46, No. 13, 2007

Figure 8. Error surfaces obtained when considering a fixed reaction order (n ) 0.65): (a) overall surface, (b) magnification around the minimum of the surface, (c) projection of Figure 8b on the E/R - ln A plane (projection of the surface with an error lower than 3.0 × 10-4), and (d) projection of Figure 8b on the E/R - ln A plane (projection of the surface with an error lower than 2.0 × 10-4).

consist of surfaces resembling folded sheets. Inspection of such surfaces reveals that they are curved in all directions as the different projections prove. The three projections on the planes (all at constant error) look elliptical, showing the correlation between parameters studied. In the case of Figure 3, this ellipse turns almost into a straight line, proving the high degree of interrelation between these parameters. The eccentricity of the ellipses observed when studying the other two combinations of parameters is lower (Figures 4 and 5), showing a less defined correlation between the preexponential factor and the reaction order and the activation energy and the reaction order. The major axis of the corresponding ellipse approximately shows the path of the minimum in the observed valley. The direction of the minor axis marks the direction more sensitive to variation of the parameters, i.e., a more marked minimum. In the case of Figure 3, this direction is close to the diagonal, but in the case of the other two plots, this direction is more parallel to the reaction order axis in both cases, proving that the correlation is more sensitive to the preexponential factor or the activation energy rather to the reaction order. Due to the clear elliptical shape of the cutting planes at different error values in Figures 4 and 5, it is possible to find somewhat different kinetic parameters that offer the same error values. For example, in Figure 4, if an approximate error of 4 × 10-5 is considered, infinite combinations (lying on a ellipse) of preexponential factor and reaction order can be found; in this case, two examples among of the infinite likely combinations can be roughly extracted from this figure: n ) 0.85 - ln A (s-1) ) 3 and n ) 1.15 - ln A (s-1) ) 9, as the points s and t mark. Effect of the Heating Rate. Figures 2-5 have been drawn by considering data at one single heating rate, but use of

different heating rates may affect the shape of these surfaces. In order to illustrate the effect of the number of heating rates considered, another approach has been considered. In this case the error surfaces have been calculated generating data up to eight heating rates by considering a constant reaction order and preexponential factor-activation energies sweeps. The width of the surface, near the minimum, can be illustrated by obtaining cutting planes at constant values of ln A or E. Figure 6a shows such planes for ln A (s-1) ) 6 and E ) 60 kJ/mol. It is possible to observe how the utilization of different heating rates provokes a certain sharpening in the average error surfaces, indicating a certain decrease in the area of potential sets of kinetic parameters yielding a similar set of conversion degree curves. Nevertheless, the effect on the absolute error, i.e.

absolute error )

∑ (Rref - Rgem)2

(8)

all data

is more marked. As is possible to observe in Figure 6b,d, the absolute error (calculated by multiplying the average error by the number of points considered) reflects a relevant sharpening in the planes when increasing the number of heating rates and in consequence the number of data points. The same effect is observed for the average error surfaces, where a certain narrowing is observed in the constant activation energy cutting planes, Figure 6c. This fact (i.e., the sharpening of the absolute or average error surfaces) can be one of the reasons that from a optimization process point of view utilization of different heating rates has been considered more desirable11,19 than utilizing just one heating rate, since as the error or other alternative objective functions considered become sharper, it is easier to find the

Ind. Eng. Chem. Res., Vol. 46, No. 13, 2007 4387

Figure 9. Experimental and calculated data from arbitrary kinetic parameters which gives an error of 3.0 × 10-4. Kinetic parameters employed to generate calculated data were 1n A (s-1) ) 36.80, E ) 259.81 kJ/mol, and n ) 0.65.

solution to the problem, i.e., the minimum of such function. As expected, the minimum error was reached in both cases for the same set of kinetic parameters. The average error will be hereafter be employed in the present work, unless otherwise indicated. Application of the Procedure to Experimental Data In order to correlate experimental data corresponding to the degradation process of a polyethylene in a thermobalance, a pseudo kinetic analysis has been employed. Although conversion curves only reflect one single-step process, it is well known that polymer degradation may occur via very complex reactions patterns involving a high number of reactions. Nevertheless, it is very common when studying polymer degradation to assume that when observing apparent single processes their kinetics could be represented by kinetic equations corresponding to one single process. In this case, we point out that the kinetic parameters calculated are actually pseudo kinetic parameters in that they only represent an average global process. Due to this fact, i.e., they are actually pseudo kinetic parameters, some of them may lack physical meaning. The experimental data have been correlated by means of eq 1 using a least-squared differences procedure in order to calculate pseudo kinetic parameters. In this case, the error calculated has been minimized using an EXCEL worksheet and the tool SOLVER. As in the case of the Matlab program, integration of the rate equation has been done by considering a high number of intervals (1000) in order to ensure an accurate value. The kinetic parameters which minimize the average error given by eq 2 are the following

1n A (s-1) = 40.74,

E = 34 135 K-1, n ) 0.64 R

The minimum average error reached is 1.6 × 10-4. Experimental and optimized conversion degree curves are represented in Figure 7. In order to establish a starting point of where to center the analysis of average surfaces error, the previous set of pseudo kinetic parameters has been employed and a surface of the type of Figure 8a has been obtained. Its analysis reveals that the

Figure 10. Projections of the error surfaces with an error beneath 3.0 × 10-4: (a) considering a constant activation energy (E ) 283.80 kJ/mol) and (b) considering a kinetic constant of ln A (s-1) = 40.74.

surface, as expected, presents a clear minimum at approximately (Figure 8b) ln A(s-1) = 40.7, E/R = 34 123 K - 1. Logically the minimum value of the error surface corresponds to the best mathematical set, but there are other combinations of pseudo kinetic parameters that can also give good results. As mentioned above, all these sets are located on the surface, below a certain value of error. For example, considering an error of 3.0 × 10-4, a correlation of the experimental data of the type of the Figure 9 can be obtained using the kinetic triplet: ln A (s-1) ) 36.80, E/R ) 31 250 K-1, n ) 0.65. The quality of the fit is acceptable, keeping in mind that data at different heating rates have been considered. Consequently, all the pseudo kinetic parameters combinations located on the surface of Figure 8a beneath this error could be also considered to offer a good correlation of experimental data. The projections of the surface at this error represent how activation energy and kinetic constant compensate each other at this level of fit quality. As occurred with generated data, experimental data show a surface whose projections constitute sharp areas, which are actually ellipses with a high eccentricity. In Figure 8c part of the projection at that level of error is shown, where the high number of likely kinetic parameter combinations is patent. If a more conservative criterion is considered (i.e., a lower error value such as, for example, 2 × 10-4), the area size decreases but there is still a

4388

Ind. Eng. Chem. Res., Vol. 46, No. 13, 2007

Figure 11. Representation of the projection of the levels of error equal or below 3 × 10-4 of the surfaces generated by considering data at each heating rate and the projection of the surface obtained by the simultaneous consideration of all the curves: (a) on the E/R - ln A plane, (b) on the n - ln A plane, and (c) on the E/R - n plane.

clear compensating effect between activation energy and kinetic constant (Figure 8d). Concerning the compensation effect between the activation energy and the kinetic constant with the reaction order, it is necessary to consider the corresponding charts. As occurred with generated data, surfaces obtained are in general similar to Figure 8a except for the fact that the corresponding projections, shown in Figure 10, are not as sharp as occurred in Figure 8c and d. In this case, the elliptical aspect of the projections is more marked and a relatively broad range of likely kinetic parameters sets can also be obtained. In contrast to Figure 4b, the relatively narrow aspect of Figure 8a is due to the fact that experimental data have certain dispersion and the minimum reached by each surface at a different heating rate is located at a different set of parameters, as possible to observe in Figure 11. In this figure the projections of the surfaces obtained for data at each heating rate on the different planes are represented with the corresponding projection of the surface generated simultaneously taking into account the three different heating rates. In general terms, the projections of the surfaces at each heating rate are centered in a different position on those planes and the areas of the projection are somewhat different, revealing that the coupling effect among kinetic parameters is different for each experimental curve. Furthermore, the projection of the overall contribution, i.e., the curve generated by considering the three different heating rates, represents an average behavior and in a

certain way the intersection among the rest of projections, a fact that ensures that the representative kinetic parameters of the center of the ellipse correspond to the best set of kinetic parameters. Concluding, analysis of the surfaces has also revealed that there exist certain compensation effects between kinetic constant and activation energy with the reaction order, although they are not as marked nor as relevant as the compensation effect between kinetic constant and activation energy. Conclusions In the present work the compensation effect among kinetic parameters has been studied. The size of the surfaces under a certain error value could be employed as a measure of the compensation effect between kinetic parameters. The methodology followed could be used to develop a useful method to evaluate compensation effect between kinetic parameters and establish, considering a certain reference error, the likely sets of the kinetic triplet which could reproduce similar conversion degree curves. In addition, evaluation of the error surfaces suggested in the present work can also be employed to support another kinetic analysis procedures in order to check whether the variation coefficient or objective function has reached an absolute minimum. The procedure suggested has been successfully applied to generated and experimental results,

Ind. Eng. Chem. Res., Vol. 46, No. 13, 2007 4389

obtaining similar conclusions and results in both cases, and could be used to standardize the kinetic analysis. Acknowledgment We are thankful for the financial support provided by the Spanish Comisio´n de Investigacio´n Cientı´fica y Tecnolo´gica de la Secretarı´a de Estado de Educacio´n, Universidades, Investigacio´n y Desarrollo and the European Community (FEDER refunds) (CICYT CTQ2004-02187), and Generalitat Valenciana (project GRUPOS03/159, ACOM06/162). Literature Cited (1) Grønli, M.; Antal, M. J.; Varhegyi, G. A Round-Robin Study of Cellulose Pyrolysis Kinetics by Thermogravimetry. Ind. Eng. Chem. Res. 1999, 38, 2238. (2) Wu, R. M.; Lee, D. J.; Chang, C. Y.; Shie, J. L. Fitting TGA data of oil sludge pyrolysis and oxidation by applying a model free approximation of the Arrhenius parameters. J. Anal. Appl. Pyrol. 2006, 76, 132. (3) Capart, R.; Khezamia, L.; Burnhamb, A. K. Assessment of various kinetic models for the pyrolysis of a microgranular cellulose. Thermochim. Acta 2004, 417, 79. (4) Slovaka, V.; Susak, P. Pitch pyrolysis kinetics from single TG curve. J. Anal. Appl. Pyrol. 2004, 72, 249. (5) Go´mez, C. J.; Varhegyi, G.; Puigjaner, L. Slow Pyrolysis of Woody Residues and an Herbaceous Biomass Crop: A Kinetic Study. Ind. Eng. Chem. Res. 2005, 44, 6650. (6) Varhegyi, G.; Szabo, P.; Jakab, E.; Till, F. Least squares criteria for the kinetic evaluation of thermoanalytical experiments. Examples from a char reactivity study. J. Anal. Appl. Pyrol. 2001, 57, 203. (7) Ranzi, E. A Wide-Range Kinetic Modeling Study of Oxidation and Combustion of Transportation Fuels and Surrogate Mixtures. Energy Fuels 2006, 20, 1024. (8) Grønli, M. G.; Varhegyi, G.; Di Blasi, C. Thermogravimetric Analysis and Devolatilization Kinetics of Wood. Ind. Eng. Chem. Res. 2002, 41, 4201. (9) Marcilla, A.; Garcia-Quesada, J. C.; Ruiz-Femenia, R. Additional considerations to the paper entitled Computational aspects of kinetic analysis. Part B: The ICTAC Kinetics Project-the decomposition kinetics of calcium

carbonate revisited, or some tips on survival in the kinetic minefield. Thermochim. Acta 2006, 445, 92. (10) Brown, M. E.; Maciejewski, M.; Vyazovkin, S.; Nomen, R.; Sempere, J.; Burnham, A.; Opfermann, J.; Strey, R.; Anderson, H. L.; Kemmler, A.; Keuleers, R.; Janssens, J.; Desseyn, H. O.; Li, C.-R.; Tang, Tong B.; Roduit, B.; Malek, J.; Mitsuhashi, T. Computational aspects of kinetic analysis Part A: The ICTAC kinetics project-data, methods and results. Thermochim. Acta 2000, 355, 125. (11) Maciejewski, M. Computational aspects of kinetic analysis. Part B: The ICTAC Kinetics Project-the decomposition kinetics of calcium carbonate revisited, or some tips on survival in the kinetic minefield. Thermochim. Acta. 2000, 355, 145. (12) Sun, D.; Wicker, L. Kinetic Compensation and the Role of Cations in Pectinesterase Catalysis. J. Agric. Food Chem. 1999, 47, 1471. (13) Poco, J. G. R.; Furlan, H.; Giudici, R. A. Discussion on Kinetic Compensation Effect and Anisotropy. J. Phys. Chem. B 2002, 106, 4873. (14) Zhang, J.; Xu, H.; Li, W. Kinetic study of NH3 decomposition over Ni nanoparticles: The role of La promoter, structure sensitivity and compensation effect. Appl. Catal. A: Gen. 2005, 296, 257. (15) Brown, M. E.; Galwey, A. K. The significance of ‘‘compensation effects” appearing in data published in ‘‘computational aspects of kinetic analysis: ICTAC project, 2000. Thermochim. Acta 2002, 387, 173. (16) Liu, N.; Zong, R.; Shu, L.; Zhou, J.; Fan, W. Kinetic Compensation Effect in Thermal Decomposition of Cellulosic Materials in Air Atmosphere. J. Appl. Polym. Sci. 89 2006, 135. (17) Ramis, X.; Salla, J. M.; Mas, C.; Manteco´n, A.; Serra, A. Kinetic Study by FTIR, TMA, and DSC of the Curing of a Mixture of DGEBA Resin and γ-Butyrolactone Catalyzed by Ytterbium Triflate. Appl. Polym. Sci. 2004, 92 (1), 381. (18) Galwey, A. K. Perennial problems and promising prospects in the kinetic analysis of nonisothermal rate data. Thermochim. Acta 2003, 407, 93. (19) Budrugeac, P.; Segal, E. On the Nonlinear Isoconversional Procedures to Evaluate the Activation Energy of Nonisothermal Reactions in Solids. E. Int. J. Chem. Kinet. 2001, 33 (10), 564.

ReceiVed for reView November 28, 2006 ReVised manuscript receiVed February 8, 2007 Accepted April 24, 2007 IE061524W