Modeling Sorption and Diffusion of Organic Sorbate in

Feb 18, 2013 - Qian Zhao† and Susan E. Burns*. † School of .... Qing Zhou , Runliang Zhu , Stephen C. Parker , Jianxi Zhu , Hongping He , Marco Mo...
0 downloads 0 Views 2MB Size
Article pubs.acs.org/est

Modeling Sorption and Diffusion of Organic Sorbate in Hexadecyltrimethylammonium-Modified Clay Nanopores − A Molecular Dynamics Simulation Study Qian Zhao† and Susan E. Burns* †

School of Civil and Environmental Engineering, 790 Atlantic Drive, N.W., Georgia Institute of Technology, Atlanta, Georgia 30332-0355, United States ABSTRACT: Organoclays are highly sorptive engineered materials that can be used as amendments in barrier systems or geosynthetic liners. The performance of confining and isolating the nonpolar organic contaminants by those barrier/lining systems is essentially controlled by the process of organic contaminant mass transport in nanopores of organoclays. In this article, we use molecular dynamics (MD) simulations to study the sorption and diffusion of organic sorbates in interlayers of sodium montmorillonite and hexadecyltrimethylammonium (HDTMA+)-modified montmorillonite clays. Simulated system consisted of the clay framework, interlayer organic cation, water, and organic sorbate. Their interactions were addressed by the combined force field of ClayFF, constant-valence force field, and SPC water model. Simulation results indicated that in HDTMA coated clay nanopores, diffusion of nonpolar species benzene was slowed because they were subjected to influence of both the pore wall and the HDTMA surfactant. This suggested the nonpolar organic compound diffusion in organophilic clays can be affected by molecular size of diffusive species, clay pore size, and organic surfactant loading. Additionally, a model that connected the diffusion rate of organic compounds in the bulk organoclay matrix with macropores and nanopores was established. The impact of intercalated organic cations on the diffusion dominated mass transport of organic compounds yielded insight into the prediction of the apparent diffusion behavior of organic compounds in organic-modified clays.



INTRODUCTION Organic contaminants in groundwater can originate from a variety of sources such as industrial manufacture, human activities, and improper waste disposal practices.1,2 Because many of these compounds are toxic and/or carcinogenic, the release of organic contaminants to the environment, along with their interaction with soil sediments, transport, and final fate is of great concern.3,4 Organic-coated clays, or organoclays, can be used as amendments in engineered barriers or geosynthetic clay liners for retardation, sorption, and isolation of organic compounds due to their strong capacity for uptake of nonpolar organic compounds.5−8 Most commonly, organic-coated clays are created by exchanging organic cations onto the clay surface, through replacement of the naturally occurring inorganic metallic cations of the base clays. When organic cations replace inorganic cations on the clay surface, organoclays usually exhibit increased saturated hydraulic conductivity compared to that of unmodified clays; however, their conductivity for water is still low (magnitude of 10−8 m/s) and may even be reduced for separate phase organic fluids.6 Consequently, the rate of advective transport remains low, and the molecular diffusion caused by concentration gradients in porous media like an organoclay can be the critical component in the overall rate of mass transport. In many instances, diffusion controls the quantity and rate of organic compound that passes through a © 2013 American Chemical Society

barrier system, especially in conditions of low hydraulic gradient.9 The organoclay capacity to uptake nonpolar solutes allows retardation of the diffusive transport of the solute and leads to a diffusion regime that can be considered as a transient process.10,11 Typically, nonpolar compound diffusion through clays with sorption is described using Fick’s law coupled with a sink term that accounts for sorption of the nonpolar solute:12,13 ε

∂q ∂C ∂ 2C = De 2 − ρ ∂t ∂t ∂x

(1)

where ρ((∂q)/(∂t)) is the sink term that corresponds to transient diffusion, ε and ρ denote the porosity and bulk density of the organoclay, respectively, and De is the effective diffusion coefficient, including the overall effect of tortuous and restricted movement of solute in the porous media:

De =

Daq εeδ

Received: Revised: Accepted: Published: 2769

(2)

τ November 8, 2012 February 12, 2013 February 18, 2013 February 18, 2013

dx.doi.org/10.1021/es3045482 | Environ. Sci. Technol. 2013, 47, 2769−2776

Environmental Science & Technology

Article

Figure 1. Interlayer diffusion process controlled by interaction between organic solute and interlayer organic surfactant.

where εe is effective porosity, δ and τ are the constrictivity factor and tortuosity factor of the porous medium, respectively, and Daq is the diffusion coefficient in bulk water. For organophilic partitioning clays,14,15 sorption of hydrophobic organic compounds in interlayer organic surfactants can be quantified by a simple distribution/partitioning coefficient Kd:

q = KdCaq

Da, i Daq, i

1 (αmacropore, i + αnanopore, iδnanopore, i) Gi

(5)

where G is a geometric factor that accounts for pore characteristics and distribution, δ is a constrictivity factor that accounts for the constricted movement of species in interlayer, and α is the molar fraction of species of interest in different types of pores. Because the microstructure of the porous media (e.g., effective pore volume, tortuosity) remains constant, Bourg and Sposito (2010) suggested the geometric factor could be canceled out by normalizing cation (i.e., Na+, Sr+, Cs+) diffusion with water diffusion. The advantage of this approach is that when the molar fraction and diffusive properties of each species in interlayer nanopores is known, the diffusive properties of them in bulk can be predicted by normalizing the water tracer. Similarly, the diffusion of nonpolar compounds in organoclays can be predicted if the sorption in interlayer pores is taken into account in the equation. For organophilic partitioning clays, organic solutes are expected to have lower average diffusion rate in interlayer pores due to sorption by the interlayer organic carbon phase (e.g., HDTMA cations). Consequently, interlayer diffusion becomes the rate limiting step for overall mass transport due to the higher concentration of organic solute in organoclay interlayer assuming no preferential transport paths for organic solute are available (Figure 1). It is the goal of this work to investigate the sorption and transport of benzene in interlayer pores of sodium montmorillonite (Na−MMT) and HDTMA ((CH3)3N(CH2)15CH3)modified organoclay (HDclay) to relate the apparent diffusion coefficient of interlayer species Dinterlayer to the bulk diffusion coefficient D in organoclay. Eq 5 is modified to take into account sorption of benzene, which is discriminated from noninteractive species like water or sodium:

(3)

where q is the mass of solute per unit dry mass of solid and Caq is the concentration of solute in solution at equilibrium [mol/ L]. Combing the transport equation, with the sorption equation and the apparent diffusion coefficient, Da, yields: ⎛ De ⎞ ∂ 2C ∂C ∂ 2C =⎜ ⎟ 2 = Da 2 ∂t ∂x ⎝ ε + Kdρ ⎠ ∂x

=

(4)

The diffusion equation for a nonpolar sorbate diffusing in an organoclay requires many parameters, some of which are difficult to obtain from experimental work. Additionally, many isotherms exhibit nonlinearity as a function of concentration, which is not addressed in this model. Most importantly, the model was derived based on the bulk properties of the porous media; consequently, the interlayer pore characteristics and the equilibrium between macropores and interlayer pores are not reflected in this equation. As a result, the modeled parameters deviate from values obtained from experiments; for example, Kd values observed in diffusion experiments are inherently lower than the value of Kd obtained from batch sorption tests.16 The overestimated K d can be problematic when designing retardation capacity and equilibrium time of transient diffusion for barrier systems. Additionally, previous studies have shown that in organophilic pores (e.g., montmorillonite with organic coatings) surface diffusion of adsorbed species can also occur in interlayer pores. This requires an extra surface diffusion term be included in Da to fit some experimental data, which causes additional difficulties in the application of this model.17−19 Recently, a diffusion model that focuses on the interlayer diffusive behavior of species and connects interlayer “compartment” pore with continuum scale has been proposed by Bourg and Sposito (2010).20 The new model described the diffusion coefficient of object species (i) as the weighted average of that in the macropore and in the interlayer nanopore:

Da,benzene Daq,benzene

=

1⎛ ⎜αmacropore,benzene G⎝ +

αnanopore,benzeneδnanopore,benzene ⎞ ⎟ R ⎠

(6)

where R is the term that accounts for the confined movement of benzene in HDclay due to sorption. On the basis of this model, molecular dynamics simulations were performed to study the interlayer diffusive behavior of water tracer and benzene in pure sodium-montmorillonite and 2770

dx.doi.org/10.1021/es3045482 | Environ. Sci. Technol. 2013, 47, 2769−2776

Environmental Science & Technology

Article

Table 1. HDclay Properties and Interlayer Diffusion Coefficientsa clay type HDclay1 HDclay2 HDclay3 Na−MMT a

cation +

HDTMA HDTMA+ HDTMA+ HDTMA+ Na+ Na+

diffusing species

water molecules per unit cell

basal spacing (Å)

interlayer diffusion coefficient (m2/s)

benzene benzene tracer benzene benzene tracer

∼5 ∼10

22.468 ± 0.045 25.211 ± 0.041

∼15 ∼22.5

28.023 ± 0.054 25.595 ± 0.050

1.40 ± 0.33 × 10−10 1.78 ± 0.37 × 10−10 1.18 ± 0.27 × 10−9 2.15 ± 0.44 × 10−10 4.11 ± 0.72 × 10−10 1.30 ± 0.19 × 10−9

Chemical formula of clay unit cell: (Cation)0.643[Al0.071Si7.929O16][Mg0.572Al3.429O4(OH)4].

periodic boundary in all three directions. A single interlayer space was modeled during the transport simulations only, and was chosen to reduce simulation times.

HDclays with different amounts of interlayer water. An advanced force field of clay was used, and each simulation was as long as 5 ns. Results of all simulations were compiled to examine the: (1) constrictivity on water tracer and organic sorbate (benzene) in the interlayer nanopores of HDclay, (2) sorption and diffusion dynamics of organic sorbate in interlayer by addressing the interactions between different phases and apparent diffusion coefficients of organic molecules, (3) relation of transport rate of organic compound in organoclay interlayer to apparent values in the bulk. On the basis of the results, the apparent bulk diffusion coefficient of organic sorbate in organophilic clays can be predicted when the diffusion coefficient of a water tracer is known.



MOLECULAR DYNAMICS SIMULATIONS Force Field. Force-field based, equilibrium molecular dynamics simulations of the HDclay interlayer were performed in the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS).21 The basis of the overall simulation was the combined force field of ClayFF22 + Constant Valence Force Field (CVFF23) + SPC water model (included in ClayFF).22,24,25 The combined force field of the three has been proven appropriate for the simulation of systems that contain clay, water, and organic phases,26−28 and for modeling the interaction between phases in the interlayer. This is especially important because the diffusion behavior of the species are addressed based on the combined force field, with the accuracy of the process having been verified through previous studies.20 The energy terms of the clay and HDTMA cations have been described in detail in previous work by the authors.29 Physical Model. The clay prototype was modified from MMT321 (southern clay30), and had a chemical formula of (Na)0.643[Al0.071Si7.929O16][Mg0.572Al3.429O4(OH)4], with a clay platelet that consisted of 28 unit cells. The initial size of the simulation cell was approximately 35.72 Å × 35.04 Å in the basal plane and approximately 22.5−28 Å in the z direction depending on the quantities of interlayer species. The net negative charge of clay platelet resulted from 18 random isomorphic substitutions, and the negative charge of the clay mineral framework was compensated by the interlayer organic cations (assuming 100% cation exchange) to maintain the electrical neutrality of the simulated system. The interlayer space consisted of species including water, HDTMA cations, and benzene. A Na−MMT system with 22.5 interlayer water molecules per unit cell (approximately 630) was also simulated for comparison with the organoclays. Three sets of simulations were performed for HDclays with interlayer water contents of 5 H2O, 10 H2O, and 15 H2O per unit cell (Table 1), which was approximately equivalent to the amount of water in one, two, and three layer hydrates in Na−MMT.20 For the transport study, the simulation box included only one interlayer space sandwiched between two half clay platelets (Figure 2) with

Figure 2. 3D view of stage 2 simulation box (periodic boundaries highlighted by black lines).

Simulated System. Each simulated system included approximately 2500−4500 atoms and the overall simulation time was more than 6 ns. The simulation time was sufficiently long to capture the diffusive behavior of interlayer species after equilibrium,20 as also verified by similar results with longer simulation time.33 In the first stage of simulation, a 2-clay 2interlayer structure was subjected to NPT and NVT time integration and had free gallery spacing equilibrium following the method described in Zhao and Burns (2012).34 In addition to benzene molecules and HDTMA cations, water molecules were also intercalated into each interlayer space. After stage 1, the simulated system was well-optimized because the interlayer species had relatively more freedom to move, and they can be used as initial configurations for the diffusion simulations.33,36 After the simulated system reached equilibrium after stage 1, a new simulation cell was created by subtracting the previous 2interlayer and 2-clay platelets to 1 interlayer and 1 clay platelet (1/2 on top and bottom). This subtracted equilibrium scenario of the sandwiched layer was then imported to LAMMPS and optimized again followed by another 500 ps of NPT time integration and 5 ns NVT time integration. During all stages, the time step was set at 1.0 fs, and for stage 2 NVT time integration trajectories were recorded every approximately 100 ps. The mean-squared displacement (MSD) of interlayer species were recorded every 1.0 ps for postanalysis. 2771

dx.doi.org/10.1021/es3045482 | Environ. Sci. Technol. 2013, 47, 2769−2776

Environmental Science & Technology

Article

Figure 3. Interlayer species of simulated system after equilibrium. (a) Pure montmorillonite, (b) interlayer water equivalent to 5 per unit cell, (c) interlayer water equivalent to 10 per unit cell, (d) interlayer water equivalent to 15 per unit cell. Note: for all species, C = gray, N = blue, H = white, O = red, Si = yellow, Na = purple, Al = pink, Mg = green. Benzene molecules are highlighted, and blue dashed line = hydrogen bond.



RESULTS AND DISCUSSION Simulated Systems at Equilibrium. Trajectories of HDTMA-organoclays and their nanopore species were obtained at the end of NVT time integration, yielding visualization of the interlayer species (Figure 3) and basal spacing at equilibrium (Table 1). It was noted that the interlayer spacing increased from HDclay1 to HDclay3 (parts b−d of Figure 3) indicating increased size of interlayer pore volume due to the incorporation of additional water molecules. Additionally, the nitrogen head groups of HDTMA cations were attracted to the clay surface, whereas the aliphatic carbon chains were driven away from the clay surface as the number of interlayer water molecules increased. The hydrogen bonds (dashed blue line, Figure 3) formed between the water and HDclay surface oxygen was increased significantly as the amount of interlayer water increased, and the first layer of water molecules bonded to clay surface had lower mobility than the other water molecules, which displaced the clay-aliphatic chain interaction. Also, it can be noticed that as the amount of interlayer water increased, water molecules developed H bonds that connected adjacent water molecules, and this H bond network phenomenon prevailed as more water molecules were introduced into the interlayer. Finally, at higher levels of interlayer water, the interlayer volume partitioned into small domains of organic phase (HD cation) versus water phase (part d of Figure 3). Within the interlayer space, the predominant mode of interaction between organic surfactant coated clay and benzene molecules was between the HD organic cation chains and benzene, as opposed to interaction between the clay surface and benzene, as would be anticipated for a partitioning clay. In contrast, interaction between the clay surface and the benzene molecules would be observed in the case of adsorptive clays29,33

(e.g., TMA clay). Because the aliphatic chain of the organic cation had relatively larger mobility than the nitrogen headgroup, the benzene molecules had mobility as the chain moved, even in the case where they were sorbed via short-range forces. Interlayer Diffusion Coefficients Calculated from MSD. The mean-squared displacement (MSD) of benzene molecules during NVT time integration in all the simulated systems was recorded, and the benzene interlayer diffusion coefficient was calculated based on the Einstein relation: Dinterlayer =

dMSD 4dt

(7)

For each simulated system, the apparent interlayer diffusion coefficients of the benzene and tracer species were obtained between 5−5000 ps of simulation time, and the mean value was calculated as the slope of MSD versus time at: 5, 6, 7,...10 ps; 20, 30,.... 100 ps; 200, 300,.... 1000 ps; 2000, 3000,...5000 ps. The diffusion coefficients plotted versus time were shown in Figure 4 (standard deviation too small to be shown on log−log scale). The calculated diffusion coefficient for all species clearly showed a hindered diffusion stage up to 100 ps of simulation time, as demonstrated in previous studies;20 consequently, the calculated diffusion coefficients after 100 ps were taken as the representative values for subsequent analysis. The calculated diffusion coefficients indicated that the diffusion rate in interlayer nanopores was a function of both solute and interlayer pore characteristics. Note that the water tracer diffused much faster than benzene in both Na−MMT and HDclays but slower than in bulk water. In all interlayer cases, benzene also diffused much more slowly than when compared to benzene in bulk aqueous solution. Other than interaction 2772

dx.doi.org/10.1021/es3045482 | Environ. Sci. Technol. 2013, 47, 2769−2776

Environmental Science & Technology

Article

rate for all benzene molecules. At the macro-scale, the increased organic surfactant loading into the clay interlayer often results in a higher uptake capacity of organic compounds,14 and a larger retardation factor for partitioning clays. In this study, for the three simulated HDclays, the simulations with lower quantity of water also demonstrated higher packing density of aliphatic carbons in the interlayer, which resulted in even slower motion of benzene molecules. Because of preferential interaction with aliphatic carbon chains of HD cations, benzene diffusion in the HDclay interlayer is not only subjected to the influence of the clay surface but also to the HD cations, with the latter becoming even more significant at intermediate to large amounts of cation loading. In addition to the possible interaction with the organic compound, the existence of HDTMA cations in MMT interlayer may have other impacts on their diffusive movement. The increased packing density of HDTMA cations in the HDclay interlayer also increased the tortuosity of the pore space and decreased the mean width of available free-path, which was reflected in a decreased rate of diffusion of benzene molecules from HDclay3 to HDclay1 (Figure 4). When the amount of interlayer water increased, HDTMA cations behaved more like pillars than compartment walls (Figure 5). For constant loading of HDTMA cations, because such influence on effective diffusion path was insignificant compared with sorption, its impact on the diffusivity was combined with influence of the pore wall and included in the constrictivity term. Consequently, HDclays with same amount of organic loading were assumed to have the same sorption factor R, which accounted for the fraction of interlayer sorbates that became relatively immobile. To separate the R factor from other factors, HDclay2 and Na−MMT were compared on the basis of their similar pore sizes. The normalized benzene diffusion coefficient for the Na−MMT interlayer to the HDclay2 interlayer yielded (Figure 6):

Figure 4. MD estimated diffusion coefficient of water tracer in Na− MMT; HDclay2 and benzene in: Na−MMT; HDclay1; HDclay2; HDclay3.

between species, the slow motion of water tracer or benzene in interlayer nanopores was caused by constraints imposed by the interlayer pores with limited scale length. For noninteracting species, when Knudsen diffusion dominated in interlayer pores, their interlayer diffusion coefficients could be described directly by a single constrictivity factor multiplied by the diffusion coefficient in aqueous solution, as suggested by Bourg and Sposito (2010):20 Dinterlayer = Daq δ

(8)

As interlayer pore size increased, solute molecules displayed an increased apparent diffusion rate corresponding to increased volume for motion. This was observed for benzene diffusion in HDclays: whereas the HD loading on the clay surface was constant, and interlayer spacing increased due to the addition of water molecules, the diffusion rate of benzene also increased. The water tracer demonstrated a higher diffusion rate than benzene due to its relatively small molecule size. Additionally, benzene molecules were subjected to interaction with the intercalated HD cations, which did not occur for the water tracer. Comparison of two interlayers with similar pore sizes (Na−MMT and HDclay2) demonstrated that benzene diffused faster in Na−MMT than in the HDclay. However, for benzene diffusion in HDclays, even when the size effect of water and benzene molecules was taken into account in terms of constrictivity, the diffusion of benzene molecules was still inherently lower than that calculated from eq 8. It is important to note that eq 8 is only valid for noninteracting species. For benzene diffusion in HDclay pores, the sorption of a fraction of benzene to the HD phases caused a reduction in the average movement of benzene, as would be anticipated from the mass transport equation. Before sorption-diffusion equilibrium was reached within macro and interlayer pores, the diffusion pattern has been proved to be transient. This effect will be discussed in the following section. Confined Movement of Benzene Caused by Sorption and Diffusion Dynamics. Comparison of the simulations of benzene molecules diffusing in pure montmorillonite with benzene molecules diffusing in HD clay clearly showed significant differences. At the nanoscale, the intercalation of the HDTMA cation increased the hydrophobic carbon content in the interlayer; consequently, a fraction of benzene molecules had the tendency to be sorbed onto the aliphatic chains via short-range forces, which resulted in a lower average diffusion

Rt =

D Na − MMT,t DHDclay2,t

(9)

Note that R is approximately 2−2.5, which suggested that a fraction of benzene molecules had very low mobility and the averaged mobility of all benzene molecules was decreased, when compared with the benzene molecules in Na−MMT. The remaining fraction of benzene molecules that desorbed or dissolved in interlayer water had faster rates of diffusion. Constrictivity Imposed by Interlayer Pores. Interlayer water and cations exhibit decreased apparent diffusion coefficients in nanopores filled with liquid compared with diffusion in bulk aqueous solution. This restricted movement of sorbate results from increased viscosity of the solvent and would be more pronounced when the pore size is small compared to the mean free path of the sorbate. Although the sorbed HDTMA cations also occupied considerable interlayer pore volume, and they imposed a drag force on benzene, it was clear that they might contribute to the restricted movement of benzene. However, in this study, the effects will be considered as sorption, which was included in the R term. Consequently, for the water tracer and benzene diffusion in clay interlayer nanopores, it was assumed that the constrictivity only corresponded to the proximity of the 2:1 clay surface. The constricitivity factor δ will be herein referred to the constraint imposed by the mineral phase only: 2773

dx.doi.org/10.1021/es3045482 | Environ. Sci. Technol. 2013, 47, 2769−2776

Environmental Science & Technology

Article

Figure 6. Sorption factor R plotted versus simulation time.

and Schultz (1970) (eq 11) (1973) 32 (eq 12):

31

and Satterfield and Colton

δ = (1 − λp)4

(11)

δ = exp( −4.6λp)

(12)

Where: λp = ((molecule diameter)/(pore diameter)) < 1. In this study, the movements of water tracer and benzene molecules in the interlayer nanopores were recorded, and the constrictivitiy of water tracer was calculated as the relative diffusion rate in the clay interlayer to that in aqueous phase, by normalizing Dinterlayer,water (t > 100 ps) with Daq,water. The constrictivitiy for benzene was calculated in the same way, except in HDclays, interlayer diffusion rate was multiplied by R term to eliminate the effect of sorption, for example RavgDinterlayer,benzene (t > 100 ps) normalized with Daq,benzene (Figure 7).

Figure 5. Top view of interlayer species density cloud (stage 2, 1−1.05 ns): blue, aliphatic carbon; black, benzene carbon.

Dinterlayer =

Figure 7. Constrictivity factor for solutes vs ratio of solutes diameter to interlayer pore diameter.

Daq δ R

(10)

Where R is 1 for noninteracting polar species (species have no significant short-range interaction with HD cations, e.g., water tracer). Empirical estimates of the constrictivity factor for both steady state and transient state diffusion have been proposed by Beck

Because HDclay had larger interlayer spacing as the number of interlayer water molecules increased, there was less constriction imposed on the diffusion of the interlayer species including both water and benzene. HDclay, with very small interlayer spacing (e.g., HDclay1), had very limited effective interlayer pore space for intermolecular collision and 2774

dx.doi.org/10.1021/es3045482 | Environ. Sci. Technol. 2013, 47, 2769−2776

Environmental Science & Technology

Article

measured effective diffusion coefficients for water tracer diffusion in GCLs with varying organobentonite amendment as observed by Lorenzetti et al. (2005).9 Additionally, the bulk diffusion rate of benzene in Na−MMT estimated from this study was 10−10 m2/s, which agreed with previous observations of benzene diffusion through geosynthetic clay liners.35 The trend of slowed diffusion of interacting species in organophilic clays was also consistent with organic sorbate diffusion in organic-rich clays.16 In summary, this study simulated benzene sorption and diffusion in the interlayer of Na−MMT and HDclays that were intercalated with different quantities of interlayer water. It was shown that the factors that controlled the diffusion of organic compounds in the interlayer of the clays included the constrictivity factor, which accounted for the Knudsen diffusion, and the sorption factor that accounted for the immobile fraction of sorbate that was taken up by the intercalated surfactant. A modified model that averaged the diffusion in the clay macropore and the interlayer pore was proposed, and it can be used for the estimate of bulk diffusion behavior of different species.

subsequently demonstrated the strongest constriction on the mobility of benzene molecules. This demonstrated that increasing the water molecules in the interlayer of organoclays increased the rate of diffusion; however, the diffusion rate of benzene through organoclay would be slowed by decreasing the interlayer pore size of the synthesized clays or by increasing the size of the diffusing molecular species. Bulk Diffusion Coefficient Estimated from MD Simulations. The bulk diffusion coefficient of each species of interest in HDclay was calculated based on eq 6. It is important to note that the equation is only applicable to organoclays with concentrated solids, as opposed to slurries or low solids content suspensions. The molar fraction in the macropore (αmacropore) was used as the only variable, and for water tracers it was calculated as the percentage of water molecules in the macropore. The bulk diffusion coefficient in Figure 8 was calculated based on an assumed geometric factor



AUTHOR INFORMATION

Corresponding Author

*Phone: (404) 894-2285, fax: (404) 894-2281, e-mail: sburns@ gatech.edu. Notes

The authors declare no competing financial interest.



REFERENCES

(1) Bedding, N. D.; McIntyre, A. E.; Perry, R.; Lester, J. N. Organic contaminants in the aquatic environment. 1. Source and occurrence. Sci. Total Environ. 1982, 25 (2), 143−167. (2) Focazio, M. J.; Kolpin, D. W.; Barnes, K. K.; Furlong, E. T.; Meyer, M. T.; Zaugg, S. D.; Barber, L. B.; Thurman, M. E. A national reconnaissance for pharmaceuticals and other organic wastewater contaminants in the United States - II) Untreated drinking water sources. Sci. Total Environ. 2008, 402 (2−3), 201−216. (3) Weber, W. J.; McGinley, P. M.; Katz, L. E. Sorption phenomena in subsurface systems - concepts, models and effects on contaminant fate and transport. Water Res. 1991, 25 (5), 499−528. (4) Pignatello, J. J.; Xing, B. S. Mechanisms of slow sorption of organic chemicals to natural particles. Environ. Sci. Technol. 1996, 30 (1), 1−11. (5) Smith, J. A.; Jaffe, P. R. Benzene transport through landfill liners containing organophilic bentonite. J. Environ. Eng.-ASCE 1994, 120 (6), 1559−1577. (6) Lo, I. M. C.; Yang, X. Y. Use of organoclay as secondary containment for gasoline storage tanks. J. Environ. Eng.-ASCE 2001, 127 (2), 154−161. (7) Yang, L. Y.; Zhou, Z.; Xiao, L.; Wang, X. R. Chemical and biological regeneration of HDTMA-modified montmorillonite after sorption with phenol. Environ. Sci. Technol. 2003, 37 (21), 5057−5061. (8) Lee, S.; Oren, A. H.; Benson, C. H.; Dovantzis, K. Organoclays as variably permeable reactive barrier media to manage NAPLs in ground water. J. Geotech. Geoenviron. Eng. 2012, 138 (2), 115−127. (9) Lorenzetti, R. J.; Bartelt-Hunt, S. L.; Burns, S. E.; Smith, J. A. Hydraulic conductivities and effective diffusion coefficients of geosynthetic clay liners with organobentonite amendments. Geotext. Geomembr. 2005, 23 (5), 385−400. (10) Fetter, C. W. Contaminant Hydrogeology; Prentice Hall: Upper Saddle River, NJ, 1999. (11) Sharma, H. D. Geoenvironmental Engineering: Site Remediation, Waste Containment, and Emerging Waste Management Technologies; Wiley: Hoboken, N.J., 2004.

Figure 8. Estimates of bulk diffusion coefficients based on simulation.

(G) equal to 4, which was equivalent to a bulk totuosity factor of 0.25. G values range from 2.4−5.6, and corresponded to totuosity factors of 0.18−0.42, typical for geomaterials. It is important to note that benzene molecules had much higher concentration in the sorbent than in aqueous phase; consequently, their concentration in macropores was much smaller than that in interlayer pores. However, because the water molecules in the interlayer pores was insignificant compared that in the macorpores, αmacropore,water was close to 1. Overall, αmacropore,benzene≪ αmacropore,water. This suggested that for transient diffusion like benzene through HDclays, the bulk diffusion rate was controlled by interlayer diffusion (αmacropore,benzene ∼ 0). The benchmark to examine the simulation results was the water tracer diffusion rate, which was compared with water tracer diffusion in an HDclay amended GCL.9 The effective diffusion coefficient for water tracer in HDclay amended GCLs as determined by diffusion testing ranged from 3.6−4.0 × 10−10 m2/s, which agrees well with the bulk value extrapolated from this MD simulation result. Furthermore, because the molar fraction of tracer could be dominant in the macropores of both MMT and HD modified MMT, the calculated bulk apparent diffusion coefficient for water tracers should demonstrate no significant difference, as it was controlled by the transport in macropores. This explained the fairly consistent values of 2775

dx.doi.org/10.1021/es3045482 | Environ. Sci. Technol. 2013, 47, 2769−2776

Environmental Science & Technology

Article

(34) Zhao, Q.; Burns, S. E. Microstructure of single chain quaternary ammonium cations intercalated into montmorillonite: A molecular dynamics study. Langmuir 2012, 28 (47), 16393−16400. (35) Lake, C. B.; Rowe, R. K. Volatile organic compound diffusion and sorption coefficients for a needle-punched GCL. Geosynth. Int. 2004, 11 (4), 257−272. (36) Heinz, H.; Paul, W.; Suter, U. W.; Binder, K. Analysis of the phase transitions in alkyl-mica by density and pressure profiles. J. Chem. Phys. 2004, 120 (8), 3847−3854.

(12) Boving, T. B.; Grathwohl, P. Tracer diffusion coefficients in sedimentary rocks: correlation to porosity and hydraulic conductivity. J. Contam. Hydrol. 2001, 53 (1−2), 85−100. (13) Bourg, I. C.; Bourg, A. C. M.; Sposito, G. Modeling diffusion and adsorption in compacted bentonite: A critical review. J. Contam. Hydrol. 2003, 61 (1−4), 293−302. (14) Redding, A. Z.; Burns, S. E.; Upson, R. T.; Anderson, E. F. Organoclay sorption of benzene as a function of total organic carbon content. J. Colloid Interface Sci. 2002, 250 (1), 261−264. (15) Yariv, S.; Cross, H. Organo-Clay Complexes and Interactions; Marcel Dekker: New York, 2002. (16) Grathwohl, P. Diffusion in Natural Porous Media: Contaminant Transport, Sorption/Desorption and Dissolution Kinetics; Kluwer Academic Publishers: Boston, 1998. (17) Thomas, M. M.; Clouse, J. A. Primary migration by diffusion through Kerogen. 1. Model experiments with organic-coated rocks. Geochim. Cosmochim. Acta 1990, 54 (10), 2775−2779. (18) Thomas, M. M.; Clouse, J. A. Primary migration by diffusion through Kerogen. 2. Hydrocarbon diffusivities in Kerogen. Geochim. Cosmochim. Acta 1990, 54 (10), 2781−2792. (19) Thomas, M. M.; Clouse, J. A. Primary migration by diffusion through Kerogen. 3. Calculation of geologic fluxes. Geochim. Cosmochim. Acta 1990, 54 (10), 2793−2797. (20) Bourg, I. C.; Sposito, G. Connecting the molecular scale to the continuum scale for diffusion processes in smectite-rich porous media. Environ. Sci. Technol. 2010, 44 (6), 2085−2091. (21) Plimpton, S. Fast parallel algorithms for short-range molecular dynamics. J. Comput. Phys. 1995, 117 (1), 1−19. (22) Cygan, R. T.; Liang, J. J.; Kalinichev, A. G. Molecular models of hydroxide, oxyhydroxide, and clay phases and the development of a general force field. J. Phys. Chem. B 2004, 108 (4), 1255−1266. (23) Dauberosguthorpe, P.; Roberts, V. A.; Osguthorpe, D. J.; Wolff, J.; Genest, M.; Hagler, A. T. Structure and energetics of ligand-binding to proteins-escherichia-coli dihydrofolate reductase trimethoprim, a drug-receptor system. Proteins: Struct., Funct., Genet. 1988, 4 (1), 31− 47. (24) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; Hermans, J. Interaction Models for Water in Relation to Protein Hydration; Riedel, Dordrecht: The Netherlands, 1981. (25) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Clarendon Press: Oxford, 1987. (26) Greathouse, J. A.; Cygan, R. T. Molecular dynamics simulation of uranyl(VI) adsorption equilibria onto an external montmorillonite surface. Phys. Chem. Chem. Phys. 2005, 7 (20), 3580−3586. (27) Liu, X.; Lu, X. C.; Wang, R. C.; Zhou, H. Q.; Xu, S. J. Molecular dynamics insight into the cointercalation of hexadecyltrimethylammonium and acetate ions into smectites. Am. Mineral. 2009, 94 (1), 143−150. (28) Zhu, R. L.; Chen, W. X.; Shapley, T. V.; Molinari, M.; Ge, F.; Parker, S. C. Sorptive characteristics of organomontmorillonite toward organic compounds: A combined LFERs and molecular dynamics simulation study. Environ. Sci. Technol. 2011, 45 (15), 6504−6510. (29) Zhao, Q.; Burns, S. E. Molecular dynamics simulation of secondary sorption behavior of montmorillonite modified by single chain quaternary ammonium cations. Environ. Sci. Technol. 2012, 46 (7), 3999−4007. (30) Heinz, H.; Koerner, H.; Anderson, K. L.; Vaia, R. A.; Farmer, B. L. Force field for mica-type silicates and dynamics of octadecylammonium chains grafted to montmorillonite. Chem. Mater. 2005, 17 (23), 5658−5669. (31) Beck, R. E.; Schultz, J. S. Hindered diffusion in microporous membranes with know pore geometry. Science 1970, 170 (3964), 1302−&. (32) Satterfi., Cn; Colton, C. K.; Pitcher, W. H. Restricted diffusion in liquids within fine pores. AIChE J. 1973, 19 (3), 628−635. (33) Zhu, R.; Hu, W.; You, Z.; Ge, F.; Tian, K. Molecular dynamics simulation of TCDD adsorption on organo-montmorillonite. J. Colloid Interface Sci. 2012, 377 (1), 328−333. 2776

dx.doi.org/10.1021/es3045482 | Environ. Sci. Technol. 2013, 47, 2769−2776