Predicting Gaseous Reaction Rates of Short Chain Chlorinated

Nov 5, 2014 - Reth , M.; Zencak , Z.; Oehme , M. First study of congener group patterns and concentrations of short- and medium-chain chlorinated para...
0 downloads 0 Views 2MB Size
Download PDF
Article pubs.acs.org/est

Predicting Gaseous Reaction Rates of Short Chain Chlorinated Paraffins with ·OH: Overcoming the Difficulty in Experimental Determination Chao Li,† Hong-Bin Xie,† Jingwen Chen,*,† Xianhai Yang,† Yifei Zhang,‡ and Xianliang Qiao† †

Key Laboratory of Industrial Ecology and Environmental Engineering (MOE), School of Environmental Science and Technology, Dalian University of Technology, Dalian 116024, China ‡ State Key Laboratory of Fine Chemicals, Dalian University of Technology, Dalian 116024, China Environ. Sci. Technol. 2014.48:13808-13816. Downloaded from pubs.acs.org by IOWA STATE UNIV on 01/10/19. For personal use only.

S Supporting Information *

ABSTRACT: Short chain chlorinated paraffins (SCCPs) are under evaluation for inclusion in the Stockholm Convention on persistent organic pollutants. However, information on their reaction rate constants with gaseous ·OH (kOH) is unavailable, limiting the evaluation of their persistence in the atmosphere. Experimental determination of kOH is confined by the unavailability of authentic chemical standards for some SCCP congeners. In this study, we evaluated and selected density functional theory (DFT) methods to predict kOH of SCCPs, by comparing the experimental kOH values of six polychlorinated alkanes (PCAs) with those calculated by the different theoretical methods. We found that the M06-2X/6-311+G(3df,2pd)//B3LYP/6-311 +G(d,p) method is time-effective and can be used to predict kOH of PCAs. Moreover, based on the calculated kOH of nine SCCPs and available experimental kOH values of 22 PCAs with low carbon chain, a quantitative structure−activity relationship (QSAR) model was developed. The molecular structural characteristics determining the ·OH reaction rate were discussed. logkOH was found to negatively correlate with the percentage of chlorine substitutions (Cl %). The DFT calculation method and the QSAR model are important alternatives to the conventional experimental determination of kOH for SCCPs, and are prospective in predicting their persistence in the atmosphere.

■

Arctic.5,8,15,21 As SCCPs represent a potential “new” category of persistent organic pollutants (POPs), several countries (e.g., the United States and Canada) and the European Union have imposed regulations on the use of SCCPs.22 Additionally, SCCPs are now under evaluation for inclusion in the Stockholm Convention on POPs. Despite the efforts to characterize the extent of SCCP contamination in the environment, little is known about their atmospheric degradation kinetics that governs the fate and persistence of SCCPs. Among the various atmospheric reactions, the reactions initiated by hydroxyl radicals (·OH) produced by atmospheric photochemical reactions, are of paramount importance. Thus, the ·OH reaction rate constant (kOH) is a key parameter for the fate assessment of pollutants in the troposphere. However, there are no kOH values available for any SCCPs. To date, several direct and indirect experimental methods have been developed to measure kOH values.23−26 However, the

INTRODUCTION Chlorinated paraffins (CPs) are commercially produced polychlorinated alkanes (PCAs) having carbon chain lengths from C10 to C30 and chlorine content from 30% to 70% by mass.1−3 Among them, short chain CPs (SCCPs) with carbon chain lengths ranging from C10 to C13 have received much attention due to their widespread occurrence and persistence in the environment,4−15 their potential for long-distance transport, 12,15 their tendency to bioaccumulate 8 and high toxicity.16,17 SCCPs show long-term toxicity to algae, aquatic invertebrates and fish at concentrations as low as 19.6, 8.9, and 3.1 μg/L, respectively.18 There is also sufficient evidence for carcinogenicity to humans from SCCPs of C12 carbon chain length group with an average chlorine content of 60%.17 The vapor pressures of SCCPs are in the range of other chlorinated organics with the same molecular weight range such as polychlorinated biphenyls and toxaphenes, suggesting that they have the potential to undergo long-range atmospheric transport through the grasshopper effect.4 SCCPs have been detected in environmental media including air,10,11,13 water,19 and sediment;9,12 in humans20 and organisms like fish,14 birds,15 marine mammals;21 and even in remote areas like Canadian Arctic lakes, the Hazen lake and the European © 2014 American Chemical Society

Received: Revised: Accepted: Published: 13808

January 7, 2014 November 2, 2014 November 5, 2014 November 5, 2014 dx.doi.org/10.1021/es504339r | Environ. Sci. Technol. 2014, 48, 13808−13816

Environmental Science & Technology

Article

Figure 1. Atom numbers of nine SCCPs (Dark gray: C; Light gray: H; Green: Cl).

It was the purpose of this study to find a solution for kOH prediction of SCCPs. By employing 6 molecules (CH3CCl3, CH2ClCHCl2, CH2ClCHClCH3, CH3(CHCl)2CH3, CH2Cl(CH2)3CH3 and CH2Cl(CH2)4CH3) with available experimental kOH values as test cases, we compared the prediction accuracy of several DFT methods. The expensive ab initio second-order Moeller−Plesset (MP2)37 method, which consistently describes all types of correlation energy including intermolecular correlation and intramolecular correlation terms, and therefore is reliable for geometry optimization,38−40 was also adopted as a reference. Based on the established method, kOH values for nine SCCPs with carbon chain ranging from C10 to C13 (Figure 1), were calculated. A QSAR model for predicting kOH of more SCCPs was developed based on the DFT calculated and available experimental kOH values for PCAs. To the best of our knowledge, this is the first study concerning kOH of SCCPs. The DFT calculation method and the QSAR model can be employed to predict kOH of other SCCPs.

experimental determination of kOH for SCCPs is difficult as the quantities and kinds of authentic chemical standards for some SCCPs, which are necessary for the determination, are very limited. Moreover, there are a large number of isomers for SCCPs,22 which could further increase the workloads and difficulties encountered in the experimental determination. Another solution is to predict kOH based on quantitative structure−activity relationship (QSAR) models, on which the Organization for Economic Co-operation and Development (OECD) has issued guidelines on model development and validation.27 Nevertheless, the utility of QSAR models is constrained by the applicability domain of the models.28,29 The currently available QSAR models on kOH do not cover SCCPs. For example, two latest QSAR models on kOH prediction cover only a limited number of PCAs with chain lengths ranging from C1 to C6.30,31 Thus, seeking a new solution for predicting kOH of SCCPs is necessary. In recent years, the increasing computational power has enabled quantum chemical calculations for molecules or systems comprising up to a few hundred atoms.32 Some studies in recent years showed that the density functional theory (DFT) can accurately predict kinetics, pathways and products for the reactions of chemicals with ·OH.33−35 For example, Zhou et al. employed DFT to calculate the reaction of 4,4′-dibromodiphenyl ether with ·OH, and found that the predicted products and the overall rate constants (kOH) at 298 K agreed well with the experimental results.33 Ci et al. employed DFT to compute the hydrogen abstraction reactions of CF3CH2CHO with ·OH, also found that the calculated kOH values were consistent with the experimental data.35 Furthermore, the theoretical approaches have the advantage that they can provide details involving reaction sites and favorable reaction pathways.33−36 As far as we know, there are no previous quantum chemical calculations on reactions of SCCPs with ·OH.

■

COMPUTATIONAL DETAILS Global Minimum Search. As the molecular structures of PCAs consist of sp3-hybridized C-X (X = H, Cl) bonds, they have a range of conformations. Thus, the global minimum search is needed to identify the optimal structures of the reactants under study, and the optimal structures were then included in further structure calculations followed by dynamic calculations of kOH. The calculations were carried out using TURBOMOLE.41 Ab initio molecular dynamics (AIMD) was used to generate reasonable gas-phase geometries by employing the BLYP functional42 along with a triple-ξ valence polarized basis set (TZVP) within the resolution-of-the-identity (RI) approximation.41 In order to sample enough conformations in a short time and increase to some extent the possibility for getting the real global minimum, the temperature was set to 700 K with the 13809

dx.doi.org/10.1021/es504339r | Environ. Sci. Technol. 2014, 48, 13808−13816

Environmental Science & Technology

Article

Figure 2. Calculated and experimental kOH values for CH3CCl3 over different temperature ranges (a), and the enlarged view of (a) at T = 200−330 K is shown in (b).

QSAR Modeling. The DFT calculated kOH values (T = 298 K) of 9 SCCPs, together with those of 22 low carbon chain PCAs (C1 ∼ C6) with available experimental kOH values,31 were employed for QSAR modeling. According to previous studies,30,31,52−55 the following quantum chemical descriptors were selected to construct the model, including the atomic Mulliken charges [the most negative net C atom charge (QC−), the most positive net H atom charge and the most negative net Cl-atom charge], polarizability (α), the energy of the highest occupied molecular orbital (EHOMO), the energy of the lowest unoccupied molecular orbital (ELUMO), chemical potential (μ = (ELUMO + EHOMO)/2),56 and absolute hardness (η = (ELUMO − EHOMO)/2).56 These descriptors were calculated from the geometries optimized by the B3LYP/6-311+G(d,p) method using the Gaussian 09 program suite. Additionally, two descriptors that describe the chlorine substitution characteristics [the number of chlorine atoms (nCl) and percentage of chlorine substitutions (Cl%)], 3D-MoRSE (3D-molecular representation of structure based on electron diffraction) and geometrical descriptors were also considered, for which the DRAGON software was employed for the calculation.57 For the modeling, the data were randomly split into a training set and a validation set with a ratio of 4:1. Stepwise multiple linear regression (MLR) analysis was employed to construct the QSAR model. The applicability domain of the QSAR model was assessed by the Euclidean distance based method, and was conducted using the AMBIT Discovery (version 0.04).58

premise that the molecules would not dissociate in the AIMD calculation (detailed in the SI). A time step of 0.96752 fs was used, and the total length of each AIMD simulation was 20 000 time steps. We selected the conformations from the AIMD run as the starting point for geometry optimization at the B3LYP/ 6-311+G(d,p) level using the Gaussian 09 program suite.43 The minimum with the lowest energy was identified as the global minimum used to investigate its reaction with ·OH. Electronic Structure Calculations. A gradual scheme was used to screen reliable theoretical methods for the geometry optimization and frequency analysis of the studied systems. For all the systems, the single point energy calculation was performed at the M06-2X/6-311+G(3df,2pd) level, since the M06-2X functional gives the best performance (mean signed errors of −0.51 kcal/mol) without increasing the computational time for the hydrogen-transfer barrier height calculation.44 As CH3CCl3 has a C3v symmetry that can reduce the computational time, its reaction with ·OH was taken as an initial test to screen the different computational methods (the functional BHandHLYP,42,45 MPW1K,46 B3LYP,42,47 TPSSh,48 B97−149 and ab initio MP237 in conjunction with the 6-311+G(d,p) basis set) by comparing the calculated zero-point corrected reaction barriers (Ea). Subsequently, the selected time-effective and reliable methods (B3LYP, TPSSh and B97-1 functionals) were evaluated with CH2ClCHCl2 and CH2ClCHClCH3 as test cases, followed by using CH3(CHCl)2CH3, CH2Cl(CH2)3CH3 and CH2Cl(CH2)4CH3 for further method validation. Finally, an optimum method (B3LYP/6-311+G(d,p)) was selected for the SCCPs. All the stationary points including reactants, products, precomplexes and product-like complexes have real frequencies, whereas the transition states (TSs) have only one imaginary frequency. The minimum energy paths (MEP) were obtained by the intrinsic reaction coordinate (IRC) theory for each reaction channel.50 Dual-Level Direct Dynamic Calculation. By means of the POLYRATE 2010-A program,51 the rate constant for each reaction channel (ki) was calculated using the improved canonical variational transition-state theory (ICVT) with small-curvature tunneling (SCT) correction51 over the temperature range 200−500 K. The frequencies of the selected points along the MEP were calculated at the B3LYP/6-311+G(d,p) level to obtain the Cartesian coordinates, gradient and Hessian matrix; and the single-point energy of these points was calculated by the M06-2X/6-311+G(3df,2pd) method. kOH was calculated by summing up the ki values of different channels.

■

RESULTS AND DISCUSSION Evaluation of the Different Calculation Methods. The reaction of PCAs with ·OH proceeds via H atom abstraction, as the abstraction of Cl atoms has much higher Ea and is endothermic.59 For the reaction CH3CCl3+·OH, we considered only one H-abstraction process from −CH3, as the three H atoms in −CH3 are equivalent (Supporting Information (SI) Figure S1). The Ea values calculated from the different methods, together with the experimental values are listed in SI Table S1. We found that the Ea values from the geometries optimized by the TPSSh, B3LYP and B97−1 functionals are close to those calculated by the more-reliable MP2/6311+G(d,p) method, and are closer to the experimental values25,60−63 than the other functionals under investigation. Thus, the TPSSh, B3LYP, and B97-1 functionals were employed for geometry optimization and frequency analysis of the other two tested PCAs. 13810

dx.doi.org/10.1021/es504339r | Environ. Sci. Technol. 2014, 48, 13808−13816

Environmental Science & Technology

Article

Figure 3. Calculated and experimental kOH values for CH2ClCHCl2 over different temperature ranges (a), and the enlarged view of (a) at T = 200− 330 K is shown in (b).

Figure 4. Calculated and experimental kOH values for CH2ClCHClCH3 over different temperature ranges (a), and the enlarged view of (a) at T = 200−320 K is shown in (b).

more universal than the B97-1 functional for predicting kOH in both the low and high temperature ranges. Furthermore, Wang et al. calculated the kOH of CH3CH2CH2Cl and CH3CHClCH3, and also concluded that the B3LYP/6-311G(d,p) method is reliable for geometries optimization.68 As the studied system involves radicals and Cl atoms, inclusion of a diffuse function to the basis set should increase the accuracy of kOH prediction. Thus, we selected the M06-2X/6-311+G(3df,2pd)//B3LYP/6311+G(d,p) method for predicting kOH values of SCCPs. The selected method was further verified with the 3 PCAs with C4−C6 carbon chain lengths, and the results are listed in SI Table S2−S5. The deviations between the calculated kOH values and the corresponding experimental values are within a factor of 0.6−1.7, 0.4−0.7, 0.6−0.7, 0.5, 0.4−0.6, and 0.5−0.6 for CH3CCl3, CH2ClCHCl2, CH2ClCHClCH3, CH3(CHCl)2CH3, CH2Cl(CH2)3CH3, and CH2Cl(CH2)4CH3, respectively, which are generally acceptable.69,70 Thus, the deviation of the DFT calculated kOH values of the PCAs is estimated to be within a factor of 0.4−1.7. Reactions of the SCCPs with ·OH. For convenience, all the atoms of the 9 SCCPs are numbered (Figure 1). Due to the C1 symmetry of the SCCPs, all the channels of H-abstraction for the SCCPs were considered. Similar to the 3 PCAs with low carbon chains, all the H-abstraction channels proceed via precomplexes, TSs, product-like complexes and products (SI Figure S3). The calculated activation free energy (ΔG) and Ea for each channel are listed in SI Tables S6−S14. By comparing the Ea values of the different H-abstraction channels, we found the channels with six-membered ring TSs formed by H-bonds either between H and O or between H and Cl atoms, tend to have low Ea and ΔG values. Generally, the conformations with cis-H-C-C-Cl substructural units facilitate the formation of the six-membered ring TSs. Thus, chlorine

For CH2ClCHCl2 and CH2ClCHClCH3, the possible reaction pathways for the ·OH reaction are described in SI Figure S1. The potential energy surface profiles for the three tested PCAs above are shown in SI Figure S2. As H-abstraction by ·OH can occur from −CHCl2 and CH2Cl−, CH2ClCHCl2 has three TSs (TS 1b , TS 2b , and TS 3b ). Similarly, CH2ClCHClCH3 has six TSs (TS1c TS2c, TS3c, TS4c, TS5c, and TS6c). The calculated and experimental25,26,60−67 kOH values for the three PCAs are shown in Figures 2−4. Several experimental studies reported the kOH values of CH3CCl3,25,60−66 which makes it possible to estimate the deviation of the experimental determinations. At 298 K and at 95% confidence level, the deviation of experimental kOH values of CH3CCl3 is within a factor of 0.4−2.8 (detailed in the SI). As for CH2ClCHCl2, only two studies reported the kOH values,25,26 and the biggest ratio between them is 1.7 (T = 295 K). By comparing the calculated with the available experimental values (Figures 2-4), it can be concluded that the TPSSh functional overestimates the kOH values of CH3CCl3 in the range T = 200−720 K except for the experimental kOH value of Chang et al., 61 and underestimates k OH values of CH2ClCHClCH3 in the range T = 200−372 K. Thus, the TPSSh functional was excluded for its inability for accurately predicting kOH of the PCAs. There are no significant differences between the kOH values calculated from the B97−1 and B3LYP functionals in the low temperature range (200−330 K). However, in the range of T > 340 K, the kOH values calculated from the B3LYP functional are generally more consistent with the experimental values than those calculated with the B97−1 functional. The kOH data in the high temperature range are important in understanding the incineration process of PCAs.25,65 Thus, the B3LYP method is 13811

dx.doi.org/10.1021/es504339r | Environ. Sci. Technol. 2014, 48, 13808−13816

Environmental Science & Technology

Article

coefficient (Q2ext) is 0.856, and the RMSEext is 0.264, indicating the model has good predictive ability. The plot of the QSAR predicted versus experimental or DFT calculated kOH values is displayed in Figure 6. The applicability

substitutional positions and steric structures of SCCPs have great effects on the reaction rates with ·OH. Beyond that, there seem no regularities for the H-abstraction channels. Thus, it is necessary to take every H atom into account when calculating the kOH values. The DFT calculated kOH values for the 9 SCCPs over the temperature range 200−500 K are shown in Figure 5 and SI

Figure 6. Plot of the QSAR predicted versus experimental or DFT calculated logkOH values for the training and validation sets. Figure 5. Plot of the DFT calculated kOH values for the 9 SCCPs versus temperature.

domain of the model is described in SI Figures S4 and S5, which show that all of the compounds are in the domain and there are no outliers for both the training and validation sets. SPH, a descriptor characterizing spherosity of organic molecules, is the most important descriptor in the model (t = 16.6 and p < 0.0001). SPH positively correlates with logkOH and explains 80.9% variance of logkOH. The shapes of organic molecules with high SPH values tend to be more spherical. Additionally, we found that SPH negatively correlates with Cl% (p < 0.001, SI Figure S6), further indicating that the chlorine substitution and steric structures of SCCPs determine their reaction rates with ·OH. As far as we know, this is the first report that SPH is an important predictor variable for logkOH prediction of PCAs. QC− represents the most negative C atom charge in an organic molecule. With a more negative QC− charge, the H atom connected to the C atom tends to have a more positive charge, leading to low reactivity toward ·OH, as · OH is a strong electrophile. The absolute hardness η measures the stability of a system. A hard molecule resists changes within itself, or in reaction with others.74 Thus, it is reasonable that PCAs with high η values tend to have low logkOH values. Implications. Based on the global average ·OH concentration75 of 9.7 × 105 molecules cm−3 and the DFT calculated kOH values, lifetimes (τ = 1/(kOH[OH]) of the nine SCCPs were calculated (SI Table S18). τ of SCCP-3 and SCCP-9 are 2 days. According to the deviation of the DFT calculated kOH values (a factor of 0.4−1.7), the deviation of the τ values is within a factor of 0.6−2.5. To understand the influence of chlorine substitution on τ values, the QSAR model was employed to predict logkOH and τ values of some other SCCPs (SI Table S19). As can be seen from SI Figure S7, the τ values of the SCCPs positively correlate with the percentage of chlorine atoms Cl% (p < 0.001), implying that SCCPs with more chlorine substitutions tend to be more widespread and persistent in the atmosphere. Figure 7 shows the variation of the τ values with the number of chlorine substitution (nCl) for the SCCPs with different carbon chain lengths. It can be observed from Figure 7 that the C10−11Cl5−8 and C12Cl6−8 groups with nCl ≤8 tend to have long lifetimes, and the C10−13 groups (especially the C13 group) with

Table S15. In general, the kOH values display a negative temperature dependence in the low temperature range. The reason can be ascribed to the negative Ea of the reactions (SI Tables S6−S14). In principle, for reaction pathways with negative Ea, the reaction rates at low temperatures are mainly controlled by the kinetics, and at high temperatures are mainly affected by the thermodynamics.36 The similar temperature dependence of kOH was also reported in the reaction of other compounds with ·OH.33,35,36,71 Driven by curiosity, we compared the DFT calculated kOH values at 298 K with those predicted by QSAR models, although the training sets of the QSAR models did not cover SCCPs. The QSAR predicted kOH values are all lower than the DFT calculated values. The deviations between the DFT calculated kOH values and those predicted by the QSAR models of Li et al.,31 Wang et al.,30 Gramatic et al.53 and Roy et al.,55 are within a factor of 3.0−6.8, 2.8−22.5, 1.2−15.9, and 40.5− 397.4, respectively. Thus, it is necessary to develop a new QSAR model that covers SCCPs in the applicability domain. QSAR Modeling. The following optimum QSAR model for logkOH prediction (T = 298 K) was obtained: log k OH = 5.057SPH + 1.167Q C− − 13.991η − 12.186

The compounds and the descriptor values involved in the QSAR model are listed in SI Table S16. For this model, there are 25 PCAs (seven SCCPs) included in the training set. The variable inflation factors (VIF) for the predictor variables are all