A New Hypothesis for the Dissolution Mechanism of Silicates

Article pubs.acs.org/JPCC

A New Hypothesis for the Dissolution Mechanism of Silicates James D. Kubicki,*,† Jorge O. Sofo,‡ Adam A. Skelton,§,⊥ and Andrei V. Bandura# †

Department of Geosciences and the Earth & Environmental Systems Institute, The Pennsylvania State University, University Park, Pennsylvania 16802, United States ‡ Department of Physics, The Pennsylvania State University, University Park, Pennsylvania 16802, United States § Department of Chemical Engineering, Vanderbilt University, Nashville, Tennessee 37240, United States ⊥ Department of Chemical and Materials Engineering, University of Dayton, Dayton, Ohio 45469, United States # Department of Chemistry, Quantum Chemistry Division, St. Petersburg State University, St. Petersburg, Russia S Supporting Information *

ABSTRACT: A novel mechanism for protonating bridging O atoms (Obr) and dissolving silica is proposed that is consistent with experimental data and quantum mechanical simulations of the α-quartz (101)/water interface. The new hypothesis is that H+-transfer occurs through internal surface H-bonds (i.e., SiOH−Obr) rather than surface water H-bonds and that increasing ionic strength, I, favors formation of these internal H-bonds, leading to a larger pre-exponential factor, A, in the Arrhenius equation, k = A exp(−ΔEa/RT), and higher rates of dissolution. Projector-augmented planewave density functional theory (DFT) molecular dynamics (MD) simulations and static energy minimizations were performed on the α-quartz (101) surface and with pure water, with Cl−, Na+, and Mg2+. Classical molecular dynamics were performed on α-quartz (101) surface and pure water only. The nature of the H-bonding of the surface silanol (SiOH) groups with the solution and with other surface atoms is examined as a test of the above hypothesis. Statistically significant increases in the percentages of internal SiOH−Obr H-bonds, as well as the possibility of Obr protonation with H-bond linkage to silanol group, are predicted by these simulations, which is consistent with the new hypothesis. This new hypothesis is discussed in relation to experimental data on silicate dissolution.

■

INTRODUCTION α-Quartz (SiO2) is one of the most common and stable minerals on the surface of the Earth, so its dissolution during chemical weathering reactions is of basic importance in geochemistry. Furthermore, quartz is an end-member composition for dissolution or corrosion of silicate minerals and glasses, making it a common material to use as a model system for understanding dissolution of natural and man-made silicates.1,2 For example, some strategies plan for long-term disposal of high-level nuclear waste in silicate glasses3 which requires knowing the mechanisms associated with silicate dissolution to predict kinetics a million years into the future. According to Fruger et al.,4 mechanistic information on silicate dissolution is necessary to confirm the use of simplified models for predicting long-term silicate glass dissolution. This comparison of mineral and glass dissolution is more appropriate than might be assumed based on the differences in mineral and glass structures, but silicate mineral and glass dissolution rates and mechanism can be equal within experimental uncertainty when the glasses are anhydrous.5,6 Consequently, detailed atomistic simulations of a quartz surface with known structure can provide insights into the molecularlevel mechanisms that comprise the silicate dissolution process in general. Rimstidt and Barnes7 studied crystalline and amorphous SiO2 dissolution over a range of hydrothermal temperatures and © 2012 American Chemical Society

began the discussion of silicate dissolution mechanisms in the geochemical community. Since that time, a large number of papers have been published on silica dissolution, but the work of Dove and co-workers is particularly relevant for the discussion in this study. By tracking dissolution rates versus the ΔG of reaction (i.e., undersaturation) while probing changes in surface topography with atomic force microscopy (AFM), Dove et al.8 developed a theory for quartz dissolution that depends upon the step retreat velocity on a crystal surface. Similar nanoscale mechanisms for dissolution have been observed and modeled for Na- and Ca-aluminosilicate minerals by Luttge and co-workers.9,10 The role of ionic strength and surface topography are linked.6 Quantum mechanical modeling of silica hydrolysis and dissolution has been studied for over 20 years.11−19 The silica− water interface has also been the subject of recent classical20−23 and DFT simulations.24−32 Adeagbo et al.19 performed periodic DFT simulations of the dissolution reaction on a quartz surface. The authors studied the α-quartz (100) surface which is significantly different from the (101) surface in that the former is composed of Q2 Si− (OH)2 (i.e., a Si atom with two Si−Obr bonds connecting it to Received: January 18, 2012 Revised: July 16, 2012 Published: July 19, 2012 17479

dx.doi.org/10.1021/jp300623v | J. Phys. Chem. C 2012, 116, 17479−17491

The Journal of Physical Chemistry C

Article

Figure 1. Quartz (101) and water interface model for DFTMD simulation (O atoms, red spheres; Si atoms, gray spheres; H atoms, light blue spheres).

and Obr on the surface using periodic density functional theory molecular dynamics (DFTMD) and classical molecular dynamics (CMD). The effects of Na+, Mg2+, and Cl− on the distribution of these H-bonds is modeled to test the hypothesis that internal surface H-bonds allow for H+-transfer to Obr atoms and that this will increase the Arrhenius pre-exponential factor and dissolution rate.

the bulk and two Si−OH bonds on the surface) whereas the latter consists of Q3 Si−OH. The polymerization state (i.e., number of Si−Obr bonds) is thought to have an effect on the rate of silicate dissolution even though the activation energy remains relatively constant,33 so this difference is significant. Furthermore, a Q2 Si−(OH)2 surface is likely to form a different H-bonding network with water compared to the Q3 Si−OH-dominated surface. The rate constants (k) change dramatically, but the activation energies (ΔEa) are similar for silica dissolution in solutions of increasing ionic strength (I).34 Previous studies of silicate dissolution mechanisms have focused on H+ transfer from solution to bridging oxygen atoms (Obr) on the silicate surface which leads to hydrolysis of Si−O−Si linkages. This mechanism does not explain the ionic strength effect of quartz dissolution because other ions are not involved in the rate-determining step. In fact, one could hypothesize that competition among H+ and other cations would result in a decreased silica dissolution rate with increasing I. Strandh et al.24 discussed the possibility of cation interaction with Obr atoms that would weaken the Si− Obr bonds, but this weakening should result in a lower ΔEa of dissolution. Wallace et al.32 examined potential entropic effects of water reorganization on the Arrhenius pre-exponential factor, but the focus of that study was reorganization of H2O molecules bonded to the adsorbed cation (i.e., Mg2+ and Ca2+) rather than intrasurface H-bonding. In this paper, we focus on exploring the role of intrasurface H-bonding between SiOH groups and bridging O atoms (Obr) on the quartz (101) surface. The role of H-bonding networks at the mineral−water interface can be missed in computational studies. Proton transfers can be enabled through intermediate H2O molecules,35 and H-bonding can at least partially stabilize configurations that would not be predicted in the absence of water.36 Rustad, Casey, and co-workers have demonstrated that reactions on oxide nanoparticles in solution can involve coordinated motions of multiple atoms, and these movements certainly are affected by the H-bonding to each atom on the surface.37−39 The observation that increased I of solutions and the type of salt involved can increase the rate of silica dissolution has also been interpreted as changes in the mineral−water interface (i.e., quartz−water H-bonding arrangements.34,40 Lastly, Washton et al.41 found a correlation between silica gel dissolution rates and the surface concentration of a species thought to be a Q3 SiOH group not H-bonded to water at the interface. This study will focus on the nature of H-bonding between the quartz (101) surface and water and among silanol groups

■

METHODS

Model Building. Although quartz exhibits a number of growth faces and will dissolve at a great variety of surfaces when ground because of its lack of cleavage, the (101) surface was chosen because it is a common crystal growth face and exhibits predominantly Q3 silanol groups (Q3 SiOH) which should be typical of silica surfaces.42,43 From the experimental α-quartz structure, cleavage of the bulk mineral to a six Si3O6-layer slab was achieved using the Surface Builder module of the Materials Studio software package (Accelrys Inc., San Diego, CA). The cleavage plane was adjusted through the model to break a minimum number of bonds and with only Q3 SiO and Q3Si groups exposed. The monoclinic surface simulation cell (4 × 4) was 14.60 Å × 9.82 Å × 39.27 Å with a slab thickness of 19.27 Å and a vacuum thickness of 20 Å (Figure 1). The angle γ between the lattice vectors a and b was 109.7°. Surface hydroxyl groups were bonded to the Q3Si and H atoms were bonded to the Q3SiO surface sites. An energy minimization of this model via relaxation of the ions with lattice parameters held fixed was performed in the Vienna Ab-initio Simulation Package.44,45 A periodic box containing 70 H2O molecules was created in Hyperchem 7 (Hypercube Inc., Gainseville FL). Ten picoseconds of MD simulation were run on this water model in Materials Studio using the Universal Force Field,46 and then the water model was placed within the vacuum gap of the SiO2 periodic slab to accomplish the model [Si72O136(OH)16](H2O)70 system. Energy minimization was carried out; then a DFTMD simulation was run for 11 ps (22 000 time steps of 0.5 fs) with VASP (Figure 2; detail discussed below). To this model structure were added H+/Cl−, Na+, and Mg2+ in separate models. The Na+, Mg2+, and Cl− were placed to random spots within the H2O molecules over surface SiOH silanol groups. One H+ was removed from the nearest silanol group for the Na+ model and 2H+ from two neighboring silanol groups for the Mg2+ model. In the case of the HCl model, the H+ was initially bonded to a random H2O molecule. Energy minimizations were performed in VASP to obtain stable starting configurations followed by MD simulations for at least 11 ps (Figure 2). 17480



Article

dimensions a = 41.259 Å, b = 43.962 Å, c = 54.843 Å. The water density was fixed at 1 g cm−1. The simulations were run for 4 ns, and particle−particle−particle mesh was used to treat the long-range electrostatics.52 The initial configurations of the classical simulations were produced by cleaving the surface slab from a perfect crystal, adding hydroxyl groups to the surface. Water was added at 1 g cm 3 density (equilibrated by running a separate simulation of only SPC water). The atoms in the bulk part of the surface slab were held fixed, but the surface atoms were allowed to relax (see refs 53, 54 for details). The ClayFF force field55 was used for the quartz surface53 which employs Lennard−Jones parameters and partial charges but no special bonding terms except between the H and O in hydroxyls. ClayFF uses the flexible single-point charge (flexible SPC) water model. ClayFF was shown to agree with DFT molecular dynamics with respect to the water structure around the hydroxyl groups (RDFs between water O and hydroxyl H and between water H and hydroxyl O).54 CMD results will be presented for H-bonding between water and hydroxyl H with surface bridging O. For information on the comparison of DFTMD with CMD for other properties see Skelton et al.54

Figure 2. Energy versus time for the α-quartz (101)−H2O model with 500 fs running average (red).

Static calculations of the monolayer-covered (101) surfaces were also performed based on the model developed in our previous work.47 To study the influence of ions on the adsorbed H2O layer structure, single NaOH, Mg(OH)2, or HCl molecules were added to both sides of the slab with the stable dense monolayer of adsorbed H2O molecules. The space for the ions was obtained by moving one or two H2O molecules from the surface to the sites above the ions. More than 10 different starting structures were constructed for each electrolyte. To this purpose, different initial positions of Na+, Mg2+, and Cl− and their counterions were chosen. The slabs with H+ at the Obr site were obtained by manual removal of an H+ from one silanol group or from one of the H2O molecules in the preliminary optimized structures and bonding it to a surface Obr next to SiO−. The lowest energy systems were selected for further analysis. All energy minimizations were made with constant cell parameters obtained in the previous α-quartz crystal structure optimization.47 Slab supercells include four primitive unit cells of the conventional base-centered rectangular lattice for the quartz (101) surface model, and they have the dimensions: a = 14.93, b = 10.06, c = 24.00, γ = 109.7°. C2 symmetry was preserved during all calculations to provide equivalence of two slab surfaces. Computational Methods. VASP was used to calculate all the structures and energies reported in this study.44,45 Calculations were performed using the projector-augmented plane wave method48,49 and the functionals of Perdew−Burke− Ernzerhof (PBE).50 The Si, O, Mg, Na, Cl, and H atom pseudopotentials are as supplied in the VASP library. An energy cutoff of 400 eV, time step of 0.5 fs, NVT ensemble with T scaling via the Nosé thermostat, and a temperature of 300K were used for MD simulations. With the simulation cell parameters fixed, a single k-point sampling was used for MD simulations and 2 × 3 × 1 k-point mesh for the static optimizations. Classical molecular dynamics (CMD) simulations were used to simulate the neutral surface of quartz with water (no ions) using the classical MD code, LAMMPS.51 A 1 fs time step was used with the velocity Verlet integration, Nóse/Hoover thermostat, in the NVT (constant number of molecules N, constant volume V, and constant temperature T) ensemble. There were 2754 atoms in the slab representing the solid and 2320 H2O molecules. An orthorhombic cell was used with

■

RESULTS Energies. In the solvated system, the potential energy leveled off after ≈6 ps of MD simulation and was relatively stable for the last 5 ps of simulation (Figure 2). A similar procedure outlined above for SiO2−H2O was carried out with the addition of a Na+, Mg2+, and H+Cl−. Figure 3 shows the equilibrated structure of the SiO2−H2O−HCl system with the HCl dissociated into H3O+ and Cl−.

Figure 3. α-Quartz (101) and water interface with added H+ and Cl− DFT-MD simulation at 300K and 0.5 fs time step. (See Figure 1 for O, Si, and H labels; Cl− ion, green sphere; H3O+ is shown in stick representation).

The water adsorption energy obtained in our previous study for dense-monolayer coverage was −59 kJ/mol. This value is confirmed in present calculations with the larger (in zdirection) unit cell. Using the static energies of the most stable systems containing adsorbed NaOH, Mg(OH)2, and HCl species, the energy of their adsorption from the gas phase to the water monolayer on the (101) α-quartz surface can be estimated (neglecting the thermal contributions) (Table 1). To this purpose, the energies of isolated NaOH, Mg(OH)2, and HCl species have been calculated using the same unit cell dimensions and symmetry. Data in Table 1 indicated that adsorption energies of NaOH and HCl molecules on the SiO2 17481



Article

the bulk to 145.3° on the surface together with broadening of the distribution toward the greater angles (Figure 4b). The Si−O−Si values for the periodic system and molecular cluster are similar to one another for the protonated Si− (OH+)−Si linkage (Table 2). The main discrepancy is found for the Si−O−Si angle calculated using MP2/6-31G* on (OH)3SiOSi(OH)3.12 This narrower angle is not predicted with other methods such as B3LYP,61 so we conclude that MP2 is over-relaxing the angle which would make the Obr more susceptible to attack by a H+.60 H-Bonding, MD Simulation Results. Examining the (101) surface in Figure 1, there are two types of terminal SiOH groups: one 1 Å farther from the bulk than the other as previously noted.47,53,54 These positions are a result of the crystal structure, not the α-quartz−water interface, but they have implications for the interactions of the surface with water. We refer to these two types as tI and tII, respectively (indicated in Figure 1). In Figure 5a the density profile of the corresponding O atoms is shown for the SiO2−H2O system. All particle distributions were averaged between the two similar surfaces existing in this model. The reference plane (z = 0) for all calculations is the averaged position of the second innermost plane of Si atoms which is the central plane in the top surface layer consisting of three Si-planes and four O-planes. The SiOHtI groups (bound to the mentioned Si atom) at 1.2 Å form weak H-bonds with O atoms on the surface and with O atoms of the H2O molecules over it, and the SiOHtII groups at 2.2 Å form H-bonds only to the H2O molecules in the space between SiO2 slabs. We refer to these weak H-bonds with surface O atoms as “intrasurface” throughout the remainder of this paper. The two SiOH groups also have slightly different mean O−H bond lengths (0.99 and 1.01 Å, respectively). Similar bond lengths were found previously using static calculations in the presence of a dense water monolayer: 0.98 and 1.01 Å for OtI− HtI and OtII−HtII, respectively.47 Figure 5b shows that the O−H RDF for SiOHtI and SiOHtII groups have similar peak values, but the SiOHtI O−H bonds have a larger concentration of shorter O−H bonds whereas the SiOHtII tail at longer O−H distances is slightly higher. The differences in these distributions is likely because the Ow−HtII interaction is more pronounced than Ow−HtI (shown in Skelton et al.53,54), because it protrudes farther into the water, pulling the OtII−HtII to greater bond lengths. Angle distribution functions (ADF) are shown in Figure 5c. The plots in Figure 5c indicate that Si−O−H angles are also

Table 1. Energy Change (kJ/mol) upon Transfer of Considered Electrolytes from Gas Phase to Model Systems species [structure labela]

calculated for the most stable system with adsorbed water monolayer

water solution, standard state, 298 K; experimental estimateb

NaOH [Na_1] Mg(OH)2 [Mg_1] HCl [Cl_1]

−251

−276

−301

−354

−64

−75

Na_1, Na+ and OH− in the adsorption layer; Mg_1, Mg 2+ coordinated to two SiO− groups; Cl_1, Cl− and H3O+ in the adsorption layer. bCalculated using data from ref 56. a

hydroxylated surface covered with water monolayer are less exothermic than corresponding heats of transfer from the gas phase to the water solutions. This can be explained by the relatively small number of H2O molecules surrounding the solvated species in the thin adsorption layer compared to bulk water solution. Structures Si−O Bulk vs Si−O Surface, Si−O−Si Angles. Periodic DFT vs MO Clusters. DFT methods can have difficulty reproducing the observed crystal structure of αquartz;53 but the methods selected in this study reproduce the observed lattice parameters reasonably accurately. Table 2 contains the Si−Obr, SiOH, and H-bond lengths and Si−O−Si angles for bulk α-quartz and the (101) hydrated model surface as calculated in our previous study.47 The values of Si−OH = 1.62−1.64 Å and Si−O−Si = 145−151° are close to the observed values and similar to previous calculations on periodic silica surfaces.57,58 Although the predicted bond lengths are comparable between these periodic DFT calculations and molecular cluster results, the Si−O−Si angles have a larger discrepancy (Table 2). The narrower Si−O−Si angles predicted for the molecular clusters should have a significant effect on calculated reactivity at this site.59,60 The Si−O bonds do not change noticeably with the addition of H2O molecules to the system (Table 2). Although Hbonding is prevalent with the H2O (discussed below), the Si− Obr and Si−OH bond lengths are within 0.01 Å of the hydrated surface. Thus, averaging upon MD trajectory for the solvated system gives 1.63 Å for Si−O bonds both on the surface and in the bulk of the slab (Figure 4a). The Si−O−Si angles change more significantly as the average angle increases from 144.2° in

Table 2. Structural Results of (101) α-Quartz Surface Obtained by Static and MD Simulationsb for Si−O Bonds and Si−O−Si Anglesa species bulk hydrated (101) surface (OH)3SiOSi(OH)3e [Si36O64(OH)16]f solvated (101) surface [Si36O64(OH)16](H2O)24f (OH)3SiOH+Si(OH)3·H2Og [Si36O62(OH)12(ObrH)2O4]{(H2O)13Na}2h

Si−Obr(H), Å

Si−OHtI, Å

Si−OHtII, Å

1.62c (1.63d)

Si−Obr−Si, deg 141c (144d)

1.65 1.62−1.64

1.65 1.64

1.62−1.64 (1.63) 1.76 1.72, 1.82

1.65 (1.64) 1.63 1.55

1.64 1.63 (1.63) 1.59

123 145−151 146−158 (145) 130 132

a

Si−OHtI and Si−OHtII are two types of silanol groups discussed in text. bThe averaged values obtained in MD simulation for [Si72O136(OH)16](H2O)70 are given in parentheses. cOptimized with experimental cell parameters. dFor the inmost layers of the slab in the SiO2−water system (MD result). eXiao and Lasaga12 MP2/6-31G*-optimized cluster results. fFrom ref 47. gThis work using MP2/6-31G* and an Obr−H+ constrained distance of 0.98 Å. hNa+-containing system with protonated bridging oxygen; data are given for bond and angles between the nearest to proton position atoms and for bonds in siloxy SiO− groups; this work (see Discussion below). 17482



Article

Figure 4. Bulk and surface O−Si bond radial distribution function (a) and Si−O−Si angle distribution function (b).

Figure 5. Density profile of hydroxyl O illustrating relative positions of SiOHtI and SiOHtII groups (a), O−H distribution function (RDF) (b), and Si−O−H angle profile (c).

estimated47 SiOtI−HtI bond stretching frequency should be approximately 3700 cm−1 and can be attributed to surface hydroxyls (in considered monolayer) that do not participate in H-bonds as H-donors, consistent with observation for non-Hbonded silanols.63 The lower frequency modes (≈3000 cm−1) are attributed to bond stretching vibrations of SiOtII−HtII in the second Si−OH that participates in H-bonds with H2O molecules, as these have longer O−H bond lengths and a broader distribution (Figure 5b). In the case of ion-added systems, only the atoms belonging to the slab surface next to positions of Cl−, Na+, or Mg2+ were taken into account during the statistical averaging. This reduces

different for two types of hydroxyls; the average value for Si− OtI−HtI is 117° and for Si−OtII−HtII is 121°. A similar value (118°) was obtained in static calculations47 for the Si−OtI−HtI angle, while the monolayer result for Si−OtII−HtII is different: 115°. This can be explained by the stronger influence of the next layer of water on the silanol group SiOHtII, which straightens the angle, compared to the SiOHtI groups that are closer to the SiO2 surface. The O−H bond distances correlate with O−H stretching vibrational frequencies, as observed for hydroxyls in molecules and was found previously for aluminosilicate cluster calculations.61 Based on the O−H distances in these simulations, the 17483



Article

Figure 6. H atom density distribution profiles for SiOHtI (a) and SiOHtII (b).

Figure 7. Radial distribution functions (a) and cumulative coordination numbers (b) for Obr−HtI and Obr−Hw interactions as predicted via DFT and classical MD simulations.

ClayFF parametrization is meant to fit average H-bonding and may not account for stronger H-bonding as reliably as the DFTMD. Addition of Na+, Mg2+, and Cl− all shift the SiOHtII group density profiles by 0.1 to 0.4 A away from the surface, and the tail and shorter distances disappear. The results are consistent with a picture of SiOHtII groups H-bonding more with the water and less with the silica surface as a function of I, which is also consistent with the idea that SiOHtI groups are forming more intrasurface H-bonds. These differences in relative size of the first and second peaks can be interpreted to be a result of increasing the percentage of intrasurface H-bonds at the interface between the α-quartz and ionic solutions. That is, the HtI approach is closer to the surface with addition of ions to the water. To check this conclusion, we have calculated the RDF for surface Obr−HtI and Obr−Hw interactions. The obtained distribution functions and cumulative coordination numbers (N gives the number of neighbors closer than a particular distance) are plotted in Figure 7a and 7b, correspondingly. The percentages of SiOHtI groups involved in intrasurface H-bonding are as follows: H2O 53%, HCl 75%, Na+ 77%, and Mg2+ 67%. The changes are not large enough to explain the observed order of magnitude increase in dissolution rates; however, these percent changes are statisti-

the number of included atoms and increases the statistical errors. Nevertheless, the calculated H atom density profile of the SiOHtI group exhibits a noticeable change relative to that of the SiO2−H2O system (Figure 6). The DFTMD density profile of the SiOHtI groups is bimodal for the SiO2−H2O system, indicating that a noticeable fraction of the intrasurface H-bond groups are converted to SiO2−H2O H-bonds from the hydrated to the solvated MD model (Figure 6a). CMD results on SiO2− H2O system are also included in Figure 6 and exhibit bimodal structure for HtI atoms similar to the DFTMD results except that the CMD density peak at ∼1 Å is smaller than that of DFTMD. The H atom density profile of the SiOHtI group becomes more unimodal on the surface nearest the Na+ as the second peak near 2.2 Å almost disappears. In the presence of Mg2+, the second peak persists, but it is reduced in intensity and shifted toward smaller distances from the surface. The density profile of the SiOHtI groups for the Cl−-containing system resembles that of the SiO2−H2O system with a reverse ratio of peak heights. The density profiles of the SiOHtII groups are unimodal (Figure 6b). The DFTMD and CMD results agree well in the case of pure water except that the CMD does not produce a tail at shorter distances. This is likely due to the fact that the SPC/ 17484



Article

Figure 8. Ion−water radial distribution functions. Na+−O RDFs (a), Mg2+−O RDFs (b), and Cl−−H RDFs (c). Lines are cumulative coordination numbers N.

density peak in Figure 6a in the CMD results is correlated with the lowering of the Obr−HtI RDF (Figure 7a). The obtained DFTMD RDFs show that the presence of ions increases the probability of H-bond formation between the HtI and surface Obr. The effect is most prominent in the case of Na+ and smaller in the case of Mg2+, which is in accordance with the density profile analysis result. This indicates that all the SiOHtI groups in the system with added NaOH are H-bonding to other atoms on the surface and not to H2O molecules, consistent with the trend in dissolution rates upon addition of NaCl and MgCl2 to solution-dissolving quartz.34 That is, the presence of Na+ and Mg2+ ions increased the H-bonding of HtI to Obr, which could increase the rate of dissolution. The current interpretation is in contrast to the model results of Strandh et al.24 who concluded that reorganization of H2O molecules around the adsorbed cation could contribute to the activation entropy changes (i.e., changes in A) observed for silica dissolution. As expected, H2O molecules are unable to strongly interact with surface Obr atoms. As can be seen on Obr−Hw functions presented in Figure 7b, H2O molecules do not approach Obr closer than 2.2 Å. DFTMD shows a small increase of corresponding H-bonds in Na+- and Mg2+-containing systems. Again, there are differences between the CMD and DFTMD, with CMD showing H-bonds at greater distances for Obr−HtI

cally limited and would benefit from simulations of larger spatial extent and longer temporal duration. Furthermore, examination of other faces will be critical. However, such DFTMD simulations are not practical, so classical MD simulations should be performed to test the effects of these factors. In addition, the cations may play a role in stabilizing the hydrolysis precursor species Si−(OH+)−Si (see Influence of Na+ and Mg2+ on the Basicity of Obr: Static Simulation Results). Due to geometrical (sterical) constraints, OHtI can form only long, weak H-bonds with the Obr. Because of this, there is no maximum in the corresponding RDF in the conventional region 1.8 to 2.0 Å. Also shown in Figure 7 are results of the CMD simulation of water in contact with quartz for H2O and HtI. There are some discrepancies between the DFTMD and CMD. The Obr−HtI interaction is weaker for CMD than for DFTMD as the density becomes nonzero at greater distances. Furthermore, the shape of the curve is different, with the CMD result showing a maximum at ∼2.8 Å. The Obr−Hw curve, however, shows a stronger interaction with the CMD curve becoming nonzero at closer distances than with DFTMD. Moreover, it rises to greater values at greater distances. The CMD result reinforces the argument that the height of the first peak (at 1 Å) in Figure 6a correlates with the intensity of interaction between HtI with Obr. The reduction of the first 17485



Article

and at closer distances for Obr−Hw. There are several reasons for discrepancies between DFTMD and CMD. First, and probably the most important reason, is that the ClayFF force field has not been parametrized for Obr−Hw or Obr−HtI interactions, and this could affect the interaction. Second, there are slight differences in the surface structure between DFTMD and CMD,54 and this can highly influence the accessibility of Hw and HtI to the Obr groups. In preliminary simulations, the surface was fixed at a different configuration (results not shown), yielding a closer interaction of both Hw and HtI to the Obr. Interestingly, Skelton et al.54 showed that both DFTMD and CMD give different surface relaxations compared to X-ray reflectivity experiments, and it is unclear which is producing the correct surface structure. The last reason for the discrepancy between DFTMD and CMD could be the sampling because the CMD simulations are larger and longer. To test this, CMD simulations were repeated for small systems (13.755 × 9.829 × 41.5 Å) for 10 ps to obtain statistics similar to those of the DFTMD simulations. The Obr− Hw and Obr−HtI RDFs are shown in Supporting Information (Figure S1). The Obr−HtI RDF shows a slight difference for the large and small systems but not enough of a difference to explain the differences between the DFTMD and CMD. Conversely, the difference between large and small CMD simulations for the Obr−Hw interaction (the small simulation showing a closer interaction) could be significant when comparing the DFTMD and CMD. To unequivocally rule out the effect of sampling, the DFTMD simulations should be run for larger time and length scales, but this is presently impractical. DFTMD and CMD showed a good comparison for the strong H-bonded interactions between hydroxyls and water,54 and one reason for the larger discrepancy for the HtI−Obr RDF is that the statistics are poor. Therefore, it is likely that H-bond distributions are sensitive to the factors discussed above. Positions of Ions and Ion−Water Radial Distribution Functions. Radial distribution functions for ion−water oxygen interactions are given in Figure 8 for Na+-, Mg2+-, and Cl−containing systems. The cumulative coordination number N(r) is also plotted in these figures (label on right y-axis). Because of slow translational diffusion, Na+ and Mg2+ cations stay in the vicinity of the deprotonated siloxy OtII atom where they were placed initially during most of the simulation time. Figure 8a (the value of N) shows that Na+ coordinates with three H2O O atoms and one surface O atom, thus having a total coordination number 4 which is less than in bulk water. The low coordination numbers of adsorbed cation can be explained by the geometric constraints for the surface side. Also, the surface hydroxyls can obstruct the orientations of H2O molecules that would be favorable for the interactions with cations. The averaged distances between Na+ and O atoms of both types are close to 2.4 Å. Only traces of OH− anion were detected in the MD trajectory (see below), so the charge compensation in simulated structures is coming from the deprotonation. However, in the case of Mg2+ (Figure 8b), one OH− is present in the vicinity of the cation, which is also bonded to one siloxy O atom. The second initially created siloxy O atom takes a H+ from one of the H2O molecules which produces OH−. The average distance from Mg2+ to the siloxy atom OtI is noticeably shorter (1.9 Å) than distances to O in OH− and H2O (2.0 Å). The total coordination number of Mg2+ in the sampled structures is the same as for Na+ (i.e., 4).

HCl exists in the dissociated form during the entire simulation. Cl− is bound by a H-bond to silanol SiOHtI and surrounded by three to four H2O molecules (Figure 8c) that also interact with Cl− by means of H-bonds. The first type of H-bond is slightly shorter (2.1 Å) than the second (2.2 Å). Species Concentrations. The calculated average water density in the space between the slab surfaces is close to that in ambient conditions (1.0 g cm−1). To estimate the concentrations of different ionic species, we considered 1.2 Å as an upper bound for O−H bonds. The SiOH group or H2O molecule was considered dissociated if the H+ moves away from the corresponding O atom farther than this limit. If a H2O or silanol O atom has an extra H+ at a distance closer than 1.2 Å, we suggest the formation of H3O+ or SiOH2+. Such estimation is a first-order approximation due to uncertainties in the upper limit of vibrating O−H bond length. Nevertheless, the statistical treatment of obtained MD results, based on such an approach, gives evidence that a dynamical equilibrium may exist between the “neutral” and “charged” species in the thin water layer (of about 0.7 Å) adjacent to the surface. We found that silanol groups SiOHtI stay in the associative state in all considered systems, except perhaps the Mg2+-containing system, where a negligible number of SiOHtI dissociation events were observed during MD simulation. The SiOHtII group has a greater tendency to dissociate, probably because the HtII−Ow Hbonding is stronger than the HtI−Ow H-bonding.54 In the water−SiO2 system, approximately 0.1% of these silanols are reversibly dissociated. One SiOHtII group stays in the dissociative state during most of the simulation in both Na+ and Mg2+ models. Only a small amount of “free” OH− was detected in the case of Na+, whereas in the case of Mg2+, the average numbers of OH− and SiO− siloxy groups coordinated with Mg2+ are 0.9 and 0.1, respectively. This is in accordance with the fact that hydrolysis of Mg2+ in water solution is much stronger than that of Na+. In the HCl-containing system, the only charged species were Cl− and H+ (H3O+). No formation of SiOH2+ groups was found in any considered systems. Influence of Na+ and Mg2+on the Basicity of Obr: Static Simulation Results. MD simulation provides evidence that ionic species can enhance intrasurface H-bonding. However, because of low translational diffusion, the ions remain near the same surface group and do not move considerably during the relatively short MD run. For this reason, we performed static DFT calculations on systems containing ions in the dense water monolayer on the (101) quartz surface with the ions placed at various locations near both types of silanol groups. Examining the different initial positions of ions in the adsorbed water layer, we have determined the lowest energy ion environments and found the conditions that favor protonation of the Obr atoms. No such stable structures were obtained previously in our study47 of water adsorption on hydroxylated (101) quartz surface in the absence of ions. All initial (dissociated) structures considered in ref 47 with protonated (at the surface level or above) Obr quickly transformed into the structures with neutral adsorbed H2O molecules via the Grotthuss H+-transfer mechanism. However, in this study the several metastable structures with protonated bridging hydroxyl (i.e., Si−(OH+)− Si) are established in the systems containing Na+ and Mg2+ ions and discussed below. This indicates an increase in Obr basicity in the presence of cations. We conclude that this is due to the extra negative charge on the Obr atom, arising because of the ionization of the silanol groups. Our results show that cations cannot directly interact with Obr because of strong repulsion 17486



Article

Table 3. Relative Energiesa of Selected Models of Fully Hydrolyzed α-Quartz Surface Covered with Dense Water Monolayer with Adsorbed Ions (the energies of unstable intermediate states are given in italics) label

structurea

Na_1

[Si36O64(OH)16] {(H2O)12Na(OH)}2 [Si36O64(OH)14O2]{(H2O)13Na}2 [Si36O62(OH)12(ObrH)2O4] {(H2O)13Na}2 [Si36O64(OH)10 (OH2)2O4] {(H2O)13Na}2 [Si36O64(OH)12O4]{(H2O)14Mg}2 [Si36O64(OH)14O2] {(H2O)13Mg(OH)}2 [Si36O62(OH)10(ObrH)2O6] {(H2O)14Mg}2 [Si36O64(OH)8 (OH2)2O6] {(H2O)14Mg}2 [Si36O64(OH)16]{(H2O)11(H3O) Cl}2 [Si36O64(OH)16]{(H2O)11(H3O) Cl}2 [Si36O62(OH)14(ObrH)2O2] {(H2O)11(H3O)Cl}2 [Si36O64(OH)16]{(H2O)12(HCl)}2

Na_2 Na_3 Na_4 Mg_1 Mg_2 Mg_3 Mg_4 Cl_1 Cl_2 Cl_3 Cl_4

ΔEb, kJ/mol

comment Na+ coordinated to two siloxy groups Si−OHtI and Si−OHtII; OH− in the adsorption layer

0

Na+ coordinated to SiO− tII group Na+ coordinated to two siloxy groups SiO−tI and SiO−tII; protonated bridging oxygen, Si−(OH+)− Si Na+ coordinated to two siloxy groups SiO−tI and SiO−tII; protonated silanol oxygenc, Si−OH2+

3 134 (95)

Mg2+ coordinated to two siloxy groups SiO−tI and SiO−tII Mg2+ coordinated to one SiO− tII group; OH− in the adsorption layer

0 14

Mg2+ coordinated to two siloxy groups SiO−tI and SiO−tII; protonated bridging oxygen, Si− (OH+)−Si Mg2+ coordinated to two siloxy groups SiO−tI and SiO−tII; protonated silanol oxygenc, Si−OH2+

112 (76)

Cl− H-bonded to two silanols SiOHtII; H3O

0

Cl− H-bonded to silanol SiOHtI; H3O

+

+

in the adsorption layer

in the adsorption layer

138 (99)

108 (72)

25

protonated bridging oxygend, Si−(OH+)−Si

130

HCl molecule in the adsorption layer

69

The chemically interconnected atoms of the “solid” phase are enclosed in square brackets, whereas the atoms of both adsorbed layers are enclosed in braces. bThe energy of relaxed structure obtained after the manual proton transfer to the nearest silanol oxygen is given in parentheses. cUnstable intermediate with protonated silanol O (optimized with restraints on proton movement). dUnstable intermediate with protonated Obr (optimized with restraints on proton movement).

a

Figure 9. The most stable structures for Na+ adsorbed in the water monolayer. Na_1, bonding to silanols SiOHtI and SiOHtII (a), and Na_2, bonding to siloxy group SiO−tII (b). (See Figure 1 for O, Si, and H labels; Na+ ion, blue sphere; participating silanols and hydroxyl ion are shown as sticks).

Na_2 where Na+ is over the siloxy OtII atom (Figure 9b). This position of Na+ is close to that simulated in MD study. Both structures, Na_1 and Na_2, have almost the same energy. In the case of Na_3 with a bridging hydroxyl, Na+ is located between two siloxy groups at separations 2.2 Å from the OtI atom and 2.6 Å from the OtII atom. The bridging H+ is connected to the nearest (opposite) silanol O atom via a short (1.6 Å) H-bond (Figure 10a). The SiObr bond length between siloxy Si and Obr atoms is considerably increased compared to the normal SiO bonds on the surface (1.82 vs 1.65 Å). To verify that the stability of the bridging H+ site is a result of the presence of Na+, we explicitly replaced Na+ with H+ at the same position. Reoptimization of the obtained structure initiates the transformation of the nearest siloxy groups to silanol groups and finally results in removal of the H+ from the bridging site. On the other hand, manual transfer of H+ from the bridging position to the nearest silanol O atom leads to an

from the next positively charged Si atoms (the minimum separation between them is about 3.3 Å). In Table 3, the energies of the selected structures are provided. The zeroenergy level corresponds to the most stable structures (Na_1, Mg_1, and Cl_1) found in each case. All mentioned structures are presented in Supporting Information in PDB format. In the cases of Na+-containing systems, the coordination number of Na+ is between 4 to 6. Coordination greater than 4 is possible when Na+ is located closer to silanol SiOHtI, whereas Na+ stays near the SiOHtII in MD simulated structures. For the same reason, the water/silanol O−Na+ distances may be larger than those found in MD simulations; they vary from 2.4 to 2.7 Å. In the structure labeled as Na_1 (Table 3, Figure 9a), Na+ is located above the two silanol O and not directly bonded to any siloxy group; charge compensation is occurring via the presence of OH− in the vicinity of the cation. Another charge compensation mechanism is represented by the structure 17487



Article

Figure 10. Structure of the metastable protonated bridging site. Na_3 structure in the presence of Na+ (a) and Mg_3 structure in the presence of Mg2+ (b). Participating siloxy groups, silanols, and hydroxyls are shown as sticks; other atoms and bonds are shown as lines. The H2O molecule that stabilizes the terminal silanol tI is shown in ball-and-sticks representation.

Figure 11. The most stable structures for Mg2+ adsorbed in the water monolayer. Mg_1, bonding to two siloxy groups SiO−tI and SiO−tII (a), and Mg_2, bonding to siloxy group SiO−tII (b). (Mg2+ ion, violet sphere; participating silanols and hydroxyl ion are shown as sticks).

Figure 12. The most stable structures for Cl− adsorbed in the water monolayer. Cl_1, H-bonding to two silanols SiOHtII (a), and Cl_2, H-bonding to silanol SiOHtI (b). (See Figure 3 for O, Si, H, and Cl− labels; participating silanols, hydroxyls, and hydronium ion are shown as sticks).

Addition of Mg(OH)2 to the adsorbed water layer results in changes similar to those observed in the Na+-containing systems. As expected, Mg2+ interacts with surface groups more strongly than Na+ cation. In the obtained structures, Mg2+ predominantly bonded to one or two siloxy O atoms. In the second case (structure Mg_1, Figure 11a), Mg2+ was located over SiOtI− and SiOtII− groups at equal distances of 1.9 and 2.0 Å from the both O atoms. In the second case Mg_2, Mg2+ was bonded to one SiOtII− and OH− at the same separation (Figure 11b). The remaining sites are occupied by H2O and/or silanol

unstable SiOH2 group that quickly relaxes via the Grotthuss H+-transfer mechanism to a more stable state which is −95 kJ/ mol lower in energy than the structure with the bridging hydroxyl. If the movement of upper H+ in intermediate OH2 group is constrained, the partially optimized structure (Na_4) can be obtained with almost the same energy (Table 3) as found for the structure Na_3. This indicates that silanol and Obr have similar proton affinities, although no stable protonated silanol (i.e., SiOH2+) could be obtained in this work. 17488



Article

interpretation of experimental data of α-quartz in salt solutions made by Karlsson et al.67 that the acid−base properties of the silica surface change with electrolytes. The similarities of computed bond lengths between molecular cluster and these periodic DFT calculations (Table 2) suggests that the molecular cluster approximation was reasonably justified in earlier studies as a practical method of simulating Si−O−Si hydrolysis. The assumption that lattice effects would significantly hinder atomic relaxations as the Obr became protonated appears to have been overestimated. Although lattice effects and solvation were taken into account in this study, calculated bond lengths and Si−O−Si angles are approximately equivalent. When one considers that the planewave basis used here and the atomic orbital basis sets used in molecular cluster calculations also influence the results, the effects of system size on structure are minor. Furthermore, the energetics of H+-transfer to Obr from the periodic DFT and molecular cluster calculations are similar.17 This fact suggests that it is the short-range structure and strength of the individual Si−O bonds that control this hydrolysis reaction. However, the H+-transfer reaction is strongly influenced by the solvation at the interface, which is the main component missing from many earlier computational studies of silica dissolution. The conclusions drawn from this study stand in contrast, however, to previous work on cation effects on quartz dissolution24,32 in that we do not ascribe the change in dissolution rate to a weakening of the Si−O−Si linkage and subsequent lowering of the activation energy barrier. The entropic effects discussed in Wallace et al.32 appear to be dominant in this case. Comparison to Experiment. Another effect not considered previously in modeling studies, however, is the possible increase in the pre-exponential factor A of the Arrhenius expression (k = A exp(−ΔEa/RT)). Most computational studies of silicate dissolution or hydrolysis have focused on the ΔEa term. In this study, the noticeable increase in the percentage of intrasurface H-bonds on the surface with Cl− and Na+ in proximity (Figure 7) is consistent with the idea that H+transfer can occur more readily through these H-bonds such that a higher concentration of intrasurface H-bonds would lead to a larger value of A. This effect also may be influenced by the topography of the surface if step defects are the dominant site for Si−O−Si hydrolysis. Both of these model results are consistent with the experimental data and interpretations of Dove and co-workers. First, Dove34 found that ionic strength and the type of electrolyte strongly influenced dissolution rates of quartz without affect ΔEa significantly. Second, Dove34 also hypothesized that hydrated cations at the silica−water interface might increase the accessibility of H2O to Si−O bonds. The results in this study suggest a modification of this interpretation in that the intrasurface SiOH groups have an increased accessibility to forming H-bonds and possibly Si−(OH+)−Si groups when cations are in the silica−water interface rather than H2O molecules directly interacting with Obr atoms. Washton et al.41 concluded that non-H-bonded Q3 groups are the rate-controlling species for silica dissolution. These DFT results show increased intrasurface H-bonding with inclusion of Na+ and Mg2+ that is consistent with this interpretation and the dissolution rate data.34 If adding salt to the solution increases the number of intrasurface H-bonds (i.e., Q3 SiOH group not H-bonded to H2O), then the rate can increase due to an increased pre-exponential factor A. This interpretation is generally consistent with the model of Dove and co-workers for crystalline and amorphous silica dissolu-

O atoms with Mg−O distances of 2.0 to 2.4 Å. The total coordination number of Mg2+ in the considered adsorption layers is generally 4, but it can be 5 if the cation is bonded to SiOtI− or it is located between SiOtI− and SiOtII−. The structure with protonated Obr (Mg_3, Figure 10b) has one extra siloxy group to retain the overall composition. As for the Na+ case, the bridging H atom is connected to the nearest silanol O atoms via a short (1.6 Å) H-bond, and the Si(siloxy)−Obr bond is rather long, 1.79 Å. Again, the constrained intermediate structure with the OH2 group obtained by the manual H+-transfer from the bridging position to the nearest silanol O atom has an energy close to that of structure with the bridging hydroxyl, thus indicating the parity of the proton affinities for silanol and Obr in the considered systems. This result is qualitatively similar to the effects observed in Wallace et al.32 that examined the effect of alkali cations on the electronic state of the Obr. No stable bridging states were found for H+ in the presence of HCl in the adsorption layer. The most favored structures (Cl_1 and similar) contain Cl− above and between the two terminal SiOHtII groups. Silanol H atoms form H-bonds with Cl− of length 2.0−2.2 Å (Figure 12a). Such reorientation of the terminal surface hydroxyl groups associated with Cl− was observed in a classical MD simulation of aqueous electrolyte solutions within slit-shaped silica nanopores by Argyris et al.64 In other (slightly less favorable) structures (Cl_2 and similar, Figure 12b), Cl− is connected via a H-bond to the H of an SiOHtI group as it was simulated by DFTMD (see Figure 3). H3O+ provides the charge compensation in all those structures. In the unstable intermediate structure Cl_3 (which corresponds to the flat shadow on the energy optimization curve for the structure Cl_1), H3O+ is bonded to the siloxy SiOtI group. This position of H3O+ is similar to the Na+ position in the structure Na_3, which presumably provides the possibility of Obr protonation, although this structure is unstable and finally relaxes to Cl_1 by H+-transfer through the Grotthuss mechanism. One structure (Cl_4) was found with HCl adsorbed in the molecular form. The energy is +70 kJ/mol higher than the energy of the most stable structure. Hence, one can conclude that HCl dissociates entirely even within the thin water monolayer adsorbed on the silica surface.

■

DISCUSSION Comparison of Periodic DFT Results to Previous Results. The H-bond interaction of H2O with Obr atoms is less energetically favorable than H2O−H2O H-bonding65 which is reflected in the observation that the current MD simulations reveal SiOH−H2O H-bonding but not Obr−H2O H-bonding (Figure 7b). However, although a minor percentage, SiOH groups were observed to rotate and weakly H-bond to Obr instead of H2O molecules during the simulations (Figures 6 and 7). Because the proton affinity of Si−OH2+ groups may be similar to that of Si−(OH+)−Si groups, as shown previously in cluster calculations66 and confirmed in the present static optimization of the adsorbed layer, the possibility of transferring a H+ from a SiOH group to the Obr can be suggested as the mechanism for forming Si−(OH+)−Si rather than directly transferring H+ from H3O+ to a Obr (e.g., refs 12, 17). Nangia and Garrison25 have also suggested that there is a critical role for intrasurface H-bonding in the dissolution of the SiO2 phase cristobalite. In this work, we also have provided evidence that Na+ and Mg2+ cations coordinated by siloxy groups can increase the basicity of the bridging oxygen, thus providing local stability of the protonated bridging site. This conclusion is similar to the 17489



■

tion.6,68 However, the results of this study do not show that reorganization of the H2O molecules or formation of OH− at the interface due to the presence of Na+ or Mg2+ will lead to increased dissolution rates32,34

■

Article

ASSOCIATED CONTENT

S Supporting Information *

Atomic coordinates for the models shown in the figures of this paper in Protein DataBank (*.pdb) format: si_001, Figure 1; si_002, Figure 3; si_004, Mg_1; si_005, Mg_2; si_006, Mg_3; si_007, Mg_4; si_008, Na_1; si_009, Na_2; si_011, Na_3; si_012, Na_4; si_013, Cl_1; si_014, Cl_2; si_015, Cl_3; si_016, Cl_4; si_017 and si_018, the SiO2−H2O−Mg2+ and SiO2−H2O−Na+ models discussed in Figures 6−8; si_021, the Obr− 465 Hw and Obr−HtI RDFs (Figure S1). This material is available free of charge via the Internet at http://pubs.acs.org.

CONCLUSIONS

The hypothesis that silicate dissolution occurs through protonation and hydrolysis of Obr atoms mediated via intrasurface H-bonding has been proposed. The hypothesis is supported by experimental evidence showing increased dissolution rates with increasing ionic strength.34 This experimental data show that ΔEas are nearly the same, so it is reasonable to propose that the Arrhenius pre-exponential factor, A, is the main source of the more rapid dissolution rates in the salt solutions. Because A is related to the activation entropy of a reaction, any phenomenon that preorganizes the system to be more predisposed toward reaction will serve to increase the rate. Observed correlations between silicate dissolution rates and adsorbed 3,3,3-(trifluoropropyl)dimethylchlorosilane (TFS) concentrations41 also support this hypothesis. TFS is thought to bond to isolated Q3 SiOH groups (i.e., silanols not H-bonded to water69), and it is these types of groups that have been proposed to be responsible for H+transfer to Obr atoms. The results of this work predict that the intrasurface Hbonding will increase as cations such as Na+ and Mg2+ are added to solution. In the silica−water interface region, the reordering of the water H-bonding network that occurs when these cations are near the surface disrupts the ability of surface silanols to form SiOH−H2O H-bonding, freeing the SiOH to form intrasurface H-bonds. Furthermore, the effect of the Na+ is predicted to be greater than that of Mg2+ in this phenomenon which is consistent with the greater increase in dissolution rate observed for NaCl over MgCl2.34 Our results also suggest that the basicity of the Obr atoms changes with addition of Na+ and Mg2+ cations, which will allow for partial stabilization of the Si− (OH+)−Si linkage.67 The latter is the suggested precursor to silicate dissolution because hydrolysis of these metastable bonds requires much less energy than that for the normal Si− O−Si linkage.11 The paradigm for oxide protonation behavior has long held that the concentration of M−OH2+ groups equals the concentration of M−O− groups at the point-of-zero charge.70 Periodic density functional theory calculations presented from this study confirm the earlier suggestion66 that the proton affinity of Obr atoms is approximately the same as SiOH groups. Consequently, if one can suggest that SiOH2+ groups begin to form on silicate surfaces, then there is no reason why Si−(OH+)−Si groups should not form as well. Our calculations demonstrate the feasibility of these metastable states and a possible mechanism for their formation. We regard these statements to be in the hypothesis stage because the DFTMD simulations are limited in spatial, temporal, and compositional extent. Extending the simulations in all three dimensions would be a useful test of this hypothesis. Although the hypothesis is borne out of experimental observations cited above, studies using spectroscopic or scattering techniques that could search for the type of intrasurface H-bonds we have proposed would be a good test of our hypothesis.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS Research sponsored by the Division of Chemical Sciences, Geosciences and Biosciences, Office of Basic Energy Sciences, U.S. Department of Energy. Computation was supported by the Research Computing and Cyberinfrastructure (RCC), a unit of The Pennsylvania State University Information and Technology Services facility. The authors also appreciate discussions and comments from Patricia Dove, Peter Cummings, Nitin Kumar, and David Wesolowski.

■

REFERENCES

(1) Dove, P. M. Rev. Mineral. 1995, 31, 235−290. (2) Klause, M.; Rothhaar, U.; Bicker, M.; Ohling, W. J. Non-Cryst. Sol. 2010, 356, 141−146. (3) Rani, N.; Shrivastava, J. P.; Bajpai, R. K. Curr. Sci. 2010, 98, 950− 954. (4) Frugier, P.; Gin, S.; Minet, Y.; Chave, T.; Bonin, B.; Godon, N.; Lartigue, J. E.; Jollivet, P.; Ayral, A.; De Windt, L.; Santarini, G. J. Nucl. Mater. 2008, 380, 8−21. (5) Hamilton, J. P.; Pantano, C. G.; Brantley, S. L. Geochim. Cosmochim. Acta 2000, 64, 2603−2615. (6) Dove, P. M.; Han, N.; Wallace, A. F.; De Yoreo, J. J. Proc. Natl. Acad. Sci. U.S.A. 2008, 105, 9903−9908. (7) Rimstidt, J. D.; Barnes, H. L. Geochim. Cosmochim. Acta 1980, 44, 1683−1699. (8) Dove, P. M.; Han, N.; De Yoreo, J. J. Proc. Nat. Acad. Sci. U.S.A. 2005, 102, 15357−15362. (9) Zhang, L.; Luttge, A. Geochim. Cosmochim. Acta 2009, 73, 6757− 6770. (10) Arvidson, R.; Luttge, A. Chem. Geol. 2010, 269, 79−88. (11) Casey, W. H.; Lasaga, A. C.; Gibbs, G. V. Geochim. Cosmochim. Acta 1990, 54, 3369−3378. (12) Xiao, Y. T.; Lasaga, A. C. Geochim. Cosmochim. Acta 1994, 58, 5379−5400. (13) Xiao, Y. T.; Lasaga, A. C. Geochim. Cosmochim. Acta 1996, 60, 2283−2295. (14) Pelmenschikov, A.; Strandh, H.; Pettersson, L. G. M.; Leszczynski, J. J. Phys. Chem. B 2000, 104, 5779−5783. (15) Pelmenschikov, A.; Leszczynski, J.; Pettersson, L. G. M. J. Phys. Chem. A 2001, 105, 9528−9532. (16) de Leeuw, N. H. J. Phys. Chem. B 2001, 105, 9747−9754. (17) Criscenti, L. J.; Kubicki, J. D.; Brantley, S. L. J. Phys. Chem. A 2006, 110, 198−206. (18) Nangia, S.; Garrison, B. J. J. Phys. Chem. A 2008, 112, 2027− 2033. Nangia, S.; Garrison, B. J. J. Phys. Chem. C 2010, 114, 2267− 2272. 17490



Article

(19) Adeagbo, W. A.; Doltsinis, N. L.; Klevakina, K.; Renner, J. Chem. Phys. Chem. 2008, 9, 994−1002. (20) Lopes, P. E. M.; Murashov, V.; Tazi, M.; Demchuk, E.; MacKerell, A. D. J. Phys. Chem. B 2006, 110, 2782−2792. (21) Du, Z.; de Leeuw, N. H. Dalton Trans. 2006, 2623−2634. (22) Mahadevan, T. S.; Garofalini, S. H. J. Phys. Chem. C 2008, 112, 1507−1515. (23) Fogarty, J. C.; Aktulga, H. M.; Grama, A. Y.; van Duin, A. C. T.; Pandit, S. A. J. Chem. Phys. 2010, 132, 174704/1−174704/10. (24) Strandh, H.; Pettersson, L. G. M.; Sjoberg, L.; Wahlgren, U. Geochim. Cosmochim. Acta 1997, 61, 2577−2587. (25) Nangia, S.; Garrison, B. J. Mol. Phys. 2009, 107, 831−843. (26) Du, Z. M.; de Leeuw, N. H. Surf. Sci. 2004, 554, 193−210. (27) Murashov, V. V. J. Phys. Chem. B 2005, 109, 4144−4151. Murashov, V. V.; Demchuk, E. J. Phys. Chem. B 2005, 109, 10835− 10841. (28) Yang, J.; Wang, E. G. Phys. Rev. B 2006, 73, 035406/1−035406/ 7. (29) Goumans, T. P. M.; Wander, A.; Brown, W. A.; Catlow, C. R. A. Phys. Chem. Chem. Phys. 2007, 9, 2146−2152. (30) Lopes, P. E. M.; Demchuk, E.; Mackerell, A. D., Jr. Int. J. Quantum Chem. 2009, 109, 50−64. (31) Leung, K.; Nielsen, I. M. B.; Criscenti, L. J. J. Am. Chem. Soc. 2009, 131, 18358−18365. (32) Wallace, A. F.; Gibbs, G. V.; Dove, P. M. J. Phys. Chem. A 2010, 114, 2534−2542. (33) Bandstra, J. Z.; Brantley, S. L. Geochim. Cosmochim. Acta 2008, 72, 2587−2600. (34) Dove, P. M. Geochim. Cosmochim. Acta 1999, 63, 3715−3727. (35) Isaev, A. N. Russ. J. Phys. Chem. A 2010, 84, 434−443. (36) Catlow, C. R. A.; Bromley, S. T.; Hamad, S.; Mora-Fonz, M.; Sokol, A. A.; Woodley, S. M. Phys. Chem. Chem. Phys. 2010, 12, 786− 811. (37) Rustad, J. R.; Hay, B. P. Geochim. Cosmochim. Acta 1995, 59, 1251−1257. Rustad, J. R.; Wasserman, E.; Felmy, A. R.; Wilke, C. J. Colloid Interface Sci. 1998, 198, 119−129. (38) Casey, W. H.; Rustad, J. R.; Spiccia, L. Chem.Eur. J. 2009, 15, 4496−4515. (39) Wander, M. C. F.; Rustad, J. R.; Casey, W. H. J. Phys. Chem. A 2010, 114, 1917−1925. (40) Icenhower, J. P.; Dove, P. M. Geochim. Cosmochim. Acta 2000, 64, 4193−4203. (41) Washton, N. M.; Brantley, S. L.; Mueller, K. T. Geochim. Cosmochim. Acta 2008, 72, 5949−5961. (42) Carroll, S. A.; Maxwell, R. S.; Bourcier, W.; Martin, S.; Hulsey, S. Geochim. Cosmochim. Acta 2002, 66, 913−926. (43) Nangia, S.; Washton, N. M.; Mueller, K. T.; Kubicki, J. D.; Garrison, B. J. J. Phys. Chem. C 2007, 111, 5169−5177. (44) Kresse, G.; Furthmüller, J. Phys. Rev. B 1996, 54, 11169−11186. (45) Kresse, G.; Furthmüller, J. Vasp the Guide; Institut für Materialphysik, Universität Wien: Vienna, Austria, 2003. (46) Rappé, A. K.; Casewit, C. J.; Colwell, K. S.; Goddard, W. A.; Skiff, W. M. J. Am. Chem. Soc. 1992, 114, 10024−10035. (47) Bandura, A. V.; Kubicki, J. D.; Sofo, J. O. J. Phys. Chem. C 2011, 115, 5756−5766. (48) Blöchl, P. E. Phys. Rev. B 1994, 50, 17953−17979. (49) Kresse, G.; Joubert, D. Phys. Rev. B 1999, 59, 1758−1775. (50) Perdew, J. P.; Burke, K.; Ernzerhof, M. Phys. Rev. Lett. 1996, 77, 3865−3868. (51) Plimpton, S. J. Comp. Phys. 1995, 117, 1−19. (52) Hockney, R. W.; Eastwood, J. W. As a Simple Demonstration of the Possible Advantage of the Computer Simulation Using Particles; IOP Publishing Ltd.: Bristol, England, 1988. (53) Skelton, A. A.; Fenter, P.; Kubicki, J. D.; Wesolowski, D. J.; Cummings, P. T. J. Phys. Chem. C 2011, 115, 2076−2088. (54) Skelton, A. A.; Wesolowski, D. J.; Cummings, P. T. Langmuir 2011, 27, 8700−8709. (55) Cygan, R. T.; Liang, J. J.; Kalinichev, A. G. J. Phys. Chem. B 2004, 108, 1255−1266.

(56) (a) Chase, M. W. J. Phys. Chem. Ref. Data, Monogr. 1998, 9, 1− 1951 (JANAF Thermochemical Tables, 4th ed.). (b) Glushko, V. P.; Medvedev, V. A.; Gurvich, L. V., Eds. Thermal Constants of Substances; Wiley: New York, 1999. (57) Demichelis, R.; Civalleri, B.; Ferrabone, M.; Dovesi, R. Int. J. Quantum Chem. 2010, 110, 406−415. (58) Tielens, F.; Gervais, C.; Lambert, J. F.; Mauri, F.; Costa, D. Chem. Mater. 2008, 20, 3336−3344. (59) Lockwood, G. K.; Garofalini, S. H. J. Chem. Phys. 2009, 131, 074703/1−074703/8. (60) Gibbs, G. V.; Wallace, A. F.; Cox, D. F.; Downs, R. T.; Ross, N. L.; Rosso, K. M. Am. Mineral. 2009, 94, 1085−1102. (61) Kubicki, J. D.; Sykes, D.; Rossman, G. R. Phys. Chem. Min. 1993, 20, 425−432. (62) Novak, A. In Structure and Bonding 18: Large Molecules; Dunitz, J. D.; Hemmerich, P.; Holm, R. H.; Ibers, J. A.; Jorgensen, C. K.; Neilands, J. B.; Reinen, D.; Williams, R. J. P., Eds.; Springer: New York, 1974; pp 177−216. (63) Fripiat, J. I. In Soluble silicates; Falcone, J. S., Ed.; ACS Symposium Series 194, American Chemcial Society: Washington, DC, 1982; pp 165−184. (64) Argyris, D.; Cole, D. R.; Striolo, A. ACS Nano 2010, 4, 2035− 2042. (65) Fubini, B.; Areán, C. O. Chem. Soc. Rev. 1999, 28, 373−381. (66) Kubicki, J. D.; Blake, G. A.; Apitz, S. E. Am. Mineral. 1996, 81, 789−799. (67) Karlsson, M.; Craven, C.; Dove, P. M.; Casey, W. H. Aquat. Geochem. 2001, 7, 13−32. (68) Dove, P. M.; Han, N. In Perspectives on Inorganic, Organic, and Biological Crystal Growth: From Fundamentals to Applications; Skowronski, M.; DeYoreo, J. J.; Wang, C. A., Eds.; Springer-Verlag: New York, 2007. (69) Van der Voort, P.; Gillis-D’Hamers, I.; Vansant, E. F. J. Chem. Soc. Faraday Trans. 1990, 86, 3751−3755. (70) Stumm, W.; Furrer, G.; Kunz, B. Croat. Chem. Acta. 1983, 56, 593−611.

17491


A New Hypothesis for the Dissolution Mechanism of Silicates

Recommend Documents