Estimation of Radiative Efficiency of Chemicals with Potentially

Dec 8, 2015 - ... warming potential estimates for the C1–C3 hydrochlorofluorocarbons (HCFCs) included in the Kigali Amendment to the Montreal Protoc...
0 downloads 0 Views 760KB Size
Download PDF
Article pubs.acs.org/est

Estimation of Radiative Efficiency of Chemicals with Potentially Significant Global Warming Potential Don Betowski,*,† Charles Bevington,‡ and Thomas C. Allison§ †

U.S. Environmental Protection Agency, National Exposure Research Laboratory, Environmental Sciences Division, P.O. Box 93478, Las Vegas, Nevada 89193-3478, United States ‡ U.S. Environmental Protection Agency, Office of Pollution Prevention and Toxics, Risk Assessment Division, 1200 Pennsylvania Avenue, N. W., Mail Code: 7408M, Washington, D.C. 20460, United States § Material Measurement Laboratory, National Institute of Standards and Technology, 100 Bureau Drive, Stop 8320, Gaithersburg, Maryland 20899-8320, United States S Supporting Information *

ABSTRACT: Halogenated chemical substances are used in a broad array of applications, and new chemical substances are continually being developed and introduced into commerce. While recent research has considerably increased our understanding of the global warming potentials (GWPs) of multiple individual chemical substances, this research inevitably lags behind the development of new chemical substances. There are currently over 200 substances known to have high GWP. Evaluation of schemes to estimate radiative efficiency (RE) based on computational chemistry are useful where no measured IR spectrum is available. This study assesses the reliability of values of RE calculated using computational chemistry techniques for 235 chemical substances against the best available values. Computed vibrational frequency data is used to estimate RE values using several Pinnock-type models, and reasonable agreement with reported values is found. Significant improvement is obtained through scaling of both vibrational frequencies and intensities. The effect of varying the computational method and basis set used to calculate the frequency data is discussed. It is found that the vibrational intensities have a strong dependence on basis set and are largely responsible for differences in computed RE values.

■

INTRODUCTION The Global Warming Potential (GWP) was first introduced by the Intergovernmental Panel on Climate Change (IPCC) as a simple metric for comparison of possible greenhouse gases in terms of their impacts on the climate system.1 There are wellestablished limitations with GWP, and other promising metrics are being explored.2 Despite its limitations, GWP continues to be widely used to assess relative impacts of various chemical substances.3 GWP is calculated using three primary inputs: molecular weight, atmospheric lifetime, and radiative forcing (RF). Radiative forcing is defined as the change in the net, downward minus upward, radiative flux (expressed in W m−2) at the tropopause or top of atmosphere due to a change in an external driver of climate change. Radiative efficiency (RE) is the radiative forcing associated with a 1 ppbv change in the atmospheric concentration of a gas and has units of W m−2 ppbv−1. Atmospheric mixing ratio information is needed to convert RE to RF. RE is commonly computed by integration of the experimental infrared (IR) spectrum, and some of the best available data, including some of the data against which the present results will be compared, arise from this method. This © XXXX American Chemical Society

article focuses on estimation of RE using quantum chemistry techniques, though other estimation schemes can be employed. For example, Young et al. used structure−activity relationships to predict the radiative efficiency of fluorinated ethers.4 A group increment scheme was employed by Kokkila et al. to predict the integrated IR intensity within the atmospheric IR window, and could be used to compute RE.5 Most estimation schemes for RE involve the use of computational chemistry. Kazakov et al. used the semiempirical PM6 method to compute the RE (and GWP) for over 56 000 compounds to screen for improved refrigerant fluids.6,7 Their use of a semiempirical method was based on the large number of compounds they needed to evaluate. Most researchers have used either the ab initio MP2 method or a density functional theory (DFT) method (typically B3LYP). The means for computing RE from computed IR spectra is based on the pioneering work of Pinnock et al., who created a Received: August 27, 2015 Revised: December 7, 2015 Accepted: December 8, 2015

A

DOI: 10.1021/acs.est.5b04154 Environ. Sci. Technol. XXXX, XXX, XXX−XXX

Article

Environmental Science & Technology Table 1. Data Sources and Common Characteristics of chemical substances that may have high GWP data sources used to identify 235 chemicals*

IPCC 5

number of chemicals Common characteristics across 235 chemicals volatility (vapor pressure, v) lifetime (atmospheric lifetime, l; OH rate constant, k; number of hydrogens) molecular weight g/mol percent of bonds in chemical that are halogen-carbon bonds and/or hydrogen−carbon bonds (%)

199 Most v > 10 Torr l > 0.5 yr or no H atoms g/mol 50%

Hodnebrog 2013

Other Data Sources**

213 Few 0.01 ≤ v ≤ 10 Torr l > 0.1 yr or k > 4.3e−13

21 Occasional or none v < 0.01 Torr l < 0.01 yr or k > 1e−12

100 < g/mol % > 10%

g/mol >500 % < 10%

*

There was substantial overlap of chemicals across data sources. **For chemicals not covered by IPCC5 or Hodnebrog 2013, includes previous IPCC reports, WMO reports, and open literature (see SI Table 1).

geometry of the chemical substance (ideally the global minimum). These are readily obtained from a vibrational frequency calculation. Such calculations were made using the B3LYP density functional theory (DFT) method. Two Popletype basis sets were used for these calculations, the smaller 631G* basis set and the larger 6-311++G** basis set. In cases where a compound contained iodine, the MIDI or the aug-ccpVDZ-PP (for I) and aug-cc-pVDZ (for all other atoms) basis sets were used in place of the smaller and larger basis sets, respectively. The geometry of each chemical substance was fully optimized at the stated level of theory, and vibrational frequencies and intensities were calculated at this geometry. The vibrational frequencies were inspected to ensure that the optimized geometry was a minimum on the potential energy surface, that is, all frequencies were positive real values. In all calculations, the global potential energy minimum structure was sought. However, especially for larger structures with many minima, it is impossible to guarantee that this structure was found in each case. All structures were carefully inspected to ensure that a chemically reasonable structure was found. In cases where multiple minima were likely to be present, optimizations were initiated from several different geometries and the lowest energy geometry was selected. The use of higher energy structures was found to increase the value of the RE, so any outliers in the data set were reexamined to ensure that no lower energy structures could be found. All calculations were carried out using the Gaussian 09 computational chemistry program.17 Additional calculations carried out using different methods and basis sets were used to characterize the range of RE values that might be produced via computational methods and are described below. Calculation of Radiative Efficiency. In this article, values of RE are computed using the values of the vibrational frequencies and the associated intensities taken from the B3LYP calculations. These values were used with the model of Pinnock et al.,8 which permits calculation of RE from the computed IR spectrum.9 This method has been shown to give results within 0.3% of a more rigorous narrowband model and has been widely applied for the computation of RE using both measured and computed IR frequency data.8 The model incorporates the discrete, frequency-dependent IR greenhouse gas (GHG) absorption and the discrete frequency dependent IR average earth sky emission. In this model, the RE per unit cross section for the average earth sky (i.e., including clouds) is parametrized in 10 cm−1 intervals (from 0 to 2500 cm−1) for a 1-part-per-billion (1 nmol/mol) by volume increase in the atmospheric concentration for the GHG1

simple, but accurate, procedure based on the Planck function and a global and annual mean atmosphere (GAM) model.8 Papasavva et al. first demonstrated the use of this method with computed IR spectra using frequencies computed at the MP2/ 6-31G** level of theory,9 building on their earlier work validating the use of quantum chemistry methods for the calculation of vibrational frequencies and intensities.10,11 Blowers et al. calculated RE values for 52 hydrofluoroether compounds at the B3LYP/6-31G* level of theory, comparing 27 values to experimental results.12 In a subsequent study, Blowers et al. computed RE values (as well as GWP) for hydrofluoroethers at the same level of theory.13 Bravo et al. compared computational methods for determining RE (and GWP) to experiments for perfluorocarbons.14 A subsequent work by Bravo and co-workers computed RE at the B3LYP/631G** level of theory for hydrofluoroethers and hydrofluoropolyethers.15 Although any saturated halogenated chemical may have significant GWP, and fully halogenated halocarbons certainly do, many chemical substances containing halogens may or may not have high GWP. The terms halogenated chemical and halocarbon are used as short hand for fluorinated, chlorinated, and brominated chemicals. Iodinated chemicals tend to have high RE but short atmospheric lifetimes limiting their GWP. The time and expense involved in measuring the RE values for all chemical substances with potential GWP would be prohibitive. Where measurement is considered part of a statutory or regulatory mandate, predictive methods to screen for likely candidates for further testing could be an important assessment and cost-saving tool. Thus, it is important to find reliable methods for estimating RE and RF and validating them. This is the goal of the present publication. Values of RE published by the IPCC, WMO, EPA, and other organizations, and especially from the recent work by Hodnebrog et al.16 who extensively reviewed the available determinations, are considered to be the best available sources in this article and form the basis for comparison to computed results. It is important to note, as discussed above, RE is only one of three parameters used to calculate GWP. The GWP of a chemical also depends on the molecular weight and atmospheric lifetime of the chemical. Where atmospheric lifetimes are short (i.e., measured in days to weeks), GWPs are likely to be low irrespective of RE. This is especially important for chemicals containing iodine and for unsaturated halogenated chemicals.

■

MATERIALS AND METHODS Quantum Chemistry Calculations. The essential quantities needed for the calculation of RE are the vibrational frequencies and associated intensities at a minimum-energy B

DOI: 10.1021/acs.est.5b04154 Environ. Sci. Technol. XXXX, XXX, XXX−XXX

Article

Environmental Science & Technology RE GHG = w ∑ σiGHGFiσ

(1)

i GHG

improvement in the values of the RE using this approach, the conventional narrow spectral feature approach will be used in all results presented in this article. Radiative Efficiency. REs were calculated for the full set of 235 substances using the B3LYP method and two different basis sets, 6-31G* and 6-311++G**, as described above, with no scaling applied to improve agreement with the best available values. The results are plotted against the best available values in Figure 1. From the figure it is clearly seen that the results

−2

−1

where RE is the RE of the GHG (W m ppbv ), w is the width of the interval, σGHG is the ith (over a particular 10 cm−1 i interval) absolute integrated absorption cross section for the GHG, and Fσi is the ith average earth sky emission factor over the same 10 cm−1 interval. In this work, three GAM models are studied: the Pinnock model (10 cm−1 interval),8 and two newer models (1 and 10 cm−1 intervals).16 Table 1 identifies the data sources used to identify a set of 235 chemical compounds with available RE values as well as common characteristics across the chemicals. Values of the RE for 235 chemicals are compiled in Table S1 (Supporting Information). In general, these chemicals have low molecular weight, are volatile, tend to be long-lived in the atmosphere, and have many carbon−halogen and carbon−hydrogen bonds. This set of chemical substances was chosen because it comprises the best available values of the RE and forms the basis for comparison with values derived from computational methods. Some care must be taken when collecting values of the RE for comparison purposes. For example, the supplementary tables from the work of Hodnebrog et al., from which many values in the comparison set were taken, give four different values of the RE.16 These values differ in the GAM profile used in their computation, in whether the results are instantaneous or constant profile, and in whether a lifetime correction was applied. Additionally, the RE reported by Hodnebrog et al.16 was increased by 10% to account for the stratospheric temperature adjustment. In this article, the same procedure is applied, that is, all RE values computed in the present study are increased by 10%, even though it generally increases the error versus the best available values.

■

RESULTS AND DISCUSSION Choice of GAM Atmospheric Model. The Pinnock model has been tested and applied extensively and has been shown to produce reliable values of the RE.9 In their paper, Hodnebrog et al.16 introduced two new models (generated using the Oslo line-by-line model18) with bandwidths of 10 and 1 cm−1. These models were shown to give values very similar to those produced by the Pinnock model in most cases (using a more up-to-date GAM model), and the present results confirm this with a maximum deviation of less than 2% between the three models for the set of 235 compounds. Results computed using the Oslo 10 cm−1 model (with the Pinnock method) are reported in this article, with the understanding that the results from the other models are similar. The computation of RE involves considering each vibrational frequency to be infinitely narrow. This single value of the vibrational frequency is then used to select an appropriate interval in the RE function, and the value of the RE function in this interval is multiplied by the corresponding vibrational intensity. The sum of these contributions forms the RE (eq 1). In reality, the vibrational frequency spectral features are broadened due to a number of effects. Calculation of the RE was tested using Gaussian broadened spectra and found that there is no strong effect on the results, consistent with the results of Bravo et al.14 and Hodnebrog et al.16 In the present results, a broadening factor of 7 cm−1 (full width at halfmaximum) produced the best results, whereas Bravo et al. recommended a width of 14 cm−1.14 As there is no significant

Figure 1. Plot of calculated versus best available values of the radiative efficiency for 235 chemicals. Ideal agreement is indicated by the solid black line. The line of best fit is indicated by the dashed blue line, and the equation of this line is given on the plot. Units of RE are W m−2 ppbv−1.

produced using the smaller basis set are significantly better than those produced using the larger basis set. To quantify this, a linear fit to the data was made for both data sets (this fit was forced to go through the origin). The slope of the fit made using the smaller basis set results is 0.883, indicating good agreement with the best available values (note that since the computed results generally overestimate the best available values, the stratospheric temperature adjustment that is applied increases the discrepancy). The slope for the larger basis set result is 0.633, indicating that the agreement with the best available values is worse. This discrepancy is attributable solely to the choice of basis set. This trend is the opposite of what one might expect (i.e., larger basis sets should yield a more accurate wave function and might be expected to produce more accurate results for derived properties), but is a trend that will be observed throughout the results presented in this article. C

DOI: 10.1021/acs.est.5b04154 Environ. Sci. Technol. XXXX, XXX, XXX−XXX

Article

Environmental Science & Technology The mean absolute error (MAE) and root mean squared error (RMS) values for the smaller basis set results are 0.044 W m−2 ppbv−1 and 0.075 W m−2 ppbv−1, respectively, when the multiplicative correction described above is applied. This procedure illustrates that multiplicative scaling is necessary (and sufficient) to correct systematic defects in the computed results. A more physically motivated procedure to scale computed results will be developed later in the article. To illustrate the effect of correcting the data in the manner described above, the corrected data for the B3LYP/6-31G* results are plotted in Figure 2 where it can be seen that most of

Figure 3. Plot of the percent error of computed radiative efficiency values as a function of the best available value for the smaller (blue bars) and larger (red bars) basis sets. Units of RE are W m−2 ppbv−1.

comparison to previous work. The PBE functional20 is another frequently used DFT method, and the B97D3 method21 incorporates a van der Waals correction. Finally, the MN12-L method of Truhlar and co-workers is included as an example of a recently developed DFT method.22 Computed values of RE for 20 chemical substances are shown in Table 2 for all five computational methods. It is clear Table 2. Statistics Measures of the Error in Computed Values of the Radiative Efficiency (W m−2 ppbv−1) versus the Best Available Values Compiled by Computational Method for Calculations Made at Various Levels of Theory with Varying Basis Sets, where n is the Number of Values, MSE is the Mean Signed Error, MAE is the Mean Absolute Error, RMS is the Root Mean Squared Error, and σ is the Standard Deviation

Figure 2. Plot of computed (x axis, calc, B3LYP/6-31G*) and best available (y axis, expt) values of the radiative efficiency (blue circles). The solid black line indicates perfect agreement, and the red shaded region indicates an uncertainty in the radiative efficiency of 30%. Units of RE are W m−2 ppbv−1.

the data lies within the 30% uncertainty that has been previously estimated for values of the radiative efficiency (some recent work indicates that this figure may be smaller16). For the larger basis set results, the MAE and RMS values are 0.038 W m−2 ppbv−1 and 0.079 W m−2 ppbv−1, respectively. This suggests that the larger error seen in the large basis set results can be effectively corrected through scaling. (In this article, scaling of both vibrational frequencies and vibrational intensities is recommended as will be seen below.) Finally, it is noted that the parameters of the scaling correction for the computed data must be derived for each combination of method and basis set used. Another means of comparing the performance of the small and large basis set results is shown in Figure 3, where the percent error of the computed results at the best available value of the RE is given. Again, the superior performance of the smaller basis set, in the absence of any correction, is seen. Effect of Computational Details. To study the effect of varying the computational method and the basis set, a set of 20 chemical substances was randomly selected from the set of 235. Calculations were made using five different computational methods and seven different basis sets. The results are described in the following subsections. Choice of Computational Method. Clearly, the choice of computational method (and basis set, as discussion in the following section) can have a significant effect on the results.19 To investigate this effect, five representative computational methods were selected. The B3LYP DFT method and the ab initio MP2 method have been used in previous studies,14 and are included as they are reasonable and popular choices and for

method

n

MSE

MAE

RMS

σ

MP2 B3LYP PBE B97D3 MN12-L

127 140 140 140 140

0.070 0.102 0.104 0.101 −0.015

0.072 0.106 0.108 0.108 0.047

0.094 0.142 0.148 0.145 0.063

0.062 0.100 0.106 0.104 0.063

from this table that the MP2 method is superior to the B3LYP method as has been seen in previous studies.14 The PBE and B97D3 methods perform about the same as B3LYP. The MN12-L method performs about as well as the considerably more computationally demanding MP2 method. This is perhaps unsurprising as the MN12-L method showed improved performance on a wide variety of properties.22 This result indicates that further study of DFT methods is warranted, and that the MN12-L method may be recommended until such a study is performed. The choice of computational method has been examined by Bravo et al.,14 who found that MP2 results were superior to B3LYP results, in accord with the present findings. However, it is important to note that the MP2 method can be significantly more computationally expensive than B3LYP, and that the use of B3LYP or another method (e.g., MN12-L as considered above) may be necessary. Choice of Basis Set. The effect of basis set was studied using seven different basis sets on all 20 chemical substances using five different computational methods. The choice of basis D

DOI: 10.1021/acs.est.5b04154 Environ. Sci. Technol. XXXX, XXX, XXX−XXX

Article

Environmental Science & Technology

employed in all of the results considered in the remainer of this article. Thus, errors in other quantities will have to be examined to explain differences from the best available values. Other vibrational frequency scaling schemes have been employed for computing RE from computational results. For example, Papasavva et al.9 and Bravo et al.14,15 employed a linear correction (e.g., y′ = my + b) to obtain improved vibrational frequencies from computed results to predict REs. The multiplicative method is used in this article as it is the most well characterized method for adjusting computed vibrational frequencies.24 Scaling of Vibrational Intensities. To understand the origin of the discrepancy between the small and large basis set RE results, the mean unsigned percent deviation (MUPD) of the vibrational frequencies and intensities was computed for frequencies in the “atmospheric window” (800 cm−1 to 1200 cm−1) with intensity greater than 0.0001 km/mol (a total of 1234 frequencies). It is found that the MUPD for the frequencies is 3%, while the MUPD for the intensities is 69%. Thus, differences in the vibrational intensities are almost entirely responsible for the large difference in agreement with the best available results between the small and large basis sets. This difference is clearly seen in the histograms presented in Figure 4. The source of this difference is due to the locality of the basis sets used. The smaller 6-31G* basis set has fewer basis functions, and these are localized near the atomic center. The larger 6-311++G** basis set adds diffuse basis functions that

sets was designed to span a range of sizes and types. The 631G*, 6-31G**, and 6-311++G** Pople basis sets, and the ccpVDZ, aug-cc-pVDZ, cc-pVTZ, and aug-cc-pVTZ Dunning basis sets were employed for the calculation of RE, and the results are presented in Table 3. The surprising result is that Table 3. Same as Table 2, but Statistics are Compiled by Basis Set. Units are W m−2 ppbv−1 basis set

n

MSE

MAE

RMS

σ

6-31G* 6-31G** 6-311++G** cc-pVDZ aug-cc-pVDZ cc-pVTZ aug-cc-pVTZ

100 100 98 98 98 97 96

0.047 0.048 0.118 0.071 0.090 0.073 0.082

0.068 0.069 0.122 0.083 0.095 0.087 0.095

0.092 0.093 0.167 0.114 0.134 0.119 0.132

0.080 0.081 0.119 0.089 0.100 0.095 0.104

smaller basis sets, 6-31G* and 6-31G**, outperform the larger and correspondingly more computationally demanding basis sets as seen previously. This trend is also seen in the Dunning basis sets where addition of diffuse functions (denoted aug-) leads to increased error versus the best available values. Between the two best performing basis sets, 6-31G* and 631G**, it is seen that the inclusion of polarization functions on the hydrogen atoms in the second basis set has very little effect on the results. Effect of Optimized Geometry. When optimizing structures for RE calculations, it is important to make sure that the lowest-energy optimized geometry is used. In several cases, it was found that improvement with the best available value was obtained when a lower energy structure was used in the calculations. A better way to compute the RE is to compute several minimum energy structures and compute the RE as a (Boltzmann) weighted average of values due to the various structures as has been done for thermodynamic quantities.23 Bravo et al. used a Boltzmann weighted average of the contributions from various conformers of n-perfluorobutane and found that the difference from the value of the RE based solely on the minimum energy geometry was less than 1%.14 This approach was not employed in the present study. It is noted that it is always important to find the lowest energy structure as this frequently has the largest effect on the computed value of the RE. When higher energy structures were used to calculate the RE, the value was always higher compared to the minimum energy structure. The best method for determining the value of the RE, when multiple conformers are present is beyond the scope of this article, but is certainly a topic that warrants further investigation. Scaling of Vibrational Frequencies. It is well-known that computed frequencies are typically too large. It is common practice to correct for this overestimation using a multiplicative scaling factor leading to significantly improved agreement with experimental spectra.24 The use of such a scaling factor in the computation of RE implies that the value of the GAM function will change as a different interval (in the Pinnock method) is selected by the modified vibrational frequency. Scaling factors of 0.960 for the 6-31G* basis set and 0.967 for the 6-311+ +G** basis set were used.25 It is found, somewhat surprisingly, that scaling of the vibrational frequency negatively impacts agreement with the accepted RE values. However, given that multiplicative scaling of vibrational frequencies increases agreement with experimental spectra,24 this procedure is

Figure 4. Histogram plot of absolute percent deviations of vibrational frequencies (top panel) and intensities (bottom panel) between small and large basis set calculations for frequencies in the atmospheric window (800 cm−1 to 1200 cm−1). E

DOI: 10.1021/acs.est.5b04154 Environ. Sci. Technol. XXXX, XXX, XXX−XXX

Article

Environmental Science & Technology

over the full set of 235 substances are broken down by chemical class in the following subsection. Analysis by Chemical Class. To further quantify the performance of quantum chemistry methods for the calculation of RE, the set of 235 substances were assigned to various groups as presented in Table 5 for REs computed using both

make significant contributions farther from the atomic center. Since the vibrational intensities are based on atomic orbital overlaps, the larger basis set generally yields larger vibrational intensities. Though it is commonplace to scale vibrational frequencies,24 there is no widely used scaling for vibrational intensities. However, if it is assumed that the vibrational intensities can be scaled using a simple multiplicative factor, then the fitting procedure employed above can be viewed as a more physically based scaling procedure, yielding improved results. Multiplicative scaling of vibrational intensities has been employed by Halls and Schlegel.26 Given that scaling of vibrational frequencies is known to improve agreement with experimentally measured spectra,24 and having seen that the major source of the difference between basis sets is largely due to the vibrational intensities, it is recommended that both types of scaling are employed to produce the highest quality results. Note that scaling factors for vibrational intensities are not readily available and must be derived for each combination of computational method and basis set, and that the quality of the final result will depend on the quality of the data set used to derive the vibrational intensity scaling factor. Analysis by Substance. The performance of quantum chemistry methods may also be analyzed by considering the errors of individual substances versus the best available values. Results of this type are presented for five quantum chemistry methods and seven basis sets for 20 substances in Table 4. It is

Table 5. Statistics for Error of Calculated Values of Radiative Efficiency from the Best Available Values Organized by Chemical Classa 6-31G*

n

MSE

MAE

RMS

σ

HCFC-22 HFC-134 HCFC-21 perfluoropropyl acetate HFE-245fa2 HFE-254cb1 difluoromethyl-2,2difluoroacetate PFC-71−18 2,2,3,3-tetrafluoro-1-propanol methylene chloride HFC-152 HFC-143a fluoromethyl carbonofluoridate (E)-HFC-1234ze PFC-14 PFC-116 3,3,4,4,5,5,6,6,7,7,7undecafluoroheptan-1-ol allyl_2,2,2-trifluoroacetate propane ethane

35 35 35 33 35 35 35

0.032 0.041 0.008 0.171 0.092 0.065 0.109

0.034 0.043 0.037 0.211 0.112 0.074 0.120

0.039 0.048 0.045 0.226 0.134 0.081 0.133

0.022 0.025 0.045 0.150 0.098 0.048 0.078

30 35 35 35 35 35 35 35 35 30

0.181 0.060 0.012 −0.002 0.058 0.096 0.051 0.126 0.109 0.250

0.241 0.060 0.031 0.015 0.073 0.098 0.056 0.132 0.128 0.250

0.259 0.061 0.038 0.017 0.082 0.112 0.059 0.165 0.139 0.263

0.188 0.014 0.037 0.017 0.059 0.058 0.031 0.108 0.089 0.085

34 35 35

0.090 0.000 0.001

0.104 0.001 0.001

0.115 0.001 0.001

0.072 0.001 0.001

a

n

MAE

RMS

MAE

RMS

hydrocarbons fluorohydrocarbons perfluorocarbons chlorohydrocarbons perchlorocarbons bromohydrocarbons contains halogen(s) contains fluorine contains chlorine contains bromine contains iodine contains nitrogen contains oxygen contains sulfur

12 39 20 5 1 2 110 189 30 11 2 1 118 1

0.004 0.067 0.136 0.026 0.078 0.024 0.071 0.122 0.077 0.045 0.051 0.252 0.152 0.099

0.005 0.088 0.149 0.040 0.078 0.030 0.096 0.193 0.096 0.056 0.058 0.252 0.233 0.099

0.006 0.167 0.267 0.026 0.082 0.022 0.144 0.238 0.121 0.078 0.143 0.753 0.287 0.151

0.009 0.207 0.284 0.039 0.082 0.028 0.188 0.309 0.146 0.090 0.145 0.753 0.358 0.151

a

Vibrational frequencies are scaled, but vibrational intensities are not. Units are W m−2 ppbv−1.

Table 4. Same as Table 2, but Statistics Are Compiled by Chemical Compounda chemical

6-311++G**

chemical classification

basis sets. It is immediately apparent that the best results are obtained for hydrocarbons. Chloro- and bromohydrocarbons also seem to be treated well at this level of theory (although the sample size is small). The fluorocarbons exhibit a considerably larger (absolute) error from the accepted values, with perfluorocarbons showing an even greater (absolute) error. Note that while the absolute errors for these compounds are larger, the RE values are also larger. This leads to relative percent errors that are comparable to compounds with smaller values of RE and is shown in Figure 3. When the results are examined by whether a particular atom is present (lower half of Table 5), it is seen that the errors are larger for chemicals that contain F, N, O, and S (note that the sample size is one for N and S). When the statistics listed in Table 5 are examined for the B3LYP/6-31G* results, the errors are reduced in almost all categories. The largest (absolute) errors are observed in chemicals containing F and O. Thus far, the results presented have employed overall scaling (Figures 2 and 3) or scaling of the vibrational frequencies only to facilitate understanding of the source(s) of discrepancies. To illustrate the efficacy of using both vibrational frequency and vibrational intensity scaling at the same time (which is recommended), the data organized by chemical class with both scalings are presented in Table 6 (the vibrational intensity scaling factor was 0.699 for the B3LYP/6-31G* results and 0.570 for the B3LYP/6-311++G** results). It is clear from comparison of Tables 5 and 6 that the vibrational intensity scaling significantly improves the results, especially for compounds containing oxygen and fluorine. For the full set of 235 compounds, using the vibrational intensity scaling reduces the MAE and RMS error from 0.112 W m−2 ppbv−1 and 0.179 W m−2 ppbv−1, respectively, to 0.045 W m−2 ppbv−1 and 0.082 W m−2 ppbv−1, respectively for the B3LYP/6-31G*

Units are W m−2 ppbv−1.

clearly seen that the results vary considerably by chemical, with more than 2 orders of magnitude difference in the MAE of the RE, for example. Simple hydrocarbons exhibit smaller errors, whereas chemical substances containing oxygen exhibit the largest errors. To study this effect more closely, the errors seen F

DOI: 10.1021/acs.est.5b04154 Environ. Sci. Technol. XXXX, XXX, XXX−XXX

Article

Environmental Science & Technology

311++G**. This results helps explain the success of earlier work based on B3LYP or MP2 calculations using the 6-31G** basis set. This is an important point, as the agreement found was fortuitous rather than a consistent characteristic of results derived from quantum chemistry computations. In general, each quantum chemistry method and basis set employed for computing RE must be evaluated for accuracy as results can differ significantly. The common practice of scaling vibrational frequencies to improve agreement with the experimentally measured IR spectrum leads to decreased agreement with the best available RE results as does the application of a stratospheric temperature adjustment. However, the vibrational intensities are found to be the main source of error in the calculated results, and scaling of both vibrational frequencies and intensities is an effective means of producing reliable RE values. A newer DFT method, MN12-L, gives results that are comparable to the computationally more expensive MP2 method, and that are superior to the other methods studied (when considering the statistical measures listed in Table 2). These findings are important for those using quantum chemistry methods to calculate RE values, and the results in this article suggest a set of best practices for producing rapid and reliable results. A set of RE values derived according to the procedures outlined in this article is provided for more than 1200 chemicals in Supporting Information Table S-2. These values could be used with current and future experimentally derived RE values to refine estimation approaches for RE and RF. These values could also be used in combination with available experimentally derived or estimated atmospheric lifetimes to consider relative global warming potential. Overall, it is found that quantum chemistry methods are appropriate for rapid and reliable calculation of RE, yielding results that are generally within the uncertainty of experimentally based approaches when appropriate scaling factors are applied.

Table 6. Statistics for Error of Calculated Values of Radiative Efficiency from the Best Available Values Organized by Chemical Classa 6-31G*

6-311++G**

chemical classification

n

MAE

RMS

MAE

RMS

hydrocarbons fluorohydrocarbons perfluorocarbons chlorohydrocarbons perchlorocarbons bromohydrocarbons contains halogen(s) contains fluorine contains chlorine contains bromine contains iodine contains nitrogen contains oxygen contains sulfur

12 39 20 5 1 2 110 189 30 11 2 1 118 1

0.006 0.027 0.025 0.033 0.104 0.016 0.029 0.051 0.077 0.038 0.115 0.083 0.060 0.095

0.008 0.033 0.044 0.047 0.104 0.018 0.043 0.090 0.096 0.052 0.118 0.083 0.108 0.095

0.005 0.030 0.030 0.040 0.139 0.010 0.032 0.049 0.042 0.041 0.034 0.076 0.055 0.124

0.007 0.046 0.046 0.055 0.139 0.014 0.048 0.083 0.050 0.058 0.044 0.076 0.097 0.124

a

Vibrational frequencies and vibrational intensities are both scaled. Units are W m−2 ppbv−1.

basis set results, and from 0.217 W m−2 ppbv−1 and 0.290 W m−2 ppbv−1, respectively, to 0.043 W m−2 ppbv−1 and 0.079 W m−2 ppbv−1, respectively for the B3LYP/6-311++G** results. Analysis by Magnitude of Radiative Efficiency. An interesting question to examine is whether computational methods are systematically better or worse for small, intermediate, and large values of the RE. To do this, three regimes were identified: RE < 0.1 W m−2 ppbv−1, 0.1 W m−2 ppbv−1 ≤ RE ≤ 0.7 W m−2 ppbv−1, and RE > 0.7 W m−2 ppbv−1. The results are presented in Table 7 where it is clear Table 7. Statistics for Calculated Values with Respect to the Magnitude of the Best Available Value of the Radiative Efficiency for the Small Basis Seta magnitude

a

n

MSE

MAE

RMS

No Vibrational Intensity Scaling RE < 0.1 31 0.014 0.1 ≤ RE ≤ 0.7 195 0.020 RE > 0.7 9 0.161

0.016 0.040 0.240

0.043 0.068 0.320

Including Vibrational Intensity Scaling RE < 0.1 31 0.008 0.1 ≤ RE ≤ 0.7 195 −0.018 RE > 0.7 9 0.009

0.013 0.042 0.208

0.036 0.061 0.253

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.est.5b04154. Sources for reference RE values are given here, in addition to RE calculated values for more than 1200 chemcial substances (PDF)

■

AUTHOR INFORMATION

Corresponding Author

*Phone: 702-798-2116; e-mail: [email protected].

Units are W m−2 ppbv−1.

Notes

The authors declare no competing financial interests.

■

that the larger values of the RE contribute strongly to the error. These results indicate that a larger absolute error is seen in the larger results, as was observed in the previous subsection. However, a larger percent error is given by the smaller results as can be seen in Figure 3. In this case, there are only nine compounds with an experimental RE greater than 0.7 W m−2 ppbv−1, and all of these compounds are highly fluorinated or perfluorinated ethers, which somewhat limits the conclusions that can be drawn from this result. In this article, the calculation of radiative efficiency (RE) using results derived from quantum chemistry methods for 235 chemical substances has been studied. It is found that a smaller basis set, 6-31G*, yields superior results to a larger basis set, 6-

ACKNOWLEDGMENTS We thank Jeff Yang and Lindsay Stanek of the U.S. EPA for thoughtful reviews of this work and Deborah Ottinger, Margaret Sheppard, Carol Hetfield also of the U.S. EPA for help in defining the scope of this work. We would also like to acknowledge the computer time at the Environmental Modeling and Visualization Laboratory provided for this project. This work is an official contribution of the United States Environmental Protection Agency and the National Institute of Standards and Technology and is not subject to copyright in the United States. The information in this G

DOI: 10.1021/acs.est.5b04154 Environ. Sci. Technol. XXXX, XXX, XXX−XXX

Article

Environmental Science & Technology

(HFEs) and hydrofluoropolyethers (HFPEs). J. Quant. Spectrosc. Radiat. Transfer 2011, 112, 1967−1977. (16) Hodnebrog, Ø.; Etminan, M.; Fuglestvedt, J. S.; Marston, G.; Myhre, G.; Nielsen, C. J.; Shine, K. P.; Wallington, T. J. Global warming potentials and radiative efficiencies of halocarbons and related compounds: a comprehensive review. Rev. Geophys. 2013, 51, 300−378. (17) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; Nakatsuji, H.; Caricato, M.; Li, X.; Hratchian, H. P.; Izmaylov, A. F.; Bloino, J.; Zheng, G.; Sonnenberg, J. L.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Vreven, T.; Montgomery, J. A., Jr.; Peralta, J. E.; Ogliaro, F.; Bearpark, M.; Heyd, J. J.; Brothers, E.; Kudin, K. N.; Staroverov, V. N.; Keith, T.; Kobayashi, R.; Normand, J.; Raghavachari, K.; Rendell, A.; Burant, J. C.; Iyengar, S. S.; Tomasi, J.; Cossi, M.; Rega, N.; Millam, J. M.; Klene, M.; Knox, J. E.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Martin, R. L.; Morokuma, K.; Zakrezewski, V. G.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Dapprich, S.; Daniels, A. D.; Farkas, O.; Foresman, J. B.; Ortiz, J. V.; Cioslowski, J.; Fox, D. J. Gaussian 09, Revision D.01, 2010. (18) Myhre, G.; Stordal, F.; Gausemel, I.; Nielsen, C. J.; Mahieu, E. Line-by-line calculations of thermal infrared radiation representative for global condition: CFC-12 as an example. J. Quant. Spectrosc. Radiat. Transfer 2006, 97, 317−331. (19) Cramer, C. J. Essentials of Computational Chemistry: Theories and Models; John Wiley & Sons, Ltd.: West Sussex, England, 2004. (20) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation made simple. Phys. Rev. Lett. 1996, 77, 77. (21) Grimme, S.; Ehrlich, S.; Goerigk, L. Effect of the damping function in dispersion corrected density functional theory. J. Comput. Chem. 2011, 32, 1456−1465. (22) Peverati, R.; Truhlar, D. G. An improved and broadly accurate local approximation to the exchange-correlation density functional: the MN12-L functional for electronic structure calculations in chemistry and physics. Phys. Chem. Chem. Phys. 2012, 14, 13171−13174. (23) Zheng, J.; Yu, T.; Papajak, E.; Alecu, I. M.; Mielke, S. L.; Truhlar, D. G. Practical methods for including torsional anharmonicity in thermochemical calculations on complex molecules: the internalcoordinate multi-structural approximation. Phys. Chem. Chem. Phys. 2011, 13, 10885−10907. (24) Irikura, K. K.; Johnson, R. D.; Kacker, R. N. Uncertainties in scaling factors for ab initio vibrational frequencies. J. Phys. Chem. A 2005, 109, 8430−8437. (25) NIST Computational Chemistry Comparison and Benchmark Database NIST Standard Reference Database Number 101 Release 16a, August 2013. http://cccbdb.nist.gov/ (accessed August 2014). (26) Halls, M. D.; Schlegel, H. B. Comparison of the performance of local gradient-corrected and hybrid density functional models in predicting infrared intensities. J. Chem. Phys. 1998, 109, 10587−10593.

document has been funded by the United States Environmental Protection Agency. It has been subjected to the Agency’s peer and administrative review and has been approved for publication. Mention of trade names or commercial products does not constitute endorsement or recommendation by EPA or NIST for use.

■

REFERENCES

(1) IPCC. 1990. Climate Change: The IPCC Scientific Assessment; Cambridge University Press: New York, 1990. (2) Manning, M.; Reisinger, A. Broader perspectives for comparing different greenhouse gases. Philos. Trans. R. Soc., A 2011, 369, 1891− 1905. (3) Pachauri, R. K.; Allen, M. R.; Barros, V. R.; Broome, J.; Cramer, W.; Christ, R.; Church, J. A.; Clarke, L.; Dahe, Q.; Dasgupta, P.; Dubash, N. K.; Edenhofer, O.; Elgizouli, I.; Field, C. B.; Foster, P.; Friedlingstein, P.; Juglestvedt, J.; Gomez-Echeverri, L.; Hallegatte, S.; Hegerl, G.; Howden, M.; Jiang, K.; Jimenez Cisneroz, B.; Kattsov, V.; Lee, H.; Mach, K. J.; Marotzke, J.; Mastrandrea, M. D.; Meyer, L.; Minx, J.; Mulugetta, Y.; O’Brien, K.; Oppenheimer, M.; Pereira, J. J.; Pichs-Madruga, R.; Plattner, G.-K.; Pörtner, H.-O.; Power, S. B.; Preston, B.; Ravindranath, N. H.; Reisinger, A.; Riahi, K.; Rusticucci, M.; Scholes, R.; Seyboth, K.; Sokona, Y.; Stavins, R.; Stocker, T. F.; Tschakert, P.; van Vuuren, D.; van Ypserle, J.-P. IPCC, 2014: Climate Change 2014 Synthesis Report. Contribution of Working Groups I, II and III to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change; IPCC: Geneva, Switzerland, 2014. (4) Young, C. J.; Hurley, M. D.; Wallington, T. J.; Mabury, S. A. Molecular structure and radiative efficiency of fluorinated ethers: A structure-activity relationship. J. Geophys. Res. 2008, 113, 113. (5) Kokkila, S. I.; Bera, P. P.; Francisco, J. S.; Lee, T. J. A group increment scheme for infrared absorption intensities of greenhouse gases. J. Mol. Struct. 2012, 1009, 89−95. (6) Kazakov, A.; McLinden, M. O.; Frenkel, M. Computational design of new refrigerant fluids based on environmental, safety, and thermodynamic characteristics. Ind. Eng. Chem. Res. 2012, 51, 12537− 12548. (7) McLinden, M. O.; Kazakov, A. F.; Brown, J. S.; Domanski, P. A. A thermodynamic analysis of refrigerants: Possibilities and tradeoffs for low-GWP refrigerants. Int. J. Refrig. 2014, 38, 80−92. (8) Pinnock, S.; Hurley, M. D.; Shine, K. P.; Wallington, T. J.; Smyth, T. J. Radiative forcing of climate by hydrochlorofluorocarbons and hydrofluorocarbons. J. Geophys. Res. 1995, 100, 23227−23238. (9) Papasavva, S.; Tai, S.; Illinger, K. H.; Kenny, J. E. Infrared intensities, atomic charges, and dipole moments in the fluoroethane series using atomic polar tensor analysis. J. Geophys. Res. 1997, 102, 13643−13650. (10) Papasavva, S.; Tai, S.; Esslinger, A.; Illinger, K. H.; Kenny, J. E. Ab initio calculations of vibrational frequencies and infrared intensities for global warming potentials of CFC substitutes: CF3CH2F (HFC134a). J. Phys. Chem. 1995, 99, 3438−3443. (11) Papasavva, S.; Illinger, K. H.; Kenny, J. E. Ab initio calculations on fluoroethanes: geometries, dipole moments, vibrational frequencies, and infrared intensities. J. Phys. Chem. 1996, 100, 10100−10110. (12) Blowers, P.; Moline, D. M.; Tetrault, K. F.; Tuchawena, S. L. Prediction of radiative forcing values of hydrofluoroethers using density functional theory results. J. Geophys. Res. 2007, 112, 112. (13) Blowers, P.; Tetrault, K. F.; Trujillo-Morehead, Y. Global warming potential predictions for hydrofluoroethers with two carbon atoms. Theor. Chem. Acc. 2008, 119, 369−381. (14) Bravo, I.; Aranda, A.; Hurley, M. D.; Marston, G.; Nutt, D. R.; Shine, K. P.; Smith, K.; Wallington, T. J. Infrared absorption spectra, radiative efficiencies, and global warming potentials of perfluorocarbons: comparison between experiment and theory. J. Geophys. Res. 2010, 115, D24317−12. (15) Bravo, I.; Marston, G.; Nutt, D. R.; Shine, K. P. Radiative efficiencies and global warming potentials using theoretically determined absorption cross-sections for several hydrofluoroethers H

DOI: 10.1021/acs.est.5b04154 Environ. Sci. Technol. XXXX, XXX, XXX−XXX