44 Chapter 3 Secondary Organic Aerosol Formation by Heterogeneous Reactions of Aldehydes and Ketones: A Quantum Mechanical Study* * This chapter is reproduced by permission from “Secondary Organic Aerosol Formation by Heterogeneous Reactions of Aldehydes and Ketones: A Quantum Mechanical Study” by C. Tong, M. Blanco, W. A. Goddard III, and J. H. Seinfeld, Environmental Science and Technology, 40, 2323-2338, 2006. © 2006 American Chemical Society 45 3.1 Abstract Experimental studies have provided convincing evidence that aerosol-phase heterogeneous chemical reactions (possibly acid-catalyzed) are involved to some extent in the formation of secondary organic aerosol (SOA). We present a stepwise procedure to determine physical properties such as heats of formation, standard entropies, Gibbs free energies of formation, and solvation energies from quantum mechanics (QM), for various short-chain aldehydes and ketones. We show that quantum mechanical gas-phase Gibbs free energies of formation compare reasonably well with the literature values with a root mean square (RMS) value of 1.83 kcal/mol for the selected compounds. These QM results are then used to determine the equilibrium constants (reported as log K) of aerosol-phase chemical reactions, including hydration reactions and aldol condensation for formaldehyde, acetaldehyde, acetone, butanal, hexanal, and glyoxal. Results are in qualitatively agreement with previous studies. In addition, the QM results for glyoxal reactions are consistent with experimental observations. To our knowledge, this is the first QM study that supports observations of atmospheric particle-phase reactions. Despite the significant uncertainties in the absolute values from the QM calculations, the results are potentially useful in determining the relative thermodynamic tendency for atmospheric aerosol-phase reactions. 46 3.2 Introduction Secondary organic aerosol (SOA) formation by gas/particle (G/P) partitioning has traditionally focused on low volatility products. The quantity of SOA formed can be estimated using absorptive or adsorptive G/P partitioning theory which assumes that this quantity is governed strongly by the vapor pressure of the compound as well as the liquid-phase activity coefficient [1-5]. Recent experimental work has suggested that the amount of SOA formed in a number of systems exceeds that based purely on G/P partitioning of low vapor pressure oxidation product [6-10]. Evidence also indicates that relatively volatile oxidation products, especially aldehydes and ketones, are being absorbed into the aerosol phase where they undergo aerosol-phase chemical reactions. The reaction products have relatively low vapor pressures compared to their parent compounds, which leads to additional partitioning from gas to particle phase, and hence, increases the organic particulate material (OPM). Aerosol-phase reactions, such as hydration, polymerization, hemiacetal/acetal formation, and aldol condensation, have been postulated as a means by which low volatility compounds can be formed thereby increasing the amount of OPM formed beyond that due to G/P partitioning of low vapor pressure gas-phase oxidation products alone. Theoretical studies by Barsanti and Pankow [11] have shown, however, that reactions such as hydration, polymerization, hemiacetal/acetal formation are not thermodynamically favorable under atmospheric conditions. Their results do suggest that aldol condensation may be thermodynamically favorable. These results seem to deviate from experimental observation [9]. On the other hand, they have shown that diol and subsequent oligomer formation are favorable for glyoxal [12], and these findings are 47 consistent with experiments [13-15] .While the experimental studies have provided convincing evidence that aerosol-phase chemical reactions (possibly acid- catalyzed) are involved to some extent in formation of SOA, uncertainty remains as to the likely aerosol-phase chemical reactions involving absorbed gas-phase organic compounds. The reactive uptake mechanism for relatively small, volatile organic compounds (short-chain aldehydes and ketones) is not well understood. Hydration is invariably the first step for volatile organics to dissolve into the particle phase, followed by various (possibly acid-catalyzed) reactions such as polymerization, hemiacetal/acetal formation, and aldol condensation. As suggested by Barsanti and Pankow [11], aldol condensation may be the most accessible reaction path for additional OPM formation. In the current study, we investigate the thermodynamic feasibility of various particle-phase heterogeneous reactions for some common atmospheric carbonyl compounds using quantum mechanical methods. In particular, we consider the hydration reaction and aldol condensation for small, short-chain aldehydes and ketones, such as formaldehyde, acetaldehyde, acetone, butanal, and hexanal. We also include glyoxal in our investigation. The relatively simple structure of glyoxal as well as its clear importance as an atmospheric oxidation product of a number of hydrocarbons makes it an excellent candidate for theoretical study. Recent studies [13-16] have shown aerosol growth by heterogeneous reactions of gas-phase glyoxal. The thermodynamic feasibility of a proposed glyoxal reaction pathway [14] is evaluated in our study. Similar to the studies by Barsanti and Pankow [11, 12], the thermodynamic analysis presented in this work is independent of the actual reaction pathway. The goal is to present a method to identify potential particle-phase reactions that may contribute to atmospheric OPM, and 48 the extent of OPM contribution if the reactions are kinetically favorable. This study does not yield any information regarding the kinetics; therefore, evidence of additional OPM formation may remain unobservable at short time scales. On the other hand, we focus only on determining the solution-phase equilibrium constant, K, which is the governing factor of the overall tendency of the reactions in the particle. Owing to the chemical complexity of atmospheric aerosols, one is often faced with the difficulty in obtaining physical properties for species for which limited experimental data exist. Here we apply novel quantum chemistry methods as an alternative predictive approach, which may reduce the number of parameters/experimental data required. 3.3 Computational Method In order to calculate equilibrium constants (K) of the reactions, the standard Gibbs free energy of a reaction (ΔG0 ) are needed and can be calculated by the standard Gibbs r free energy of formation (ΔG0). It is also related to the equilibrium constant (K) f according to the fundamental equation "G0 =$v "G0 =#RTlnK (3.1) r j f,j j where v is the stoichiometric coefficient for j in the reaction, ΔG0 is the standard i f,j ! Gibbs free energy of formation for j. ΔG0 can be determined using gas-phase heats of f,j formations and standard entropies (ΔG0 = ΔH0 - TΔS0). Then, as illustrated by Figure f,j f,j j 3.1, the free energy of reaction in aqueous solution, ΔG0 (aq), is related to the gas-phase r free energy of reaction, ΔG0 (gas) by adding the solvation energies of the species, ΔG0 . r s "G0(aq)="G0(gas, 1M)+#v "G0 (3.2) r r j s j ! 49 All the necessary quantities (ΔH0 , ΔS0 and ΔG0 ) can be determined by QM and f,j j s the procedure will be described below. It should be noted here that proper standard state conditions should be used for the free energy calculations in equation 3.2. The standard state for gas-phase reactions is 1 atm at 298K while the standard state for aqueous solution is 1 M at 298K. A brief description of the standard state conversion is provided in the Supporting information. "G0(gas) A + B # # #r ## $ C (gas) (gas) (gas) "G0(A) "G0(C) s s "G0(B) s "G0(aq) A + B # # #r ## $ C (aq) (aq) (aq) Figure 3.1: Thermodynamic cycle for computation of energy changes of reaction in the gas phase and solution. Adapted from Cramer (2004) [17] ! 3.3.1 Gas-phase Standard Heats of Formation. QM calculated standard heats of formation would be the first step to obtain free energies of formation for equation (3.1). QM heats of formation at 298K is given by the equation: "H0 (QM)= E + E +"H +"H* (3.3) f,j elec zpe (0#298K) 298K To obtain the ground state energy (E ), the molecules were optimized using elec ! DFT/X3LYP, with a fairly large basis set aug-cc-pVTZ(-f) [18] in the gas phase. For molecules with different conformations, geometry optimization using a smaller basic set 6-31g** [19] was carried out to select 4 (or 5) relatively stable conformations. Further optimizations were then performed for the selected structures at the higher level to calculate E . The vibrational frequency calculations are performed for the most stable elec 50 structure at the HF/631g** level in the gas phase with a scaling factor 0.8992. This scaling factor was determined by comparing the theoretical harmonic vibrational frequencies with the corresponding experimental values utilizing a total of 122 molecules (1066 individual vibrations) and a least-squares approach [20]. The frequency calculations ascertain the structure of the molecules and provide zero point vibrational energy (E ), the thermal vibrational, rotational, and translational enthalpy from 0 K to zpe 298 K (ΔH ). The last term in equation 3.3, ΔH* , allows for corrections between (0-298K) 298K the theoretical heats of formation, referenced to the electrons and nuclei separated an infinite distance, and experimental values, referenced to the elements at standard temperature and pressures. A common reference is the enthalpies of formation for the neutral atoms in the gas phase, so: n "H* =$n (h0#hQM) (3.4) 298K i i i i=1 where n is the number of elements in the compound, n is the number of atoms of each i ! element. h0 is the experimental atomic heats of formation in the gas phase [21], and hQM i i is the theoretical value: 3 1 hQM = Eelec + RT"PV = Eelec + RT (3.5) i i 2 i 2 where Eelec is the quantum electronic energy, 3/2 RT is the translational energy at 298K i ! and the ideal gas law is applied for PV, R is the gas constant, and T is the temperature. Large discrepancies can be found between the QM calculated heats of formation using equation 3.3 – 3.5 and the experimental heats of formation. Improvement can be made by applying the correction scheme such as the J2 model of Dunietz et al. [22] to the 51 enthalpies. The J2 model is based on the generalized valence bond-localized Møller - Plesset method (GVB-LMP2), and it uses a three parameter correction term composed of σ bond and π bond parameters and an additional parameter to account for the difference of lone pairs between the molecule and the separated atoms. Recently, Blanco and Goddard [23] developed a correction scheme closely following the J2 corrections [22], but their “chemical bond” scheme for high level corrections (CBHLC) applies to DFT quantum calculations and only consists of the two σ bond and π bond parameters. Details of the CBHLC scheme are given in the Supporting Information. All calculations were performed using the Jaguar 6.0 package [24]. QM calculations were carried out using DFT with X3LYP. X3LYP is an extended hybrid density functional that has been shown to be an accurate and practical theoretical method [25-27], with a particularly accurate estimation of van der Waals interactions. 3.3.2 Standard Free Energies of Formation and Equilibrium Constants Adding the entropy term to the QM gas-phase heats of formation will give us the standard Gibbs free energy, ΔG0 , at 298K: f % n ( "G0 ="H0 (QM)#T’ S (QM)#$nS * (3.6) f,j f,j j,298K i i,298K & ) i=1 where S (QM) is the entropy of the compound at 298K and S is the entropy for j,298K j,298K ! the elements in their reference states [21]. At this point, we can obtain the gas-phase Gibbs free energies of reaction using equation 3.1. In order to obtain solution-phase energies, the solvation effect is accounted by the free energy of solvation, (ΔG 0), as s shown in Figure 3.1. 52 3.3.3 Solution phase energy QM determined solvation energies, ΔG 0, for both the parent compounds and s reaction products are shown in Table 3.2. Solvation energy (ΔG 0) describes the s interaction of a solute with a surrounding solvent. Change in entropies due to conformational changes of the molecules from gas phase to aqueous phase is minimized by re-optimization of the compounds at X3LYP/cc-pVTZ(-f) using the implicit continuum solvent model [28-30]. In a continuum model, solute atoms are treated explicitly and the solvent is represented as a continuum dielectric medium. The electrostatic contribution to the solvation free energy is computed using the Poisson- Boltzmann (PB) method (see Supporting Information). The PB equation is valid under conditions where dissolved electrolytes are present in the solvent, but the current implementation of the QM solvation model is confined to zero ionic strength. The effect on the equilibrium constant for reactions in solution with dissolved electrolytes remains unknown because the change in solvation effect due to ionic strength of the solvent varies with species. Although the solvation model is limited to solution zero ionic strength, we can apply the solvation model to any solvent with the proper choice of dielectric constant. Water is chosen in this study because it is often the most important component in aerosols. 3.4 Results and Discussion Calculated QM free energies of formation in the gas phase, ΔG0 (QM), and f estimated free energies of formation obtained from the group contribution “Joback method” [31], ΔG0 (Joback). Both estimations are compared to the gas-phase literature f values [32], and the differences (δG) are calculated (see Table 3.1). 53 Table 3.1: Quantum calculated free energies of formation. The literature free energies of formation and group contribution (Joback) method estimated values are included for comparison. Absolute errors (δG) are also shown. All values are in kcal/mol ΔG0 (exp)a ΔG0 (Joback) δGc ΔG0 (QM) δGd f f f [31] Water -54.60 ± 0.11 --b -- -53.92 0.68 Formaldehyde -24.51 ± 0.74 --b -- -24.96 -0.45 Acetaldehyde -31.84 ± 0.96 -31.90 -0.06 -33.68 -1.84 Acetone -36.14 ± 1.08 -36.91 -0.77 -39.29 -3.16 Butanal -27.78 ± 0.83 -27.88 -0.10 -29.06 -1.28 Hexanal -23.90 ± 0.72 -23.86 0.04 -23.72 0.18 Glyoxal -45.31 -55.67 -10.36 -45.23 0.07 Ethylene glycol -72.08 ± 2.16 -73.49 -1.40 -69.41 2.67 Hydroxyacetone -68.07 ± 6.81 -69.59 -1.52 -70.04 -1.97 2,4-pentanedione -64.37 ± 0.64 -63.68 0.69 -61.77 2.60 Glutaraldehyde -48.58 ± 2.43 -49.64 -1.06 -46.61 1.97 Cyclopropane -58.49 ± 1.75 -71.74 -13.25 -56.80 1.69 carboxylic acid a. Gas-phase literature values and uncertainties of free energies of formation, ΔGf0(exp) are obtained from DIPPR database [32]. b. The group contribution Joback method [31] lacks sufficient groups to estimated ΔGf0 for water and formaldehyde. c. δG = ΔGf0(Joback) – ΔGf0(exp) d. δG = ΔGf0(QM) – ΔGf0(exp). Group contribution methods are widely used because of their efficiency and accuracy. Barsanti and Pankow [11] found that the estimated ΔG0 by the Joback method f agrees with the values from Yaws [33] to ±0.7 kcal/mol (±3 kJ/mol) for a set of compounds including alcohols, aldehydes and ketones. However, as for all predictive methods, group contribution methods are limited by the experimental data that were used in the parameterization, which may be problematic for multi-functional compounds where experimental data are scarce. For the selected 12 compounds, the root mean square (RMS) value for errors between ΔG0(Joback) and ΔG0(exp) is 5.80 kcal/mol (vs. 0.90 f f