Environ. Sci. Technol. 2006, 40, 2298-2303
Molecular Simulations of Benzene and PAH Interactions with Soot JAMES D. KUBICKI† Department of Geosciences and The Earth and Environmental Systems Institute, The Pennsylvania State University, University Park, Pennsylvania 16802
Molecular mechanics simulations and ab initio calculations were performed to investigate the mechanism of PAHsoot adsorption. Partitioning of benzene, naphthalene, fluorene, phenanthrene, anthracene, pyrene, fluoranthene, benzo[a]anthracene, benzo[k]fluoranthene, benzo[a]pyrene, and benzo[g,h,i]perylene between water and soot was modeled with classical and quantum mechanical calculations in order to determine a method for predicting log(Kd) values. In both cases, the predicted mechanism of adsorption is interaction of the π-electrons in the PAH and soot (i.e., π-π van der Waals forces). Solvation energies, the energy difference between the solute in the gas phase and in the model aqueous phase, calculated with molecular mechanics did not follow the observed solubilities of the PAHs. Molecular dynamics simulations overestimate the favorability of PAHs in the aqueous phase. Hence, the partitioning between the aqueous phase and soot does not accurately correlate with observed log(Kd) values. Models of PAH adsorption using structures from molecular mechanics and energies from ab initio calculations do produce water-soot partitioning energies that correlate well with observed log(Kd) values. The log(Kd) values for benzene, anthracene, fluorene, and benzo[a]pyrene were predicted based on the correlation between calculated partitioning energies and observed log(Kd) values. Results presented here suggest that partitioning of PAHs onto soot should depend on the size of the PAH, the planarity of PAH molecule, and the aromaticity of the compound. The methodology developed by this research can be used to predict PAH Kd values that have not yet been measured.
Introduction The observation that particles of soot or black carbon can influence the transport, bioavailability, and fate of organic contaminants (1), such as polycyclic aromatic hydrocarbons (PAHs), has led to research on the interactions of these two components of soils and sediments (2-5). Although difficult to detect in nature due to their small grain size, amorphous character, and indistinct composition, soot particles can be the dominant phase containing PAHs, especially when the total organic carbon (TOC) of a soil or sediment is low (6). Bucheli and Gustafsson (2) conducted a study in which a series of PAHs (naphthalene, fluoranthene, phenanthrene, and pyrene) were adsorbed to a standard diesel soot (NIST SRM 1650) to determine distribution coefficients, log(Kd), of these PAHs between soot and water. The research presented † Phone: (814)865-3951; fax: (814)863-7823; e-mail: kubicki@ geosc.psu.edu.
2298
9
ENVIRONMENTAL SCIENCE & TECHNOLOGY / VOL. 40, NO. 7, 2006
here is motivated by these experiments to determine PAH log(Kd) values. Two purposes are served by performing molecular simulations that parallel the experimental adsorption studies. First, if the molecular simulations can mimic the observed behavior of the experimental system, then the simulation methods used acquire a significant level of validation. With a validated simulation method, the molecular-level interactions of PAHs and soot can be studied. This issue is important because there are many different types of “black carbon” in the environment (e.g., wood soot, diesel soot, charcoalssee Jonker and Koelmans (3) for experiments on partitioning to various types of black carbon). If each type has significantly different log(Kd) values, then the black carbon at each site will need to be characterized and the corresponding adsorption isotherm determined. Alternatively, if PAH partitioning to black carbon is controlled by similar chemical interactions in each case, then generalized log(Kd) values may be employed to predict long-term behavior. For example, Jonker and Koelmans (3) measured PAH sorbent-water distribution coefficients as log(Ks) values for a variety of black carbon types. Phenanthrene has values of 6.24, 5.27, 5.31, 5.54, 6.57, 6.19, and 5.68 for traffic soot, oil soot, wood soot, coal soot, coal, charcoal, and diesel soot. Thus, excluding coal and charcoal, which are not classified as “soot” because they do not form in the gas phase, the measured range is 5.27-6.24. This is a remarkably narrow range. The differences should be examined further to look for compositional or structural differences that can explain this variation of 1 log unit, but a similar mechanism is likely to be responsible for partitioning in each case. Note that the models created for this study do not attempt to mimic charcoal, coal, or free-phase oils. Partitioning from these phases will be different than from soot (3, 4). Furthermore, the modeling presented here is applicable to samples where the concentration of PAHs are well below aqueous solubility. When PAH concentrations are high, a free-phase oil may form (4), and the partitioning behavior of PAHs onto soot cannot be distinguished from that of sedimentary organic matter (5). A related issue is the mechanism of partitioning versus adsorption. PAHs measured in sediments are commonly thought to be equivalent, but often biodegradation experiments have revealed components that are highly resistant to decomposition even though they are labile in the free state (7, 8). Consequently, the question of PAH sources can be raised. If the PAHs partition onto soot from the aqueous phase (e.g., from oil spills), then they may be located at the surface and relatively available. On the other hand, if the PAHs are coformed with the soot during combustion (9, 10), they may be located inside the soot particle and highly recalcitrant, e.g., diesel exhaust washed into storm drains (11). However, even adsorbed PAHs on the surface of a soot are likely to be resistant to biodegradation because the sootwater partition coefficients are high (2, 3). A second purpose of performing molecular modeling studies of PAH-soot interactions is that determining distribution coefficients via experiment may require as long as 3 months in the laboratory (2). In contrast, the partitioning energies predicted in this research require on the order of 1 day of computer time to generate. Given the large number compounds that could potentially require log(Kd) values, obtaining a faster route to predicting PAH partitioning onto soot could save time and money. Ideally, one would like to calculate ∆G for adsorption because ln(K) ) -∆G/RT in general. However, calculations 10.1021/es051083s CCC: $33.50
2006 American Chemical Society Published on Web 03/04/2006
of ∆G using computational chemistry techniques is highly computationally demanding for one system because a statistically significant sampling of conformation space must be sampled in order to calculate the Helmholtz or Gibbs free energies. Performing these calculations on a large variety of different compounds would be impractical using quantum mechanical approaches. Continuum approaches, such as those applied to estimating ∆G for solvation, are not likely to be appropriate for modeling adsorption where half of the molecule is in contact with the solvent and the other half in contact with the substrate. Hence, for modeling PAH-soot interactions one is faced with a choice: use molecular mechanics methods that are quick but perhaps inaccurate or use quantum mechanics and make approximations about the adsorption behavior of PAHs to soot. In this paper, both approaches are investigated. The first part of this paper discusses an attempt to use molecular mechanics in molecular dynamics simulations in order to model the energetics of PAH adsorption to soot. Although the force field chosen does a reasonably accurate job of reproducing organic structures, the solvation and adsorption energetics calculated here are qualitatively wrong. This is not to say that no force field will ever be able to be used for this task, but it is a warning against using a deceptively simple route to determining the molecular behavior of PAHs and other contaminants in the environment. The second half of the paper describes a quantum mechanical method for predicting ln(Kd) values. Instead of attempting to calculate ∆G, the thermodynamic property ∆H is estimated instead. The justification of this approach can be seen by examining the following basic thermodynamic relationships. To begin, ∆G ) ∆H - T∆S and ∆H ) ∆E + P∆V. The assumption here is that the P∆V term is negligible at 1 atm pressure for these adsorption reactions; hence ∆H ≈ ∆E. Also, ln(K) ) -∆G/RT ) -∆H/RT + ∆S/R. When one plots ∆E (as an approximation to ∆H) versus ln(K), the slope should be equal to (-1/RT) and the intercept equal to ∆S/R. This is why the computed ∆E’s are plotted against ln(Kd). The assumption is that ∆S is relatively constant over this range of adsorption reactions. This assumption is not justified a priori but will be justifiable if the correlation is linear over a significant range of ln(Kd) values. If ∆S varies significantly, then the correlation would deviate from linearity. Partition energies, ∆EPart, were then determined using the computed average potential energies for each model system via
∆EAds ) E(Soot + PAH) - E(Soot) - E(PAH) ∆ESolv ) E(Water + PAH) - E(Water) - E(PAH) ∆EPart ) ∆EAds - ∆ESolv This formulation takes into account the solvation energies of the compounds as well as the energy of interaction with the substrate to determine partitioning behavior. It is critical to remember that it is the relative energy difference between the dissolved and adsorbed state that determines partitioning, not just the energy of interaction between the substrate and adsorbate.
Methods Molecular Mechanics. The molecular modeling package Cerius2 (12) was employed for simulations of larger scale models. Potential energies were derived from molecular dynamics (MD) simulations that compute the position and velocity of all particles in the system as a function of time. A temperature of 300 K (27 °C) set by the Hoover thermostat, and a time step of 1 × 10-15 s was used for all simulations.
TABLE 1. Average Potential Energies and Standard Deviations in Parentheses from COMPASS-Based MD Simulations Used to Calculate Solvation, Adsorption, and Partition Energies for Benzene and PAHs Partitioning between Water and Soot model soot water benzene benzene + water benzene + soot naphthalene naphthalene + water naphthalene + soot fluorene fluorene + water fluorene + soot anthracene anthracene + water anthracene + soot phenanthrene phenanthrene + water phenanthrene + soot pyrene pyrene + water pyrene + soot fluoranthene fluoranthene + water fluoranthene + soot
potential energy (kcal/mol) ∆ESolv ∆EAds ∆EPart 8155(17) -1777(18) 12(2) -1778(19) 8151(18) 123(4) -1640(19) 8252(18) 135(5) -1652(16) 8265(16) 234(3) -1525(20) 8368(17) 235(3) -1539(19) 8358(15) 338(5) -1450(21) 8477(20) 356(8) -1412(17) 8483(18)
-3 -13 -16 -30 +14 -26 -15 -10 -25 -39 +18 -21 -35 +3 -32 -5 -11 -16 -37 +9 -28
The force field COMPASS (13) was used to calculate atomic interactions. Charges were represented by the default atomic charges in COMPASS. Simulations were performed within 3-D periodic boundary conditions. Three model types were generated for each compound studied: an isolated gas-phase molecule, the molecule surrounded by 207 H2O molecules (at a density of approximately 1 g/cm3), and the molecule in juxtaposition with a soot model (C338H84O30) based on the hexane soot structure proposed by Akhter et al. (14, 15). The water and soot model systems were also simulated without a PAH present to obtain the potential energies of these systems. Simulations were then run for 20 000 time steps and the potential energies averaged over the last 10 000 time steps of the simulation. Ten picoseconds is considered a relatively short MD simulation; however, comparison of longer runs showed that the potential energies in these models were determined within this time frame. For example, a 100 ps simulation of a phenanthrene-water model system was run to test the effect of time on calculated energies. The longer run resulted in an average potential energy of -1563 ( 20 kcal/mol for steps 30-100 ps compared to -1562 ( 20 kcal/mol for steps 1020 ps (Table 1), which was considered an insignificant change. Quantum Mechanics. Molecular orbital calculations were applied to smaller model systems. The “soot” in this case was a single coronene (C24H12) molecule. This simplification is justifiable if the soot-PAH adsorption occurs mainly through π-π interactions. The scheme above was used to calculate adsorption, solvation, and partitioning energies with quantum mechanical calculations as well, but the details of computing the component terms (e.g., E(Soot + PAH)) were altered. The solvation energy was computed as the difference between the isolated gas-phase PAH molecule and the molecule enclosed in a dielectric continuum with a dielectric constant of 78.4 (16). Cavitation, dispersion, and repulsion energies were computed with Gaussian 03 (17) and added to the solvation energy to account for the loss of H-bonding between H2O molecules as the solute displaces water in the volume it occupies (18). (Note: this factor is explicitly included in VOL. 40, NO. 7, 2006 / ENVIRONMENTAL SCIENCE & TECHNOLOGY
9
2299
TABLE 2. Partitioning Energies Required to Move Each Compound from Octanol into Water Calculated with the MP2/6-31G(d) Method and the Dielectric Continuum Solvation Model IEFPCM (16)a
benzene naphthalene fluorene anthracene phenanthrene pyrene fluoranthene BaP BaA BkF BghiPe
EGas hartree
EAq hartree
EOct hartree
∆ECav kJ/mol
∆EOct/Water kJ/mol
-231.4572 -384.6136 -499.7716 -537.7637 -537.7732 -613.7824 -613.7161 -766.9432 -690.9264 -766.9070 -842.9153
-231.4632 -384.6221 -499.7812 -537.7747 -537.7842 -613.7936 -613.7279 -766.9515 -690.9398 -766.9212 -842.9280
-231.4622 -384.6207 -499.7796 -537.7729 -537.7824 -613.7917 -613.7258 -766.9491 -690.9373 -766.9187 -842.9257
+13 +18 +21 +27 +23 +24 +24 +29 +28 +30 +29
+10 +14 +18 +22 +18 +19 +18 +23 +23 +23 +23
a E Aq and EOct (including electrostatic solute-solvent, cavitation, dispersion, and repulsion energies) calculated with Gaussian 03, ref 14. ∆EOct/Water ) (EAq - EOct) × 2625.5 kJ/hartree. Correlation with the observed ln(Kow) values (2444) is shown in Figure 2a. BaP ) benzo[a]pyrene, BaA ) benzo[a]anthracene, BkF ) benzo[k]fluoranthene, and BghiPe ) benzo[g,h,i]perylene.
the MD simulations of PAHs in water because H2O-H2O H-bonds are disrupted by the presence of the solute.) Structures for each PAH-soot dimer were obtained via energy minimizations using the COMPASS force field (13) in Cerius2. Second-order Møller-Plesset, MP2, (19) calculations with the 6-31G(d) basis set (20) were used to obtain energies. Basis set superposition errors (that account for energy gained via juxtaposition of orbitals from a second molecule) were estimated by the counterpoise methods of Boys and Bernardi (21).
Results and Discussion Results of PAH-soot partitioning energies derived from MD simulations using the COMPASS force field have been reported previously (22). These results suffered from large fluctuations in the potential energy due to the fact periodic boundary conditions were not included in these simulations. The surface effects on the water nanodroplets were too large to allow derivation of a significant correlation between model energies and experimental partitioning. The calculations discussed above attempted to solve this problem by including periodic boundary conditions, thus eliminating surface effects and decreasing the standard deviation in the computed potential energy. This change in method decreased the computed potential energy fluctuations. The typical standard deviation of potential energy in the 3-D periodic simulations was on the order of 80 kJ/mol compared to 4000 kJ/mol for the simulations without boundary conditions (22). Considering the potential energy fluctuations should be smaller in the nanodroplet model based purely on the number of atoms present (≈650 periodic vs ≈3100 nonperiodic), this difference clearly demonstrates the need for 3-D periodic boundary conditions in assessing thermodynamic parameters from the MD simulations. Although the 3-D periodic MD simulations resulted in smaller standard deviations, no correlation resulted between computed ∆EPart and the experimental log(Kd) values predicted with the classical MD simulations using the COMPASS force field (Table 1). The problem is due to an inadequate description of the H2O-PAH interactions. Fluorene and pyrene are predicted to have negative ∆ESolv (Table 1) even though these are hydrophobic compounds. The relative trends among the compounds do not reproduce observations either. For example, the difference between ∆EAds and ∆ESolv (i.e., ∆EPart) for pyrene is -5 kJ/mol, whereas for naphthalene it is -30 kJ/mol (Table 1). The log(Kd) for pyrene is 2 log units higher than that of naphthalene (2), so the order of these computed results is reversed. Analysis of the simulations suggests that the force field is overestimating H2O-π interactions and underestimating H2O-H2O H-bonding. 2300
9
ENVIRONMENTAL SCIENCE & TECHNOLOGY / VOL. 40, NO. 7, 2006
FIGURE 1. Phenanthrene-coronene model as an example of MP2/ 6-31G(d) calculations of PAH-soot adsorption energy. This failure prompted a search for another method to predict log(Kd) values. Previous work with density functional theory (DFT) methods did not produce a minimum energy configuration for PAH adsorption to soot (22) because DFT does not adequately describe dispersion forces (23). Dispersion forces are probably dominant in bonding nonpolar organic contaminants to soot. Consequently, the more computationally intensive method MP2 was employed. This method for describing electron correlation effects is more effective at representing dispersion forces and should be relatively accurate in this case. Because MP2 calculations are so time-consuming, a dielectric continuum model for solvation (16) was chosen rather than explicit solvation by H2O molecules, and a model PAH-coronene system (Figure 1) was used to calculate the energy of adsorption. Coronene was substituted as a truncated representation of the soot surface with the assumption the π-π van der Waals forces dominate PAH/soot attraction. A test of the MP2/6-31G(d) level of theory and the IEFPCM solvation model was conducted to verify that it would produce realistic values for distribution coefficients. The distribution coefficients between water and octanol, Kow, are widely used in environmental chemistry, so this experimental data was used as the basis for testing the methodology used here (24). However, it should be noted that there can be significant uncertainties in measured ln(Kow) values (25). Although there is some scatter, Figure 2a shows that the calculated energy differences, ∆EOct-Water (Table 2), correlate strongly with the observed ln(Kow) values for the eight compounds included in this soot partitioning study. A simpler method, calculating the surface area of the compound, also produces an excellent correlation for a suite of PAH compounds (Figure 2b). With the exception of dibenzoanthracene, an extremely quick calculation of molecular surface area can be used to predict
FIGURE 3. Calculated changes in energy for PAHs in the aqueous phase adsorbing onto the surface of a model soot are plotted against experimental log(Kd) values ((2), upper line; (3), lower line). Four experimental points for naphthalene, fluorene, phenanthrene, and pyrene were used to generate the correlation vs the (2) data. The points for benzene, anthracene, fluoranthene, and benzo[a]pyrene were generated by predicting the log(Kd) values from the theoretical soot partitioning energies and the correlation line based on the four experimental values. The lower line is based on the log(Ks) values supplied in (3) for “added PAHs” to diesel soot (SRM 1650). Calculated ∆EPart values include BSSE. Calculated log(Kd) values are 2.37 (benzene), 4.60 (naphthalene), 5.80 (fluorene), 6.37 (anthracene), 6.30 (phenanthrene), 6.90 (pyrene), 6.89 (fluoranthene), 7.56 (benzo[a]anthracene), 8.18 (benzo[k]fluoranthene), 8.21 (benzo[a]pyrene), and 8.51 (benzo[g,h,i]perylene).
FIGURE 2. (a) Correlation of the calculated energy difference of the eight compounds included in this study between the aqueous and octanol phases. More positive energy values indicate a stronger preference for the octanol phase (i.e., more hydrophobic compounds have a larger ∆EOct-Water and ln(Kow) (24). (b) Correlation between the calculated surface areas of a suite of compounds (benzene, naphthalene, acenaphthylene, acenaphthene, fluorene, anthracene, phenanthrene, fluoranthene, pyrene, chrysene, benzo[a]pyrene, dibenzoanthracene, benzoperylene, and indenopyrene) and observed ln(Kow) values (24). The surface areas were calculated with the PM3-SM3 method in AMSOL 5.4 (28). ln(Kow) values based on this correlation. Such a simple approach may be useful for estimating Kow’s and for suggesting when experimental data may be inconsistent, such as the dibenzoanthracene point in Figure 2b. However, the more computationally demanding approach of calculating ∆E values will be more widely applicable to compounds that have more specific interactions with solvent and substrate than PAH partitioning between octanol and water. With this in mind, the ∆EPart was correlated with the observed log(Kd) values taken from (2). (Note that for Kow, ln(Kow) was used because thermodynamics suggests this form based on K ) exp(-∆G/RT). However, Bucheli and Gustasfsson (2) use log(Kd), so these values are used directly here. A simple conversion can change the log(Kd) values to ln(Kd) and vice versa.) As can be seen in Figure 3, the difference between the MP2/6-31G(d) adsorption and solvation free energies (Table 3) correlates strongly with the experimental log(Kd) values of (2). The ∆EPart values include BSSE. This correction is necessary to account for the artificial energy decrease generated as basis set overlap allows extra relaxation
of the electron density that was not possible for the isolated molecules (26, 27). Boys and Bernardi (21) developed a method to estimate this BSSE called the counterpoise method, and it was included here. When BSSE was not included, however, the ability to generate a strong correlation between observed log(Kd) and calculated ∆EPart values was not affected, and the predicted log(Kd) values differed by less than 0.3. Correlations of the calculated ∆EPart with the log(Ks) (the terminology used in ref 3 for log(Kd) onto soot) values for “added” PAHs reported in ref 3 were also good. Correlations for added PAHs from (3) were also performed for the SRM 1650 diesel soot. The correlation is good and the two lines are roughly parallel. The log(Kd) values are systematically lower than those in ref 2; hence, there is an offset between the two correlation lines. In addition, the predicted value for benzo[a]pyrene log(Kd)sa compound not used in ref 2sfrom the correlation with the data in ref 2 was 8.21. This is in the middle of the range of values measured for traffic, oil, wood, coal, and diesel soots measured in ref 3 which range from 7.41 (coal soot) to 8.87 (traffic soot). The predicted log(Kd) value for anthracene is 6.37 which is in good agreement with the measured value of 6.28 for native anthracene from (3) for SRM 1650 (the substrate used in ref 2). The oil, wood, and coal soot values are significantly lower as well as the added PAH values for soots except for the traffic soot (log(Kd) ) 6.27). Comparison of calculated ∆EPart versus observed log(Ks) for anthracene, phenanthrene, fluoranthene, pyrene, and benzo[a]pyrene with oil, wood, and coal soots (3) resulted in R2 values of 0.94, 0.95, and 0.78, respectively. Two suspect data points lowered the quality of the fit. The coal soot log(Ks) values for both fluoranthene and pyrene are less than that for anthracene (5.73 and 5.71 vs 5.74). Although only three points of comparison were available, the correlation between observed log(Ks) values for added PAHs (anthracene, phenanthrene, and fluorene) onto SRM 1650 diesel soot and calculated ∆EPart values gave an R2 of 0.99. VOL. 40, NO. 7, 2006 / ENVIRONMENTAL SCIENCE & TECHNOLOGY
9
2301
TABLE 3. MP2/6-31G(d) Calculations in COMPASS Force Field Minimum Energy Configurations of Isolated Molecules and in Association with Coronenea,b model
EGas (hartree)c
EAq (hartree)c
Eads (hartree)c
∆ESolv (kJ/mol)
ECav (kJ/mol)
∆EAds (kJ/mol)
∆EPart (kJ/mol)
benzene naphthalene fluorene anthracene phenanthrene pyrene fluoranthene benzo[a]pyrene benzo[a]anthracene benzo[k]fluoranthene benzo[g,h,i]perylene coronene
-231.4572 -384.6136 -499.7716 -537.7637 -537.7732 -613.7824 -613.7161 -766.9432 -690.9264 -766.9070 -842.9153 -918.9777
-231.4632 -384.6221 -499.7812 -537.7747 -537.7842 -613.7936 -613.7279 -766.9515 -690.9398 -766.9212 -842.9280
-1150.4514 -1303.6166 -1418.7781 -1456.7723 -1456.7817 -1532.7946 -1532.7284 -1685.9560 -1609.9408 -1685.9234 -1761.9354
-16 -22 -25 -29 -29 -29 -31 -22 -35 -37 -33
+59 +84 +101 +109 +105 +111 +112 +133 +128 +138 +136
-43 (-30) -66 (-48) -76 (-57) -81 (-59) -81 (-59) -91 (-67) -91 (-72) -92 (-65) -96 (-70) -102 (-75) -111 (-86)
-86 (-73) -128 (-110) -152 (-133) -161 (-139) -157 (-135) -173 (-149) -172 (-153) -197 (-170) -189 (-163) -203 (-176) -214 (-189)
a Solvation performed on MP2/6-31G(d) energy minima of isolated molecules in a dielectric continuum using IEFPCM, ref 16. Values of ∆E Ads in parentheses include BSSE, ref 21. ECav is the cavitation energy in water only. b ∆ESolv ) EAq - EGas; ∆EAds ) EAds - EGas - Coronene(g); ∆EPart ) c ∆EAds - ∆ESolv - ECav. -1 hartree/molecule - 2625.5 kJ/mol.
The fact that the model partitioning energies correlate strongly with the experimental log(Kd) values suggests that the main assumption made in constructing the modelss approximating the soot surface with a coronene molecules was justified. Consequently, one can interpret from the correlation in Figure 3 that the main mechanism of interaction between these PAHs and soot in the experiments is van der Waals forces between the π-electrons within the aromatic rings of both PAH and soot. With van der Waals forces between π-electrons driving the adsorption interaction, the larger, more planar, and more aromatic the adsorbate, the stronger the partitioning should be. This conclusion is consistent with the interpretation of Jonker and Koelmans (3); however, the current study only modeled surface adsorption and did not take into account partitioning that may occur into micropores. Movement of PAHs into micropores could enhance log(Kd) because both sides of the PAH can form π-π interactions with the soot. Although the alignment of a PAH molecule on an aromatic soot surface can explain the stronger affinity of these compounds for soot compared to that of generic natural organic matter (1), it cannot explain the long-term preservation of PAHs in soils and sediments (e.g., 29). Because PAHs adsorbed onto a soot surface by this mechanism would be in equilibrium with pore water, a significant kinetic barrier to their eventual water solubility and biodegradation would not exist (unless physically sequestered with a soil or sediment aggregate). Consequently, one can infer that the long-term presence of PAHs associated with soot in soils and sediments is due to their location inside the soot particle which probably occurs at the time of formation of the soot (9). Jonker and Koelmans (3) demonstrated this effect by comparing native and deuterated PAH log(Kd) values for a variety of soots. The native PAH log(Kd) values were generally lower than the deuterated PAH values which were only added by adsorption. From these data, Jonker and Koelmans (3) concluded that a fraction of native PAHs is occluded within soot. If this is the case, then one can also conclude that the long-term presence of PAHs in soils and sediments is more likely to be associated with combustion processes than with spills of petroleum products. These correlations were used to predict log(Kd) values for benzene, anthracene, fluoranthene, and benzo[a]pyrene using the MP2/6-31G(d) calculated partitioning energies. Measuring distribution coefficients can be difficult in the laboratory, especially at the extreme values. For example, Bucheli and Gustafsson (2) note problems with volatility in determining the log(Kd) for naphthalene. Thus, it is worthwhile to have a computational method for estimating these 2302
9
ENVIRONMENTAL SCIENCE & TECHNOLOGY / VOL. 40, NO. 7, 2006
parameters. On the other end of the scale, compounds such as benzo[a]pyrene may have extremely low water solubilities, so it becomes difficult to accurately obtain their concentrations dissolved in water. The approach discussed here is not subject to these limitations. Classical molecular mechanics simulations can fail to predict thermodynamic properties in two ways. First, if an unrealistic physical model is used, such as the nanodroplets employed in an earlier study (22), then the results will not reproduce experiment. In this study, the MD simulations were inaccurate due to inaccuracies in the force field that describes the potential energy of atomic interactions. Comparison to ab initio calculations can be used to estimate these errors and determined the components of the force field responsible for them (e.g., H-bonding interactions). Although MP2 calculations are time-consuming, computing values for the ∆EPart can be a reliable method for predicting partitioning behavior of PAHs in the environment. Future experimental work could measure these values and test the reliability of the model predictions presented here.
Acknowledgments This research was funded by the NSF Grant “Stony BrookBNL Collaboration to Establish a Center for Environmental Molecular Sciences”. Computation was supported in part by the Materials Simulation Center, a Penn-State MRSEC and MRI facility, and the Center for Environmental Kinetics Analysis, an NSF/DOE Environmental Molecular Sciences Institute. The author thanks the reviewers for comments that improved the content and presentation of the manuscript.
Literature Cited (1) Gustafsson, O.; Haghseta, F.; Chan, C.; MacFarlane, J.; Gschwend, P. M. Quantification of the dilute sedimentary soot phase: Implications for PAH speciation and bioavailability. Environ. Sci. Technol. 1997, 31 (1), 203-209. (2) Bucheli, T. D.; Gustafsson, O. Quantification of the soot-water distribution coefficient of PAHs provides mechanistic basis for enhanced sorption observations. Environ. Sci. Technol. 2000, 34 (24), 5144-5151. (3) Jonker, M. T. O.; Koelmans, A. A. Sorption of polycyclic aromatic hydrocarbons and polychlorinated biphenyls to soot and sootlike materials in the aqueous environment: Mechanistic considerations. Environ. Sci. Technol. 2002, 36 (17), 3725-3734. (4) Hong, L.; Ghosh, U.; Mahajan, T.; Zare, R. N.; Luthy, R. G. PAH sorption mechanism and partitioning behavior in lampblackimpacted soils from former oil-gas plant sites. Environ. Sci. Technol. 2003, 37 (16), 3625-3634. (5) Nguyen, T. H.; Sabbah, I.; Ball, W. P. Sorption nonlinearity for organic contaminants with diesel soot: Method development
(6)
(7) (8)
(9) (10)
(11)
(12) (13)
(14) (15) (16)
(17) (18)
and isotherm interpretation. Environ. Sci. Technol. 2004, 38 (13), 3595-3603. Allen-King, R. M.; Grathwohl, P.; Ball, W. P. New modeling paradigms for the sorption of hydrophobic organic chemicals to heterogeneous carbonaceous matter in soils, sediments, and rocks. Adv. Water Res. 2002, 25 (8-12), 985-1016. Hatzinger, P. B.; Alexander, M. Effect of aging of chemicals in soil on their biodegradability and extractability. Environ. Sci. Technol. 1995, 29 (2), 537-545. Apitz, S. E.; Arias, E.; Clawson, S. A.; Lin, E. W.; Melcher, R. J.; Hemmingsen, B. B. The development of a sterile, PAH-spiked, aged marine sediment for biodegradation experiments: chemical results. Org. Geochem. 1999, 30 (8B), 891-900. Dobbins, R. A.; Fletcher, R. A.; Lu, W. Laser microprobe analysis of soot precursor particles and carbonaceous soot. Combust. Flame 1995, 100 (1-2), 301-309. Reddy, C. M.; Pearson, A.; Xu, L.; McNichol, A. P.; Benner, B. A.; Wise, S. A.; Klouda, G. A.; Currie, L. A.; Eglinton, T. I. Radiocarbon as a tool to apportion the sources of polycyclic aromatic hydrocarbons and black carbon in environmental samples. Environ. Sci. Technol. 2002, 36 (8), 1774-1782. Wade, T. L.; Velinsky, D. J.; Reinharz, E.; Schlekat, C. E. Tidal river sediments in the Washington, D. C. area. II. Distribution and sources of organic contaminants. Estuaries 1994, 17 (2), 321-333. Cerius2, version 4.9; Accelrys Inc.: San Diego, CA, 2003. Sun, H. COMPASS: An ab initio force-field optimized for condensed-phase applicationssOverview with details on alkane and benzene compounds. J. Phys. Chem. B 1998, 102 (38), 73387364. Akhter, M. S.; Chughtai, A. R.; Smith, D. M. The structure of hexane soot I: Spectroscopic studies. Appl. Spectrosc. 1985, 39 (1), 143-153. Akhter, M. S.; Chughtai, A. R.; Smith, D. M. Structure of hexane soot II: Extraction studies. Appl. Spectrosc. 1985, 39 (1), 154167. Cance`s, E.; Mennucci, B.; Tomasi, J. A new integral equation formalism for the polarizable continuum model: Theoretical background and applications to isotropic and anisotropic dielectrics. J. Chem. Phys. 1997, 107 (8), 3032-3041. Frisch, M. J. T.; et al. Gaussian 03, revision C.01 ed.; Gaussian, Inc.: Wallingford, CT, 2004. Keith, T. A.; Frisch, M. J. Inclusion of explicit solvent molecules in a self-consistent-reaction field model of solvation; ACS Symp. Ser. 1994, 569, 22-35.
(19) Møller, C.; Plesset, M. S. Note on an approximation treatment for many-electron systems. Phys. Rev. 1934, 46, 618-622. (20) Ditchfield, R.; Hehre, W. J.; Pople, J. A. Self-consistent molecularorbital methods. IX. An extended Gaussian-type basis for molecular-orbital studies of organic molecules. J. Chem. Phys. 1971, 54 (2), 724-728. (21) Boys, S. F.; Bernardi, F. The calculation of small molecular interactions by the differences of separate total energies. Some procedures with reduced errors. Mol. Phys. 1970, 19 (4), 553566. (22) Kubicki, J. D. Molecular modeling of soot and interactions with polycyclic aromatic hydrocarbons. Abstr. Pap.sAm. Chem. Soc. 2000, 220, ENVR 24, Part 1. (23) Pelmenschikov, A.; Leszczynski, J. Adsorption of 1,3,5-trinitrobenzene on the siloxane sites of clay minerals: Ab initio calculations of molecular models. J. Phys. Chem. B 1999, 103 (33), 6886-6890. (24) Zander, M. Physical and chemical properties of polycyclic aromatic hydrocarbons. In Handbook of Polyaromatic Hydrocarbons; Bjorseth, A., Ed.; Marcel Dekker: New York, 1983; pp 1-25. (25) Pontolillo, J.; Eganhouse, R. P. The search for reliable aqueous solubility (Sw) and octanol-water partition coefficient (K+) data for hydrophobic organic compounds: DDT and DDE as a case study; USGS Water-Resources Investigations Report 01-4201; 2001; pp 1-51. (26) Simon, S.; Duran, M.; Dannenberg, J. J. How does basis set superposition error change the potential surfaces for hydrogen bonded dimers? J. Chem. Phys. 1996, 105 (24), 11024-11031. (27) Ruuska, H.; Pakkanen, T. A. Ab initio study of interlayer interaction of graphite: Benzene-coronene and coronene dimer two-layer models. J. Phys. Chem. B 2001, 105 (39), 9541-9547. (28) Cramer, C. J.; Truhlar, D. G. PM3-SM3: A general parametrization for including aqueous solvation effects in the PM3 molecular orbital model. J. Comput. Chem. 1992, 13 (9), 1089-1097. (29) McGroddy, S. E.; Farrington, J. W. Sediment porewater partitioning of polycyclic aromatic hydrocarbons in three cores from Boston Harbor, Massachusetts. Environ. Sci. Technol. 1995, 29 (6), 1542-1550.
Received for review June 8, 2005. Revised manuscript received February 1, 2006. Accepted February 3, 2006. ES051083S
VOL. 40, NO. 7, 2006 / ENVIRONMENTAL SCIENCE & TECHNOLOGY
9
2303
