A Systematic Theoretical Kinetics Analysis for the Waddington Mechanism in the Low-Temperature Oxidation of Butene and Butanol Isomers

The Waddington mechanism, or the Waddington-type reaction pathway, is crucial for low-temperature oxidation of both alkenes and alcohols. In this study, the Waddington mechanism in the oxidation chemistry of butene and butanol isomers was systematically investigated. Fundamental quantum chemical calculations were conducted for the rate constants and thermodynamic properties of the reactions and species in this mechanism. Calculations were performed using two different ab initio solvers: Gaussian 09 and Orca 4.0.0, and two different kinetic solvers: PAPR and MultiWell, comprehensively. Temperature- and pressure-dependent rate constants were performed based on the transition state theory, associated with the Rice Ramsperger Kassel Marcus and master equation theories. Temperature-dependent thermochemistry (enthalpies of formation, entropy, and heat capacity) of all major species was also conducted, based on the statistical thermodynamics. Of the two types of reaction, dissociation reactions were significantly faster than isomerization reactions, while the rate constants of both reactions converged toward higher temperatures. In comparison, between two ab initio solvers, the barrier height difference among all isomerization and dissociation reactions was about 2 and 0.5 kcal/mol, respectively, resulting in less than 50%, and a factor of 2–10 differences for the predicted rate coefficients of the two reaction types, respectively. Comparing the two kinetic solvers, the rate constants of the isomerization reactions showed less than a 32% difference, while the rate of one dissociation reaction (P1 ↔ WDT12) exhibited 1–2 orders of magnitude discrepancy. Compared with results from the literature, both reaction rate coefficients (R4 and R5 reaction systems) and species’ thermochemistry (all closed shell molecules and open shell radicals R4 and R5) showed good agreement with the corresponding values obtained from the literature. All calculated results can be directly used for the chemical kinetic model development of butene and butanol isomer oxidation.


INTRODUCTION
The Waddington-type reaction pathway was proposed by Knox 1 in 1965 in "A new mechanism for the low temperature oxidation of hydrocarbons in the gas phase". It was subsequently tested experimentally by Ray and Waddington et al. 2−4 in 1971 and 1973 when they investigated the low-temperature (184−400°C ) oxidation of a series of alkenes, including ethylene, propene, 1-butene, 2-butene, isobutene, 2-methyl-2-butene, and 2,3dimethyl-2-butene, etc. It was formally named the Waddington mechanism by Wilk and Stark et al. 5,6 in 1989 and in 1995 in their study of the low-and intermediate-temperature chemistry of propene. This reaction pathway in alkene oxidation is explained in Figure 1, where R 1 −R 4 represent any alkyl group or H atom. The hydroxyl (ȮH) radical, added to an alkene, forms a β-alkylhydroxy radical, followed by the addition of molecular oxygen to form an alkyl-hydroxy-peroxy radical. This radical then undergoes a six-member ring isomerization by abstracting an H atom from the ȮH moiety to form an alkyl-hydroperoxide-alkoxy radical, followed by C−C and O−O bond fission, leading to the formation of formaldehyde or ketone plus the ȮH radical. In this reaction route, six-member ring isomerization and subsequent dissociation reactions (circled by dashed red lines) are treated as a typical Waddington mechanism. In this case, the β-alkylhydroxy radical, resulting from the ȮH radical addition to an alkene, can also result from the H atom abstraction from the βpositon of an alcohol molecule, making this reaction pathway important for both alkene and alcohol oxidation chemistry.
Review papers and a textbook published by Zador,7 Battin-Leclerc, 8,9 and Sarathy et al. 10 and more recent papers on the detailed kinetic model development of butanol isomers, 11 isobutene, 12 2-butene, 13 and 1-butene, 14 have also noted the importance of this reaction class in the oxidation chemistry of both olefins and alcohols. To predict the rate coefficient of the Waddington mechanism, Benson 15 proposed an empirical method to estimate the activation energy for an isomerization reaction from the sum of two contributions: (i) the activation energy for H atom abstraction from the molecule by analogous radicals and (ii) the strain energy involved in the cyclic transition state. More fundamental quantum calculations have been performed mainly for the C2−C4 hydrocarbon oxidation systems. 16−27 All these Waddington mechanism-related studies are summarized in Table S1 of the Supporting Information 1; however, there is still a lack of advanced theoretical treatment of the rate coefficient of this reaction class, and few studies have been concerned specifically with the effect of the reactivity of the isomeric molecule structures.
In light of these considerations, a systematic study of the Waddington mechanism in the low-temperature chemistry of butene and butanol isomers was conducted. In general, the alkyl-hydroxy-peroxy radicals were generated either from ȮH radical addition to butenes, followed by O 2 addition, or from β-positon H atom abstraction of butanols, followed by O 2 addition, as shown in Figure 2. Fundamental quantum chemical calculations were performed for the rate constants of the Waddington mechanism for the C4 reaction systems, as shown in Figure 3.

COMPUTATIONAL METHODS
Geometries, frequencies, and zero-point energies (ZPEs), as well as the hindered rotation treatments for lower frequency modes, were determined at the M06-2X 28 level of theory using the 6-311++G(d,p) 29,30 basis set. Vibrational frequencies and ZPEs were scaled by 0.983 and 0.9698, respectively. 28 Intrinsic reaction coordinate (IRC) calculations 31 for the optimized transition states were performed to verify that the given transition states were those expected, which connected the desired reactants and products. The optimized geometries of all species and TSs are summarized in Supporting Information 1; their vibrational frequencies can be found in the input files in Supporting Information 4.
Single-point energies (SPEs) were determined for the M06-2X geometries with the CCSD(T) 32 method with the cc-pVTZ and cc-pVQZ 29,33 basis sets, using the Gaussian 09 solver. 34 To obtain accurate energies at a relatively lower cost, the DLPNO− CCSD(T)/cc-pVTZ and DLPNO−CCSD(T)/cc-pVQZ 35 levels of theory, with the TightPNO truncation threshold, were also adopted using the Orca 4.0.0 solver. 36 The resulting SPEs were extrapolated to the complete basis set limit (CBS) level using the following formula (where method = CCSD(T) and DLPNO− CCSD(T)): 37,38 Temperature-and pressure-dependent rate constants were calculated based on the transition state theory (TST), associated with the Rice Ramsperger Kassel Marcus (RRKM) and master equation (ME) theories, 39 at pressures of 0.01−100 atm over the temperature range of 298.15−2000 K. The Lennard-Jones (L-J) potential 40 predicted the interaction between the reactant and N 2 bath gas. L-J parameters were empirically estimated as σ = 3.61 Å and ε = 68 cm −1 for N 2 and σ = 5.05 Å and ε = 1449 cm −1 for C 4 H 9 Ȯ3 radicals, using the OneDMin module packed in the PAPR solver. 41 The collisional energy transfer was represented by a single-parameter exponential-down model with ⟨ΔE down ⟩ = 200 × (T/300) 0.75 cm −1 . 42,42 The torsional modes of the −CH 3 group, −OH group, and −COOH groups were treated as 1-D hindered rotors, with hindrance potentials evaluated at the M06-2X/6-311++G(d,p) level of theory. Quantum mechanical tunneling was taken into account for an  The Journal of Physical Chemistry A pubs.acs.org/JPCA Article unsymmetrical Eckart barrier model, 43 and the calculated rate coefficients were fitted to a modified Arrhenius expression as In the quantum chemical method for the thermodynamic properties calculation, zero Kelvin energies (ZKEs) were obtained using combined compound methods: a combined compound method G3/G4/CBS-APNO, 44−46 which yielded results approaching chemical accuracy (arbitrarily, ≈4 kJ mol −1 ) when benchmarked against the enthalpy of formation from active thermochemical tables (ATcT). 47−49 Thereafter, the atomization method was utilized to derive the enthalpies of formation at 0 K. Basically, a molecule or radical was divided into its component atoms via the reaction in which the theoretical atomization energy at 0 K (TAE 0 ) can be calculated by where H 0 is the enthalpy of formation at 0 K calculated using each compound method. Thereafter, the enthalpy of formation of the species (Δ f H 0 ) was calculated knowing the TAE 0 and the standard formation enthalpies of the component atoms in their gaseous state from the ATcT, 4−6 shown in Table 1, via To clarify: All density function theory (DFT) calculations, CCSD(T) calculations, and G3/G4/CBS-APNO calculations were performed using Gaussian 09; 50 DLPNO−CCSD(T) calculations were conducted using the Orca 4.0.0. 36 Highpressure limit (HPL) rate coefficients were obtained using the MultiWell solver, 51,52 while temperature-and pressure-dependent rate coefficients were performed with the PAPR solver. 41 Thermodynamic parameters were obtained using the MultiWell solver and fitted to the NASA polynomial using the Fitdat utility in ANSYS CHEMKIN-PRO. 53 3. RESULTS AND DISCUSSION 3.1. C 4 H 9 Ȯ3 Potential Energy Surfaces and Barrier Heights. Corresponding to the reaction systems shown in Figure 3, the C 4 H 9 Ȯ3 potential energy surfaces (PESs) were conducted for 15 species and 10 transition states (TSs) in total, shown in Figure 4. Notably, the T1 diagnostic values for all the species and TSs on the PES were below 0.03, which proved the reliability of using the single reference treatment; detailed results have been summarized in Table S2 of Supporting Information 1. The temperature-and pressure-dependent rate constants of all isomerization, dissociation, and chemical activation reactions (15 reactions) and thermodynamic properties of all species (15    Table 2 shows the calculated barrier heights using two ab initio solvers (Gaussian 09 and Orca 4.0.0) for all isomerization and dissociation reactions. In addition, theoretical results were summarized from several studies in the literature: Chen et al., 20 Sun et al., 21 Welz et al., 24 and Lizardo-Huerta et al. 26 The independent theories adopted in the various studies resulted in different values of the barrier height. It should be noted that (a) results from Chen et al. and Sun et al. originated from the same research group but showed a discrepancy of about 5 kcal/mol; (b) Welz et al. used two different levels of theory achieving two widely disparate results for the barrier height of the R1 ↔ P1 reaction (about 8 kcal/mol). The T1 diagnostic for the SPE calculation of the transition state (R1P1-TS) for this reaction using the RQCISD(T) method was reported as 0.096, indicating that a multireference method was required for the treatment of electronic configurations. This study adopted an advanced level of theory (CCSD(T)/cc-pVTZ, QZ//M06-2X/6-311++G-(d,p)) in comparison to studies in the literature, which generally resulted in very consistent barrier height values: Gaussian predicted barriers were about 2 kcal/mol higher for all isomerization reactions other than Orca, and both solvers generated identical barrier height values for all dissociation reactions (less than 0.6 kcal/mol difference).
3.2. Rate Constant Comparisons. All calculated rate constants were fitted to the modified Arrhenius expressions in the corresponding calculated temperature range of each reaction, and every result is summarized in Supporting Information 2, which can be directly incorporated into a chemical kinetic model. All reactions can be compared in many different aspects/dimensions by comparison among two ab initio solvers, two kinetic solvers, various pressures, different PESs or reaction channels, different reaction types (isomerization, dissociation, and chemical activation), etc. The following comparisons were considered to be the most important: • HPL rate constants calculated by two kinetic solvers (MultiWell and PAPR) for two reaction types (isomerization and dissociation). • HPL rate constants systematically calculated by two ab initio solvers (Gaussian and Orca) for all five reactions systems (R1−R5 reaction systems), including comparison against results in the literature. Figure 5 shows kinetic solver comparisons for the HPL rate constants of two reaction types: (a) isomerization and (b) dissociation; different colors correspond to various reactions, and different line types represent results obtained using various solvers. Overall, the rate constants calculated by the PAPR solver were consistently faster than from the MultiWell solver. Agreement was very good for the five isomerization reactions, with less than a 32% difference for the entire temperature range. In all the dissociation reactions, all the rate constants converged toward higher temperatures and became less temperature dependent. Two solvers also showed reasonably good agreement (except the P1 ↔ WDT12 reaction); the difference was found to be 1−2 orders of magnitude, depending on temperature. The very first publication compared the HPL rate constant calculated using the Multiwell and PAPR solvers by Li et al. 54 In that paper, some discrepancies were also found between the two solvers. The authors are aware that PAPR is a recently developed solver, indicating that the solver may still require some refinements, compared to MultiWell, The discrepancies obtained in this work constitute useful information which can be communicated to the code developers. Note that all ab initio results were calculated using the Gaussian solver for comparison; the corresponding results for the Orca solver are summarized in Figure S2 in Supporting Information 1.

Ab Initio Solvers and Comparison of Literature
Results. Figure 6 shows comparison of the ab initio solvers for the HPL rate constants of two reaction types (isomerization and dissociation) in all five reaction systems; different colors and lines correspond to various reactions and results from different solvers, respectively, in each figure. For the R4 and R5 reaction systems, theoretical results from Sun et al. 21 and Lizardo-Huerta et al. 26 were also plotted using dotted and dashed dotted lines, respectively. First, the rate constants of the dissociation reactions were significantly faster than those of the isomerization reactions for each reaction system (or PES), by 1−16 orders of magnitude, depending on the reactions, temperatures, and solvers. However, such differences decreased or converged with the increasing temperature, indicating that the differences were in energies (barrier heights) while the pre-exponential factor (vibration frequencies and partition functions) were consistent. The Journal of Physical Chemistry A pubs.acs.org/JPCA Article Second, compared to the Gaussian solver, the Orca solver predicted higher and lower rate constants for the dissociation and isomerization reactions, respectively, due to the different barrier heights predicted by the two solvers (summarized in Table 1). Third, compared to the results in the literature, the rates of the dissociation reactions calculated in this study fell between those of Sun and Lizardo-Huerta, while they agreed well with the rates for the isomerization reactions in both literature studiesespecially for the R5 ↔ P5 reaction. Note that the MultiWell solver was used to generate all the rate constants in this section's comparison, and the corresponding results for the PAPR solver are summarized in Figure S3 in Supporting Information 1.

Thermochemistry Calculation.
In this section, three well-known thermochemistry databases were selected to validate the calculated thermochemical properties in this work:  Table 3 compares thermochemistry for all aldehyde and ketone molecules in this study: propanal (C 2 H 5 CHO), acetone (CH 3 COCH 3 ), acetaldehyde (CH 3 CHO), and formaldehyde (CH 2 O). Excellent agreement was obtained for the 298 K enthalpies of formation (Δ f H Θ ), 298 K entropies (S Θ ), and heat capacities (C p ) at selected temperatures (differences were within 0.5 kcal mol −1 and 0.5 cal K −1 mol −1 , respectively). This indicated the reliability of the method/approach adopted for computing thermochemical values with a much lower calculation load. Table 4 summarizes the thermochemical properties for all C 4 H 9 O 3 radicals calculated using the MultiWell solver. Thermochemistry results for R4 and R5 were obtained using the PAPR solver and the theoretical study of Sun et al.; 21 they are listed in the Table 3. Excellent agreement can be seen for the enthalpies of formation, entropies, and heat capacities of R4 and R5 of the three sources. All thermodynamic properties were fitted to the NASA polynomial format and are summarized in Supporting Information 2, which can be directly incorporated into the chemical kinetic model.

IMPLICATIONS FOR KINETIC MODEL DEVELOPMENT
One of the major applications for accurate rate constants and thermochemistry is in the development of detailed chemical kinetic models. In this section, the HPL rate coefficients and thermodynamic properties calculated were incorporated into the AramcoMech 2.0 model 13 using the Gaussian and MultiWell   53 under the condition of φ = 1.0, p = 20 atm, and T = 650− 1300 K. Figure 7 shows a comparison between the original model and the model incorporated with new rates. The IDT was affected mainly at low to intermediate temperatures (600−900 K), quantitatively; at 750 K, the IDT of isobutene was increased by a factor of about two. This emphasized the requirement for further investigation on kinetic model development or improvement, which is beyond the scope of this study.

CONCLUSIONS
This study provided a systematic theoretical treatment of the Waddington mechanism, associated with C4 hydrocarbons (butene isomers) and oxygenates (butanol isomers) oxidation. The temperature-and pressure-dependent rate coefficients for isomerization, decomposition, and chemical activation reactions were investigated using RRKM/ME analyses. The thermodynamic properties of all C 4 H 9 Ȯ3 radicals were also evaluated. Two ab initio solversGaussian and Orcawere employed to conduct the geometry optimization, vibrational frequency, rotational constant, and single-point energy calculation, etc. Moreover, two kinetic solvers, MultiWell and PAPR, were also adopted to calculate the rate constants and thermochemistry, based on the transition state theory (TST) and statistical thermodynamics, respectively. All calculated results were comprehensively cross-compared. Major conclusions are summarized as follows: • Two ab initio solvers predicted similar barrier heights for the five isomerization and five dissociation reactions, with a difference of about 2 and 0.5 kcal/mol, respectively. • Based on these differences, the HPL rate constants of the dissociation reactions showed little difference (less than 50%), while those of the isomerization reactions showed much greater discrepancies (a factor of 2−10). • Dissociation reactions were significantly faster than isomerization reactions, while the rate constants of both reactions converged toward higher temperatures. • Two kinetic solvers generated an overlapping HPL rate constant of isomerization reactions (less than 32% difference); however, the rate of the P1 ↔ WDT12 reaction exhibited discrepancies of 1−2 orders of magnitude. • Calculated HPL rate coefficients for the R4 and R5 reaction systems agreed well with results from the corresponding literature. • Thermochemical properties calculated in this study showed excellent agreement with results in the literature for both closed shell molecules and open shell radicals. From an application perspective, all calculated results in this study can be used to develop the chemical kinetic model of butene and butanol isomer oxidation and to estimate the rate rules for the Waddington reaction pathways of larger olefin and alcohol oxidation.
Studies related to the Waddington mechanism; kinetics and ab initio solvers; T1 diagnostics and geometries of all species and transitions states; conformer distributions (PDF) All fitted rate coefficients and thermodynamic properties (TXT) Species glossary (PDF) All input and output files for MultiWell and PAPR solvers (ZIP)