Analysis of acetic acid gas phase reactivity: Rate constant estimation and kinetic simulations

The gas phase reactivity of acetic acid was investigated combining first principle calculations with kinetic simulations. Rate constants for the unimolecular decomposition of acetic acid were determined integrating the 1D master equation over a Potential Energy Surface (PES) investigated at the M06-2X/aug-cc-pVTZ level. Energies were computed at the CCSD(T)/aug-cc-pVTZ level using a basis set size correction factor determined at the DF-MP2/aug-cc-pVQZ level. Three decomposition channels were considered: CO 2 + CH 4 , CH 2 CO + H 2 O, and CH 3 + COOH. Rate constants were computed in the 700–2100 K and 0.1–100 atm temperature and pressure ranges. The simulations show that the reaction is in fall off above 1200 K at pressures smaller than 10 atm. Successively, the PESs for acetic acid H-abstraction by H, OH, OOH, O 2 , and CH 3 were investigated at the same level of theory. Rate constants were computed accounting explicitly for the formation of entrance and exit van der Waals wells and their collisional stabilization. Energy barriers were determined at the CASPT2 level for H-abstraction by OH of the acidic H, since it has a strong mul- tireference character. The calculated rate constant is in good agreement with experiments and supports the experimental finding that at low temperatures it is pressure dependent. The calculated rate constants were used to update the POLIMI kinetic model and to simulate the pyrolysis and combustion of acetic acid. It was found that acetic acid decomposition and the formation of its direct decomposition products can be reasonably predicted. The formation of secondary products, such as H 2 and C 2 hydrocarbons, is underpredicted. This suggests that reaction routes not incorporated in the model may be active. Some hypotheses are formulated on which these may be.


Introduction
Acetic acid is found in significant concentrations in the troposphere and it is known to affect considerably the environmental chemistry. Its origin is both biogenic and anthropogenic. Biomass combustion and vehicle emissions are believed to be among the most important sources that can be directly related to human activities. Moreover, oxygenated species are very abundant in the tar released from biomass pyrolysis. Acetic acid is the major acidic component of bio-oils derived from biomass fast pyrolysis [1] . These renewable fuels are acidic liquids (pH 2-3), largely differing from petroleum fuels with regard to both their physical properties and chemical composition. The importance of acetic acid on the reactivity of pyrolysis products of cellulose was recently discussed by Fukutome et al. [2] . The modeling of each individual constituent of the bio-oils is not feasible, therefore these complex mixtures are characterized in terms of a limited number of representative chemical components (surrogate mixture). Acetic acid is the reference species used to represent carboxylic acids in bio-oil surrogates [3] . Therefore, a proper knowledge about its pyrolysis and combustion behavior is necessary to characterize the combustion of biomass pyrolysis oils. To control the formation and the successive reactivity of acetic acid it is necessary to understand its formation and decomposition pathways. In the present work we focus on the decomposition pathways of acetic acid and therefore we limit the literature analysis to this aspect.
Acetic acid can either decompose directly through a unimolecular process or through Habstraction reactions. The unimolecular decomposition of acetic acid has been the subject of both theoretical and experimental studies [4][5][6][7][8][9][10][11][12][13] . The first kinetic study was performed in 1949 by Bamford and Dewar using a flow method in a quartz tube in the 773-1173 K temperature range [4] . They found that the main decomposition reactions are: CH 3 COOH → CH 4 + CO 2 (R1) CH 3 COOH → CH 2 CO + H 2 O (R2) Successively Blake and Jackson investigated the same reaction at atmospheric pressure in the 733-868 K temperature range [5] . In both cases rate constants were determined fitting the evolution profiles of the products, thus not accounting for secondary chemistry. The first study in combustion relevant conditions (1300-1950 K) was performed by Mackie and Doolan using a single pulse shock tube [6] . The experimental results were interpreted using a kinetic model, which was used to fit R1 and R2 rate constants. A similar branching of about 0.5 was found between the two reaction channels. Two successive experimental studies by Saito et al. [7] and Butkoskaya et al. [8] found k R2 / k R1 branching ratios of 1 and 2, respectively. Recently, Elwardany et al. [9] studied the same system using a shock tube in the 1230-1821 K temperature range. The k R2 / k R1 branching ratio was about 2, while rate constants were about one order of magnitude faster than those of Mackie et al. Several theoretical stud-ies have investigated the decomposition mechanism of acetic acid. The potential energy surface (PES) was first studied by Nguyen et al. [10] and Duan and Page [11] , and more recently by Clark et al. [12] and Sun et al. [13] . All these studies agree on the PES qualitatively, though there is some difference in the k R2 / k R1 branching , predicted to be about 2 in the earlier studies and about 1 by Clark et al.
The second mechanism of decomposition of acetic acid is through bimolecular H-abstraction reactions. While H-abstraction from OH has been the subject of many studies in the literature, both experimental [14][15][16][17][18] and theoretical [17][18][19] , only Mendes et al. [19] have systematically investigated H-abstraction by other radicals, such as CH 3 , H, OOH, and O, which can significantly contribute to acetic acid reactivity. H-abstraction by O 2 has never been theoretically or experimentally investigated. In addition, most studies on abstraction reactions by OH were performed at room temperature. There is thus still some uncertainty concerning these reaction channels, and in particular for the complicated reaction between OH and acetic acid, which has evidenced some peculiar pressure dependence at room temperature [16] .
The purpose of the present work is to reexamine acetic acid reactivity combining first principles rate constant estimations and kinetic simulations. First, pressure dependent rate constants for the unimolecular decomposition of acetic acid were determined solving the multi-well master equation on a PES determined from high level ab initio calculations. Then rate constants for H-abstraction by H, OH, OOH, CH 3 , and O 2 were determined. All electronic structure calculations and rate estimations were performed at a level of theory that is higher than that presently reported in the literature. The capability of the updated acetic acid module of the POLIMI kinetic mechanism to predict pyrolysis, laminar burning velocities, and speciation in low pressure flames is then discussed.

Theoretical calculations
The unimolecular decomposition rate of acetic acid was determined performing master equation (ME) simulations on an ab initio PES. The PES was initially investigated at the M06-2X/6-311 + G(d,p) level and then refined at the M06-2X/aug-cc-pVTZ level, used also to compute the Hessian matrix. Energies were determined at the CCSD(T)/augcc-pVTZ level and corrected with the difference between Density-Fitting (DF) MP2 energies computed using aug-cc-pVQZ and aug-cc-pVTZ basis sets [20] . Torsional motions were treated in the 1D hindered rotor approximation using a M06-2X/6-311 + G(d,p) rotational PES. Torsional vibrational frequencies were projected out from the Hessian as suggested by Sharma et al. [21] , so that no human intervention is necessary. Quantum tunneling was accounted for using the Eckart model. In addition to the reaction channels reported in the literature two homolytic decomposition routes were considered: H loss from the methyl group and decomposition to methyl and COOH. The transition state density of states (DOS) for all the reaction channels was determined using variational transition state theory. For reactions proceeding through saddle points variational effects were accounted for determining the DOS along the intrinsic reaction coordinate pathway. The homolytic decomposition reaction energies and vibrational frequencies along the reaction pathway were computed at the multireference CASPT2 level with the breaking bond lengths spaced by 0.2 Å . Calculations were performed using a (4e,4o) active space that includes the CH 3 -COOH σ and σ * orbitals and the π and π * orbitals of the carboxyl group. The resulting energy profile was rescaled to match the reaction energy computed at the CCSD(T) level. ME calculations were performed at 0.1, 1, 10, and 100 atm in the 700-2100 K temperature range. The collisional energy transfer was described using a single exponential energy transfer model with a temperature dependent E down of 260 ×(T/300) 0.85 cm −1 . The calculated rate constants were fitted to the modified Arrhenius form for each considered pressure. The same approach, made exception for CASPT2 calculations, was used to investigate the decomposition kinetics of CH 2 COOH, which can be produced either from the decomposition of acetic acid or, more easily, through H abstraction from the acetic acid methyl group. Computational details are reported as supplemental information.
The rate constants for H-abstraction reactions by OH, H, CH 3 , OOH, and O 2 were computed solving a two wells ME whose stationary points were the reactant and product van der Waals (vdW) wells and the transition state by which they are connected. Entrance and exit channels from the PES were computed using phase space theory. The advantage of this approach is that it allows to set properly the energy of the reactants as a lower bound for the reacting flux and contextually to account for the impact of the depth of the vdW wells on tunneling. Contributions of reactive fluxes from vdW wells that get collisionally stabilized can as well be computed properly, as their formation rate and reactivity is included in the ME solution. In addition, the rate of formation of the vdW wells sets un upper limit to the rate constant. For consistency with the study of the unimolecular reaction channels, geometries, frequencies, and energies of all stationary point were determined at the M06-2X/aug-cc-pVTZ level, while energies were computed at the CCSD(T) level using the protocol described above. Since both reactants and saddle points can have a considerable number of internal torsional degrees of freedom (four in the case of abstraction from OOH), the search for the minimum energy configuration of the well and the saddle point was performed starting from an initial guess geometry, specified in z-matrix format, whose torsional angles were then modified using a random value in the 0-360 range to generate uncorrelated guess structures. Up to 15 independent optimizations were performed to search for the minimum energy saddle point structure. A 2D hindered rotor model was used to describe H abstraction by OH from the methyl group [22] . In one case, for acidic H abstraction from methyl by OH, energies were computed at the CASPT2 level, as CCSD(T) T1 diagnostics evidenced a significant multireference character. Details of the used active space are reported in the results section.
All M06-2X calculations were performed using the G09 suite of codes [23] , while CCSD(T), DF-MP2, and CASPT2 calculations were performed using Molpro 2010 [20] . ME input files were generated using a new code, EStokTP, designed to perform automatically the investigation of the torsional conformation space, to project torsional motions from Hessians, and determine 1-D PES for rotors [24] . ME simulations were performed using MESS [22] . A summary of the adopted computational procedures is reported in Table S1.

Development of acetic acid kinetic model and kinetic simulations
The rate constants computed in this study were used to update the acetic acid sub-mechanism in the POLIMI kinetic model. Moving toward a unification of core mechanisms, the POLIMI mechanism was recently revised coupling the H 2 /O 2 and C1/C2 from Metcalfe et al. [25] , C3 from Burke et al. [26] , and heavier fuels from Ranzi et al. [27] . The thermochemical properties were adopted, when available, from the ATcT database of Ruscic or from Burcat's database [28] . The kinetic mechanism, with thermodynamic and transport properties, is available as Supplemental Material (SM) to this paper. All the kinetic model simulations presented in this work have been carried out using the solvers for ideal reactors and laminar 1D flames implemented in the OpenSMOKE ++ software [29] . Mixtureaverage diffusion coefficients were used in the flame simulations, including also thermal diffusion for the species transport equations. A large number of grid points were adopted to ensure grid-insensitive results (above 300).

Acetic acid unimolecular decomposition: temperature and pressure dependence
The PES for acetic acid decomposition is shown in Fig. 1 . It consists of two wells, CH 3 COOH and CH 2 C(OH) 2 , connected by one transition state, and of four decomposition channels.
The calculated energy barriers are in good agreement with those determined by Clark et al. [12] at the M06-2X/6-311 + G(2df,p) level, from which they differ by less than 1 kcal/mol, made exception for isomerization to CH 2 C(OH) 2 , for which our barrier is 1.4 kcal/mol higher. The difference is more significant with respect to Duan and Page [11] , who calculated, using the G2 method on CASSCF structures, a similar energy barrier of 71.8 kcal/mol for decomposition to CH 4 + CO 2 , but a much higher energy barrier for decomposition to CH 2 CO + H 2 O both from W1 (76.4 kcal/mol) and W2 (72.6 kcal/mol). It is likely that the difference with our data is due to the relatively low level of theory at which structures were determined, which does not account for dynamic correlation. Homolytic decomposition to CH 3 and COOH, included in the PES for the first time in this work, has a reaction energy of 91.8 kcal/mol, in good agreement with the ATcT value of 92.5 kcal/mol [28] .
Rate constants calculated at 1 atm assuming that W2 concentration is at a pseudo steady state with respect to W1 are compared with literature values in Fig. 2 for decomposition to CH 4 + CO 2 and in Fig. S1 for decomposition to CH 2 CO + H 2 O. Branching ratios to CH 4 + CO 2 and CH 3 + COOH are reported in Fig. 3 . Fig. 2. Comparison between rate constants for acetic acid decomposition to CH 4 + CO 2 calculated in this work and literature values [6,9,12,13] . Fig. 3. Comparison between branching ratios ( k CH 4 + CO 2 / k total) for acetic acid decomposition to CH 4 + CO 2 calculated in this work and literature values (both experimental and theoretical data determined at 1 atm) [6,9,12,13] .
The rate constants values here calculated are intermediate between the two experimental measurements of Mackie and Doolan [6] and Elwardany et al. [9] , though they are more similar to the former. It should be noted that in both experimental works the rates were fitted using a kinetic mechanism, so that uncertainties in the mechanism are reflected in the estimates. The agreement with the high pressure predictions of Clark et al. is quite good for both the CH 4 + CO 2 and the CH 2 CO + H 2 O reaction channels. On the other side, the agreement with the recent calculations of Sun et al. [13] is good only for the CH 4 + CO 2 channel, while the rates calculated by Sun et al . for the CH 3 + COOH and the CH 2 CO + H 2 O channels are about a factor of 2.5 higher and 3 slower than our rates, respectively.
The difference for what concerns the CH 3 + COOH channel is most likely determined by the fact that Sun et al. used phase space theory to determine rate constants of barrierless channel, which is a level of theory lower than the one used in the present calculations and it is effective mostly in providing higher limits for rate constants, rather than accurate estimates.
The comparison between literature and calculated branching ratios ( Fig. 3 )  The predicted pressure dependence of the rate of decomposition to CH 4 + CO 2 is reported in Fig. S2 and shows that fall off effects start being significant above 1200 K. Rate constants calculated at different pressures and interpolated in the modified Arrhenius form are reported in Tables S2 and  S3, while the rate constants and PES calculated for the unimolecular decomposition of CH 2 COOH are reported in Tables S4 and S5 and Fig. S3. On the basis of the level of theory used to calculate the rate constants and the comparison with literature data we estimate that the uncertainty of our rate constants is about factor of 2.

Abstraction rates
Rate constants for H-abstractions by OH, H, CH 3 , OOH, and 3 O 2 from the methyl and carbonyl groups of acetic acid were computed as described in the method section. The calculated rates are generally within a factor of 2 of those determined by Mendes et al. [19] (see Fig. S5a-i for a direct comparison). The difference is most likely determined by the higher level of theory used in the present work: ME coupled to variational vs conventional transition state theory to determine the rates, energy barriers estimated at the CCSD(T) level using large basis sets, and hindered rotor treatment. We estimate that the uncertainty of our rate constants is a factor of 2 or better. Two main deviations are found with respect to Mendes et al. estimates. The first is for OOH abstraction from the carboxyl group, where our rate is smaller than Mendes et al. by a factor increasing from 2 to 4 with temperature. A factor of 1.5 is determined by variational effects, not accounted for by Mendes. The remaining discrepancy, as the computed energy barriers are similar (23.3 kcal/mol in the present work vs 24.9 kcal/mol, which explains the better agreement at lower temperatures), is most likely due to the estimation of the hindered rotors DOS and vibrational frequencies. This is a significantly endothermic reaction (25.7 kcal/mol) with a submerged barrier and a product like transition state. Because of that, several low vibrational frequencies are found at the saddle point, indicating that the transition state is rather loose. The manual removal of the hindered rotor rotational frequency from the transition state frequencies may be complicated in such cases as it may be mixed with other internal motions. The projection technique we use in this work removes such ambiguity.
The second reaction for which a considerable difference with Mendes et al. is found is Habstraction by OH of the acidic H. Though this reaction has been the subject of several theoretical investigation [17][18][19] , it is interesting to notice that it was never explicitly reported that it has a strong multireference character. The CCSD(T) T1diagnostic we computed for the saddle point for the radical hydrogen abstraction (TSA) is quite high, 0.142. In addition, we found that, similarly to what found for formic acid [30] , the carboxyl H atom can be abstracted through a second mechanism, which involves a proton coupled electron transfer mechanism. We were able to find also this second saddle point (TSB), which has as well a high T1 diagnostic of 0.046. The two saddle point structures are compared in the SM. To evaluate reliably energy barriers we performed multireference calculations at the CASPT2 level for this specific reaction channel. The multireference nature of this reaction is determined both by the electronic degeneracy of the p orbitals of the OH reactants and by two resonance states of the CH 3 COO product. To account properly for the change in the electronic structure that takes place going from reactants to products the minimum active space must be relatively large. The active space we used consists of 15 electrons in 12 orbitals and includes: the lone pair and the radical center of the OH radical, the σ and σ * and the π and π * orbitals of the CH 3 C = OOH double bond, the σ and σ * orbitals of the CH 3 COO −H bond, the σ and σ * orbitals of CH 3 CO −OH, and the two lone pairs of the acetic acid oxygens. As the active space is rather large, it was not possible to determine the saddle point structure at this level of theory. Therefore M06-2X/aug-cc-pVTZ structures were used in the calculations. The energy barriers, computed at the CASPT2 level on DFT geometries with respect to the reactants kept at 10 Å using an IPEA shift of 0.25, are 2.21 and 0.71 kcal/mol for TSA and TSB, respectively. The calculated total rate constant for H-abstraction by OH from CH 3 COOH is compared with experimental and theoretical estimates in Fig. 4 , while saddle point geometries and the reaction PES are reported in Fig. S4.  Fig. 4. Comparison between experimental [14,15] , points and calculated [19, this work, 1 atm], lines, rate constants for the total H-abstraction from CH 3 COOH. The well skipping label means that the contribution to the rate constant from the collisionally stabilized van der Waals well was neglected.
Rate constants were computed both including and neglecting (well skipping) the contribution from the collisionally stabilized vdW reactant well summing TSA and TSB reaction fluxes. It was thus found that the inclusion of reacting fluxes from the vdW well is of key relevance in order to describe properly the low temperature dependence of the rate constant and it supports the claim that the rate is pressure dependent at low temperatures [16] . The rapid increase of the rate constant above 500 K is determined both by reacting fluxes through TSA becoming faster than those through TSB and by the increase of relevance of H-abstraction from the methyl group. The dominant reaction channel is Habstraction of the acidic hydrogen up to 850 K, above H-abstraction from methyl becomes slightly dominant.
Rate Constants interpolated in the modified Arrhenius form, structures, energy barriers, hindered rotor PES, vibrational frequencies of all stationary points are reported in Table S6 for each investigated reaction channel.

Kinetic simulations
The POLIMI kinetic model, updated with the acetic acid rate constants here determined, was used to simulate the shock tube experimental data of Mackie et al. for the thermal decomposition of acetic acid in the 1300-1950 K temperature range. The experimental data are compared with those calculated using the present mechanism, the Aramco mechanism [25,26] , and the mechanism proposed by Christensen and Konnov [31] in Fig. 5 .
The agreement is good for what concerns the prediction of acetic acid decomposition and the formation of CO 2 , while CH 4 composition profiles are slightly underestimated, CH 2 CO is overestimated, and CO, H 2 , as well as C 2 species profiles are underestimated. With respect to the literature, the present model performs generally better than the Christensen et al. model and, for what concerns CO, H 2 , and C 2 species slightly worse than the Aramco mechanism. The different performance of the Aramco model for the mentioned species is mostly due to its inclusion of a fast acetic acid decomposition channel to CO + OH + CH 3 . At 1600 K the rate of this channel is comparable to that of reactions R1 and R2. This is in contrast with our estimate of the rate of decomposition of acetic acid to CH 3 + COOH, which at 1600 K is a factor of 3 slower than the rate used in Aramco and the rate of the other decomposition channels. The discrepancy is most likely due to the fact that Aramco uses high pressure rate constants, while our ME simulations show that at 1600 K and 1 atm the homolytic decomposition channel is in significant fall off.
The reaction flux analysis, reported in Fig. S6, shows that at 1600 K several decomposition pathways are active for acetic acid: 36% decomposes to H 2 O and CH 2 CO, and a similar amount to CH 4 and CO 2 , which are thus primary decomposition products. The remaining decomposition pathways are H-abstraction (20%) and homolytic decomposition to CH 3 and COOH. CO is formed through several channels, mostly involving the formation of CH 2 CO as an intermediate species. The CH 2 COOH radical plays an important role in controlling the system reactivity in the present kinetic model, since recombining with H it can either terminate the radical chain forming CH 3 COOH, or propagate decomposing to CH 3 and COOH or OH and CH 2 CO. The dominant CH 2 COOH reaction channel is decomposition to CH 3 and COOH, which requires overcoming a barrier of 41.9 kcal/mol to isomerize to the reactive CH 3 COO intermediate. The bimolecular reactivity of CH 2 COOH is at present mostly unknown, made exception for the H addition reaction, which is the backward process of CH 3 COOH decomposition and for which it was then possible to determine channel specific temperature and pressure dependent rate constants solving the ME on the CH 3 COOH PES. The calculated rates are reported in Table S2. Additional channels involving CH 2 COOH that may play an important role in controlling the system reactivity are the addition of CH 3 , OH, and OOH to CH 2 COOH. In particular, CH 3 addition followed by C 2 H 5 + COOH decomposition may enhance the rate of formation of C2 species and decrease the CH 2 CO concentration, thus providing a possible explanation of the reason why these species are overestimated and underestimated, respectively, in the present kinetic simulations. In addition we note that CH 2 CO py- Fig. 5. Comparison between experimental data (symbols [6] ) and model simulations (lines) for the shock tube thermal decomposition of 5% acetic acid in Ar, at 1 atm and τ -0.1 ms. Red lines represent fuel conversion, black lines represent intermediates and products. Fig. 6. Comparison between experimental data (symbols [33] ) and model simulations (lines) for acetic acid combustion in laminar premixed flames. rolysis and combustion are well predicted by the present model (Fig. S7). Also, the simulations of the CH 3 COOH pyrolysis study of Sebbar et al. [32] show that at lower temperatures (1023 K) there is good agreement between model and experiment for what concerns acetic acid decomposition and ketene formation (Fig. S8).
Laminar flame speeds predicted by the POLIMI model, compared with experiments [32] in Fig. S9, are on the average overestimated by about 4 cm/s. This is a significant improvement with respect to the original POLIMI model, which overestimated the experimental data by 9-10 cm/s, and with respect to the POLIMI model using the acetic acid decomposition kinetics implemented in the Aramco model (Fig. S10). The sensitivity analysis reported in the SM shows this is independent on acetic acid decomposition kinetics.
Finally, the POLIMI updated model was used to simulate the combustion of acetic acid. Experimental data and simulations are compared in Fig. 6 . As found for acetic acid pyrolysis, the rate of decomposition of acetic acid is well predicted by the model, as well as CO 2 , H 2 O, and CO profiles. Ketene is underpredicted, while methane is overpredicted. These deviations are not consistent with what found in pyrolysis conditions. Additional simulation results are reported in Fig. S11, while the results of the simulation of the experimental results of Elwardany et al. [9] are reported in Fig. S12.

Conclusions
Rate constants for the decomposition of acetic acid through unimolecular and H-abstraction reac-tions were investigated theoretically. The calculated elementary reaction rates are in good agreement with experimental data. Simulation of acetic acid pyrolysis and combustion showed that the adopted kinetic model is able to predict well acetic acid decomposition and the formation of its main primary products in all the investigated systems. The only exception is ketene, whose formation is well predicted at low temperatures, but underpredicted in combustion and slightly overpredicted in pyrolysis. In addition, the formation of secondary species, in particular C 2 species and CO and H 2 in pyrolysis, is underpredicted. The system reactivity, in particular the flame speed, is sensitive to the kinetics of the CH 2 COOH reaction intermediate, which is formed after H-abstraction from the methyl group of acetic acid, and to ketene kinetics. Bimolecular reactions involving CH 2 COOH and other radicals that are present in high concentration, such as CH 3 and OH, may play an important role in controlling the system reactivity and influence the flame speed, which is currently overpredicted.
Further experimental studies would also be beneficial, since the two most recent literature experimental CH 3 COOH pyrolysis studies cannot be modeled using the same set of decomposition rates, which must be modified by a factor of 3 or more to fit the data.