Quantitative Mechanistic Modeling of Silica Solubility and Precipitation

ARTICLE pubs.acs.org/JPCC

Quantitative Mechanistic Modeling of Silica Solubility and Precipitation during the Initial Period of Zeolite Synthesis Claire E. White,†,‡ John L. Provis,*,† Thomas Proffen,‡ and Jannie S. J. van Deventer† † ‡

Department of Chemical & Biomolecular Engineering, University of Melbourne, Victoria 3010, Australia Lujan Neutron Scattering Center, Los Alamos National Laboratory, New Mexico 87545, United States ABSTRACT: The mechanistic details of the structural changes occurring during the initial stages of silica polymerization to form gels or zeolites remain largely unknown due mainly to the complexity of sol
gel synthesis processes. Previous simulation studies have applied simple lattice models to qualitatively replicate spontaneous silica nanoparticle formation using fitted interaction energy parameters to replicate the behavior of real systems. This study moves for the first time to the use of quantum chemical-based interaction (dimerization) energies, determined through density functional theory computations, in a coarse-grained Monte Carlo simulation of the initial stages of gel/cluster formation in sodium silicate systems across a range of concentrations. The use of accurate dimerization energies as model inputs enables semiquantitatively accurate results to be obtained, as determined by comparisons with 29Si nuclear magnetic resonance data. The most concentrated system simulated (9.36 m silica concentration) undergoes Ostwald ripening, whereby a single large cluster forms, indicative of colloidal silica formation. Furthermore, as the extent of the reaction progresses, this initial nonequilibrated cluster (known as the primary amorphous phase in the zeolite synthesis literature) is progressively transformed to the secondary amorphous phase via structural rearrangements, without significant changes in the overall size of the cluster. These results are in good agreement with previous experimental studies on the nature of the amorphous phase(s) present during silicate gelation and prior to zeolite crystallization. Hence, this investigation demonstrates for the first time the successful application of multiscale simulation methodology to the coarse-grained Monte Carlo approach of a sol
gel process, revealing important quantitative mechanistic information regarding complex silicate reaction processes.

I. INTRODUCTION The use of synthetic zeolites in various industrial and technological situations has become more and more commonplace in many areas of industry due to their molecular sieving and catalytic properties.1
3 Zeolites are silica-based framework structures, and most are synthesized hydrothermally. Furthermore, during the hydrothermal synthesis of zeolites, the system undergoes a sol
gel transition, and it is during these initial stages of zeolite synthesis that aqueous silica chemistry becomes extremely relevant. Understanding silica chemistry (both structures and mechanisms), and the manner in which individual monomeric silica species impact the formation of zeolites, is of paramount importance in zeolite research and development. However, investigations in this area have been complicated by: (i) the number of individual processes occurring concurrently in alkaline silicate systems, and (ii) the current limitations in the implementation of full multiscale models. Hence, here we explore the mechanisms responsible for the early stages of silica zeolite formation using a molecular multiscale model of an alkaline silicate solution, to bridge the gap between the atomistic and mesoscale levels and to provide quantitatively accurate simulations of the behavior of silicate species in solution. Several previous theoretical and experimental investigations of aqueous silica chemistry under alkaline conditions2
14 have r 2011 American Chemical Society

enabled clarification of the individual processes, which are responsible, in large part, for the overall behavior of silica in zeolite formation. However, the atomistic theoretical models generally cannot replicate the mechanisms occurring at the mesoscale level in zeolites predominately due to the massive computational requirements associated with atomistic modeling across the necessary range of length scales. Population balance and dynamic Monte Carlo simulations have been developed based on reaction rate expressions14
16 but do not in general provide the necessary spatially resolved 3D insight into mechanisms and processes. Jorge et al.17,18 made significant advances in reducing much of this complexity by implementing a simplified lattice-based coarse-grained Monte Carlo (CGMC) model to simulate nanoparticle self-assembly, using an effective attraction force between lattice sites, electrostatic interactions by the use of orientationdependent terms, and size exclusion effects using second-neighbor interactions.17,18 In a later article from the same group, Jin et al. extended the model to simulate silica nanoparticle growth using a body-centered cubic (bcc) lattice to mimic the tetrahedral structure of silica.19 The results of the CGMC simulations agreed Received: January 20, 2011 Revised: April 1, 2011 Published: April 29, 2011 9879

dx.doi.org/10.1021/jp2006217 | J. Phys. Chem. C 2011, 115, 9879–9888

The Journal of Physical Chemistry C qualitatively with experimental results, including the spontaneous formation of silica nanoparticles under conditions where they are seen experimentally. However, their model required empirical calibration of parameter values, rather than taking absolute interaction energies as inputs. Thus, in the search for quantitatively accurate results from a CGMC model, it would appear logical that absolute interaction energy parameters should be utilized. Such parameters can be obtained via first principles calculations of monomer and dimer silica species, and the model should incorporate the full range of deprotonated silicate monomers and dimers that can form in alkaline solutions.20 This extension is not straightforward because the addition of extra species adds much complexity to the model in terms of the determination of bonding environments and in maintaining mass and charge balances. Nevertheless, the results of such a simulation will greatly benefit the current understanding of how and why silicate gels and zeolites form. Hence, in applying interaction energies obtained from density functional theory (DFT) computations, this investigation bridges the gap between atomistic and mesoscale modeling, making it a true multiscale model. It should be mentioned that this study focuses on modeling the solubility of silica in alkaline silicate systems, which is different from the previous investigations of Jorge et al. and Jin et al., where nanoparticle self-assembly was modeled.17
19 In this investigation, different silica concentrations have been simulated and their structures analyzed, revealing that the model is able to quantitatively replicate silica solubility in alkaline aqueous solutions. This article is organized as follows: in section 2 the model is reported, together with the methods used to describe the individual silica mechanistic events occurring during the initial stages of zeolite synthesis. In section 3, the simulation technique employed is outlined, and section 4 presents the results obtained, along with a discussion of the implications of these results and the scope for potential future work. Section 5 outlines the conclusions drawn from this investigation.

II. MODEL DESCRIPTION It is known that the atomic structures of silicate gels are complex, and therefore it is foreseen that a relatively large simulation system size is required to simulate the molecular structure. Furthermore, the number of molecular events required for the system to progress through silica precipitation and obtain a (metastable/kinetically limited) pseudoequilibrium state is predicted to be high, and therefore coarse-graining from an atomistic to a molecular (monomeric species) representation will be required for an equilibrium state to be approached at a reasonable computational cost.21,22 There is obviously a price associated with coarse-graining, as the results obtained are expected to be lower in precision due to the approximations required in the coarse-graining procedure. However, the accuracy of the results obtained depends largely on the quality of the input parameters and the use of values derived from detailed quantum chemical calculations20 may enable quantitatively accurate results to be generated. The model described here is a simple stochastic lattice model on a cubic lattice. The recent investigation of Jin et al.19 revealed that, for silicate systems, a body-centered cubic (bcc) lattice system is the most realistic lattice model for replicating the bonding environment present in pure silicate systems (without other species that are known to participate in higher bonding

ARTICLE

environments, such as calcium). However, earlier investigations by the same research group proved that the primitive cubic lattice does provide sensible results for silicate-based nanoparticle systems. Although it has been shown that the bcc lattice does better replicate silicate bonding in pure silicate systems, the cubic lattice has been employed in this current investigation so that minimal changes are required to extend to more complex systems, such as calcium silicate hydrates and geopolymers (alkali-activated aluminosilicates). Each lattice site in the primitive cubic lattice can be occupied by a single monomeric silicate species, and each silicate site can be bonded to at most four other neighboring sites, to represent the tetrahedral coordination of silicate gels within the constraints of a cubic lattice. As will be discussed in detail below, the water molecules and sodium cations present have been explicitly represented in the quantum chemical calculations in our previous investigation20 and therefore are not included explicitly in the Monte Carlo simulations; this assumption is an essential aspect of the coarse-graining procedure. However, knowledge of their existence is pivotal, particularly with regard to the calculation of initial species concentrations and the corresponding occupied lattice site percentages. There are six neighboring sites defined for each lattice site in the Monte Carlo lattice; sites must share a face in the cubic lattice to be considered neighbors; sharing an edge or a corner is not sufficient to enable bonding. This assumption is essential in ensuring the validity of the interaction energies obtained from DFT for structures with particular bond lengths and geometries; if sites separated by different distances were allowed to bond to each other in the model, the use of the same energy parameters would be invalid. If a lattice site is occupied by a silicate species, it can participate in bonding with its occupied neighbor sites. The model does not describe any change in the reactivity of a site depending on whether or not it is already bonded to other sites; the energy change associated with a site bonding to another site is calculated using the Gibbs free energy of dimerization of two monomeric species,20 irrespective of the existing bonding environment of the sites involved. This does not precisely represent the case for a real system because, as was discussed in detail by Sefcík and McCormick,13 the pKa of a silicate site is known to depend on its connectivity, and so it would be expected that the Gibbs energy of bond formation would also change. If this was a significant limitation in the simulation methodology, it should become apparent in the results obtained, particularly during comparison with previous literature related to coordination environments, and this will be discussed later in this investigation. In the following sections, the various chemical processes incorporated into the coarse-grained MC model are outlined. A. Polymerization of Species. The model presented here is designed to simulate the formation of the silicate gel from individual monomers via polymerization reactions. The most fundamental reaction in silica sol
gel chemistry is the dimerization of two neutral silicic acid species, Si(OH)4, which react to produce a dimer, (HO)3SiOSi(OH)3, and a water molecule. Hence, the cornerstone of the MC model presented in this investigation is the use of the energetics of this and other dimerization reactions involving silicate monomers. As was discussed in our previous investigation,20 there have been several previous experimental and theoretical investigations into silica dimerization reactions, reporting structural configurations and energetics. However, there has been a void in the literature regarding the 9880

dx.doi.org/10.1021/jp2006217 |J. Phys. Chem. C 2011, 115, 9879–9888

The Journal of Physical Chemistry C

ARTICLE

Table 1. Gibbs Free Energy of Dimerization Reactions (ΔGreaction) Occurring in Alkaline Silicate Solution at a pH of 11, as Determined in Our Previous Investigation20b monomer speciesa

M

M
3 Na 3 3H2O

M2
3 2Na 3 6H2O

M


1.8


9.3


5.3

M
3 Na 3 3H2O M2
3 2Na 3 6H2O


0.9

8.1 35.0

M = Si(OH)4, M
= SiO(OH)3
, M2
= SiO2(OH)22
. b All values are in units of kJ/mol. a

Table 2. Equilibrium Constants for Silicate Speciation, Obtained from the Literature As Noted species

pKa1

pKa2

pKa3

Q0 (monomer)

9.5a

12.6a

15.7a

1

Q (dimer) Q

2

Q3

b

13.25

c

d

9.85

11.2

b

16.65b

13.1

11.2c

a Sefčík and McCormick.13 b Sefčík and McCormick,13 where pKa1 is the average of the first and second pKa values of the dimer (pKa,d1 and pKa, 2 2 3 4 3 d ), pKa is the average of pKa,d and pKa,d in that article, and pKa is the average of pKa,d5 and pKa,d6. c Rimer et al.12 d Caullet and Guth,4 where pKa2 is the average of the values of pKa4, pKa5, and pKa6 for a trimer.

energetics and structural changes associated with dimerization reactions involving more highly deprotonated silica species, in particular doubly deprotonated silicate monomers. Our previous investigation reported the calculation of the energetics of all possible silica dimerization reactions, which can take place in a typical zeolite synthesis under alkaline conditions.20 These dimerization energies form the heart of this MC model, as they dictate polymerization. Table 1 summarizes the condensation reaction energies for dimerization of silicate monomers obtained in our previous work.20 The total energy of the system is calculated by summing the relevant dimerization energies for all occupied sites in the lattice according to the bonding environment present. B. Silica Deprotonation. The influence of pH on the system has been taken into account in this model first during density functional calculation of the dimerization reaction energies (by the use of the COnductor-like Screening MOdel (COSMO) with an appropriately selected dielectric constant),20 and second in the MC model calculation procedure by updating the distribution of deprotonated silica species regularly during the simulation runs to match the experimentally determined distribution reported in the literature. Simulations were carried out at a constant pH of 11. To update the distribution of silica species, the connectivity of each species is determined (Q0, Q1, Q2, Q3, or Q4), where a Qn species has n bonded neighbors, up to a maximum allowable number of four bonded neighbors. The pKa values used in the model depend on connectivity and are shown in Table 2. Every 30 iterations during the MC simulation, the distribution of deprotonated sites on silica species is determined. If the current distribution is different from the calculated equilibrium distribution (determined using the values in Table 2), then sites are selected at random and their protonation state altered to rectify the difference.

III. SIMULATION DETAILS The Monte Carlo simulations have been conducted in the canonical ensemble (NVT; constant number of occupied sites,

system volume, and temperature). This ensemble consists of a 125 000 site cubic lattice (503 sites) with periodic boundary conditions in all three Cartesian coordinates. Larger lattice sizes were investigated, and it was determined that this size (503 sites) provided a good balance between minimizing the artificial effects induced in small lattice simulations (refer to Figure 6 in ref 17) and obtaining convergence in a reasonable amount of time. It should be mentioned that the size of simulation used here is less than the size reported by Somvarsky and Dusek necessary to accurately replicate the infinite system (106
107 site simulations were the closest to replicating the real system), where kinetic Monte Carlo simulations were used to model network formation.25,23 However, the investigation by Herrmann et al.,24 which modeled irreversible gelation using a simple lattice approach ranging from 153 up to 603 sites, was capable of replicating the gelation phenomena using these smaller lattice sizes. Hence, the choice of 503 for the lattice size was deemed appropriate for the current investigation. The code was created in Microsoft Visual Studio 2008 Express Edition, and the lattice is depicted by a 3D char array. To discern between bonded and nonbonded occupied sites within the MC lattice, each occupied site also has a parameter defining whether it is bonded or nonbonded. The random number generator used is the minimal standard random number generator reported by Park and Miller,26 with the seed given by the Visual Studio 2008 Express Edition random number generator. Two functions are implemented to model polymerization and dissolution in the system, and are trialed alternately in the MC simulations. The first, the swap function, randomly selects one occupied site (any silicate species, bonded or nonbonded) and another site of any type (any site in the lattice, occupied or unoccupied). The sites are swapped, including their bonding information, if the swap lowers the total energy of the lattice. If it raises the total energy, the move may be accepted with probability X if it raises the energy, where X is calculated based on the Boltzmann factor associated with the configurational change (biased sampling),27 as described by eq 1 (where kB is the Boltzmann constant and T is the temperature), and ΔE is the change in energy that would occur if the move were carried out, as given in eq 2. This swap function enables the effect of a nonzero temperature to be modeled during the simulation, by allowing species to move. It should be noted that the total energy of the lattice is defined relative to the initial energy (the energy of the system at iteration zero), so there is not direct implementation of temperature effects via kinetic energy calculations. The other purpose of the swap function is to avoid the system falling into a local minimum and structurally not being able to get out of the configuration. Hence, this function helps enable the system to reach equilibrium. X ¼ e
ΔE=kB T

ð1Þ

ΔE ¼ Eafter
Ebefore

ð2Þ

The second function, the bond function, selects an occupied site at random (either bonded or nonbonded) and then checks to see if the site has a viable neighbor to bond with, as defined previously. If, by allowing the sites to bond, the total system energy is lowered, then the bonding is accepted. Otherwise the bonding move is once again accepted with probability X based on the Boltzmann factor.27 It should be noted that (i) bonding in the MC model is not directional, and therefore, if an occupied 9881

dx.doi.org/10.1021/jp2006217 |J. Phys. Chem. C 2011, 115, 9879–9888

The Journal of Physical Chemistry C

ARTICLE

Table 3. Compositions of Several of the Sodium Silicate Solutions Simulated, and the Corresponding Silica Concentrations in Parts Per Million (ppm) and Molality Units Lattice Site Percentage (%)

a

[SiO2] (ppm)

m (mol/kg)

water þ Naa

silica 0.15

1128

0.019

99.85

11 250

0.187

98.50

1.50

112 500

1.87

85.00

15.00

562 500

9.36

25.00

75.00

Water and sodium are represented by unoccupied sites in the lattice.

site is defined as bonded, then all the neighboring positions are checked, and any neighboring occupied site defined as bonded is bonded to this site, and (ii) the locations of negative charges do not change once bonding occurs (i.e., when two singly deprotonated sites bond, each site remains singly deprotonated). It is also worth noting the absence of an unbond function; such an event is implicitly included in the swap function, whereby a site can move to elsewhere in the lattice, and if there are no bonded neighbors at this new location the site reverts back to nonbonded. As mentioned previously, the two functions are called in an alternate fashion (i.e., swap, bond, swap, bond, and so on), so that 50% of trialed events are carried out using the swap function, the other 50% using the bond function. Through this stochastic process where silicate monomers can bond or swap sites, the system will move toward equilibrium. Oligomers (or polymerized structures) may form, which are characteristic of the formation of silicate gels, including zeolite precursor gels. Analysis of Model Outcomes. For each simulation run, the energy of the system is recorded at regular intervals (every 1 to 200 iterations, depending on the length of the simulation). One iteration is defined here as one swap or bond trial. Systems are allowed to equilibrate for up to several million iterations, until the energy profile reaches a constant value (to within a user defined convergence criterion, usually ∼ (5 hartree, where 1 Ha ≈ 2625.5 kJ/mol, for the more concentrated silicate systems). During simulations, the configuration of the box contents is also recorded periodically, usually every 10 000 to 100 000 iterations. These files (containing species information and x,y,z coordinates) are directly readable by visualization software; Materials Studio Visualizer (Accelrys, San Diego, USA) is used in this investigation, although various other visualization packages could also be used. More detailed analyses of cluster distributions during simulations have been carried out using purpose-designed algorithms, where a flood technique is implemented to search the occupancy and bond environment of all neighboring sites and thus identify all sites in a cluster. This program also records the local coordination environment (Qn), enabling direct comparison with 29Si nuclear magnetic resonance (NMR) experimental data. The MC simulation technique has been applied here to sodium silicate solutions of various concentrations. Fourteen different concentrations of silicate in sodium hydroxide solutions have been investigated; Table 3 describes four of these solutions, which span the range of concentrations investigated. The highest concentration (9.36 m silica, corresponding to approximately 562 500 ppm SiO2) has molar ratios Na2O/SiO2/H2O of 1:2:11, which is a highly concentrated and viscous solution, as previously reported by Provis et al.28 Hence, it is expected that the simulation for this concentration should produce

predominantly polymerized species; the lower SiO2 concentrations would be expected to have a correspondingly lower extent of polymerization. All simulations are initialized by filling the specified percentage of lattice sites with unbonded (monomeric) silicate units, meaning that there are no polymerized sites present at the start of the equilibration process; bonding is achieved only through the implementation of the bond event in the MC simulation.

IV. RESULTS AND DISCUSSION Energy Evolution and Cluster Size. Coarse-grained MC simulations of sodium silicate solutions have been carried out with various silica concentrations (including those shown in Table 3), allowing the systems to equilibrate for 400 000 to 6 000 000 iterations in each run. The energy evolution profiles of the simulations depicted in Table 3 are displayed in Figure 1. The total energy of the 0.019 m silica system (part a of Figure 1) remains relatively unchanged throughout the simulation, indicating that there are no significant structural rearrangements occurring in this system; the system was initialized with a small number of monomeric sites in the box and remains in this state throughout the simulation. Furthermore, the cluster analysis of the final lattice configuration after 400 000 iterations, as reported in Figure 2, shows that no clusters are formed. Once the silica content is increased to 0.187 m there are small structural changes occurring, as visible in the energy profile in part b of Figure 1, where there is a small decrease in energy during the first period of the simulation (over the first 10 000 iterations). This change in energy, and the associated small changes in structural bonding, can be confirmed by cluster analysis results, which show that 5.7% of the initial silicate monomers have formed clusters (Figure 2). Analysis of the cluster size distribution after 400 000 iterations shows that these clusters are predominantly dimers, with only a few trimers formed (Figure 3). Part c of Figure 1 displays the energy profile for the system with 1.87 m silica concentration. In this figure, it can be seen that there are obvious changes in the energy during the simulation, with an overall decrease in energy as the system evolves toward equilibrium. Cluster analysis of the final system configuration (after 400 000 iterations) reveals that 37.4% of the silicate monomers are situated in clusters (Figure 2). Figure 3 displays the cluster size distribution, which reveals that these clusters are predominately dimers and trimers, with only a small percentage of oligomers ranging in size between four and ten sites. Hence, in the 1.87 m system there has been a significant amount of silicate polymerization occurring to enable this system to reach equilibrium. When the silica concentration is increased to 9.36 m, cluster analysis reveals that practically all of the silica has polymerized (99.2%), as seen in Figure 2. Part d of Figure 1 displays the energy profile for this system, where 6 000 000 iterations are required for equilibrium to be reached. This system has an extremely highsilica concentration, which corresponds to an extremely viscous liquid, and therefore it would be expected that very large oligomers would form. As can be seen in Figure 3, this system predominantly consists of a single silicate particle containing roughly 100 000 silica sites, along with a small number of uncombined monomeric species. Confirming this observation, data for a silica concentration of 4.68 m are also shown in Figure 3; the final cluster in this case is approximately 38 000 units in size, and is in equilibrium with a solution containing 15% 9882

dx.doi.org/10.1021/jp2006217 |J. Phys. Chem. C 2011, 115, 9879–9888

The Journal of Physical Chemistry C

ARTICLE

Figure 1. Energy evolution profiles for sodium silicate solutions with silica concentrations (a) 0.019 m, (b) 0.187 m, (c) 1.87 m, and (d) 9.36 m. For each system, the zero-value for the relative energy is taken as the energy of the initial configuration (iteration 0). Energies are given in units of hartrees, where 1 Ha ≈ 2625.5 kJ/mol, to give a more convenient vertical axis scaling.

Figure 2. Percentage of silica species present as clusters at the end of each simulation.

of the total silica as monomers, and a small number (5% of the total silica) of small oligomeric species (less than 32 sites in size). Silicate Speciation (Qn) and Colloidal Silica. The results presented in this investigation are able to be directly compared with previous experimental and theoretical trends to provide validation that the DFT/CGMC simulation method does qualitatively replicate silica solubility. The investigation of Provis et al.28 studied speciation in sodium silicate solutions of high alkalinity using 29Si NMR spectroscopy and the Pitzer activity coefficient model,29 together with extensive comparisons with

Figure 3. Cluster size distributions for a selection of the silica systems investigated after reaching equilibrium.

previous literature. Figure 4 displays the results from ref 28, which are compared directly with the results obtained in the current study. Provis et al.28 excluded Q4 sites from their model and 29Si NMR investigation, as the concentrations of these sites have been shown to be strongly and sometimes unpredictably dependent on the silica source,30,31 and their quantification by NMR is problematic due to the much longer relaxation delays required to fully sample these site populations compared to the lower-coordination sites.32,33 However, in this investigation there is a substantial fraction of Q4 sites present, particularly at high-silica concentrations, as seen in Figure 5. 9883

dx.doi.org/10.1021/jp2006217 |J. Phys. Chem. C 2011, 115, 9879–9888

The Journal of Physical Chemistry C

Figure 4. Silica speciation for all the systems investigated after reaching equilibrium (lines), compared with experimental results from Provis et al. (markers).28 Solid line and square markers represent Q0 species, dashed line and triangle markers represent the sum of Q1
Q3 species.

Figure 5. Percentage of silica sites existing as Q4 sites for the range of concentrations investigated after each simulation reached equilibrium.

Therefore, even though previous presentations of experimental results report only Q0 to Q3 speciation, the current investigation also reports the amount of Q4 sites in all systems. Figure 4 displays the percentage of monomeric species (Q0 sites) present across the range of silica concentrations modeled using the DFT/ CGMC methodology, together with the total percentage of oligomeric species excluding Q4 (i.e., the sum of all Q1 to Q3 environments). These results are shown in Figure 4 in conjunction with the experimental results of Provis et al.28 There is generally good agreement between the simulation and experimental results, especially at high-silica concentrations. In fact, the agreement is surprisingly good considering the number of approximations inherent in the coarse-graining procedure, and the fact that the model contains no fitted parameters. There are evident deviations from the experimental data in Figure 4 between silica concentrations of 2
6 m, which could be artifacts of the superficial bonding environment implemented in this study (primitive cubic lattice), preventing the occurrence of important silicate species, such as cyclic trimers. Instead, this simulation methodology allows for the formation of 4-rings. Furthermore, the absence of interaction parameters beyond the nearest neighbor coordination shell, together with the assumption that the reactivity of a site is not affected by existing bonds, does simplify the simulation methodology and therefore influence the results obtained. Nevertheless, even with these evident

ARTICLE

limitations, the results highlight the value of the DFT/CGMC methodology, and show that the use of genuine interaction energy parameters in a simulation such as this, although computationally demanding as a means of parameter determination, provide the opportunity to quantitatively replicate the chemistry controlling monomeric and oligomeric speciation in these solutions. The change in slope of the simulated Q0 and Q1
3 curves in Figure 4 coincides with a transition from the lattice consisting of small and medium size clusters (less than ∼300 sites in size) to one large cluster (>30 000 sites in size) together with smaller clusters of less than 50 sites in size. This transition happens between 3.28 and 4.68 m. Hence, this DFT/CGMC simulation methodology predicts that the gelation of sodium-silicate solutions occurs between these two silica concentrations. To further understand the outcomes obtained through the DFT/CGMC simulation methodology, direct comparison with experimental data of the percentages of Q0 through Q3 can be carried out. For a sodium silicate solution with molar ratios SiO2/Na2O = 0.5 and H2O/Na2O = 11 (i.e., the same SiO2/Na2O ratio as in the simulations, and a silicate concentration of 2.52 m), which is between the 2.34 and 2.81 m solutions studied in this investigation, the percentage of monomers was experimentally determined to be 62%, with 23% Q1, 15% Q2, and 0% Q3 sites.28 For the DFT/CGMC simulations, linear interpolation to the same concentration gives 50% monomers, 24% Q1 and 17% Q2, with the remainder in higher-coordination . Hence, there is generally excellent agreement with the results of Provis et al.,28 especially for the Q1 and Q2 sites. There are some discrepancies visible in Figure 4 between the experimental data and the results from this study, especially at the midrange silica concentrations (2
6 m), but the agreement with experiment is in general very good indeed. There are two possible explanations for the presence of a significant fraction of Q4 sites in the model results presented here. The first explanation concerns colloidal silica, and whether this DFT/CGMC methodology is in fact replicating the formation of colloidal and nanocrystalline silica precipitates, which are typically formed during the early stages of zeolite synthesis,34
36 as well as during extended room temperature storage of sodium silicate solutions.37,38 Furthermore, as mentioned by Cundy and Cox,3 the primary amorphous phase is colloidal, and therefore invisible to the naked eye. Hence, it is highly probable that the concentrated silica systems studied in this investigation are in fact replicating the formation of these colloidal precipitates. This phenomenon can be directly seen in Figure 3, where between the systems with silica concentrations of 1.87 and 4.68 m, a large silica precipitate begins to form in preference to small oligmeric species. The second explanation for the significant amount of Q4 species present at high-silica concentration is related to the tendency for the DFT/CGMC methodology to overestimate the extent of bonding due to the absence of second-nearest neighbor interactions. In previous investigations using CGMC to analyze silica polymerization, repulsive second-nearest neighbor interactions were included among the empirically determined interaction parameters.17
19 Obtaining accurate absolute second-nearest neighbor interaction parameters using density functional modeling is by no means a trivial task, and has not been carried out here. Also worth noting is the additional simplification in the bond event implemented here, where the reactivity of a site in the model does not change as a result of pre-existing bonds. Nevertheless, if these aspects posed a significant 9884

dx.doi.org/10.1021/jp2006217 |J. Phys. Chem. C 2011, 115, 9879–9888

The Journal of Physical Chemistry C

ARTICLE

Figure 6. 2D slices of the MC model 3D lattice for the simulation at 1.87 m silica concentration, showing the clusters present (nonbonded monomers have been omitted for visual clarity; sites that appear isolated in the single slices are bonded to sites in the slices above or below). Yellow spheres represent neutral bonded silicate sites, and green spheres are singly bonded deprotonated silicate sites.

limitation in the DFT/CGMC methodology then there would not be such good agreement with the experimental monomeric speciation results as is evident from Figure 4. It should also be mentioned that quantitative analysis of silica speciation, especially Q4 sites, using 29Si NMR is not straightforward. Several previous investigations have not included the Q4 sites during quantification due to the long relaxation times required and the relatively ill-defined Q4 peak.39,40 Furthermore, attempts to quantify the Q4 sites have been hindered by the overlap of this peak with instrumental contributions of amorphous silica.41 The investigation of Bahlmann et al. reported a sizable percentage of Q4 (approximately 20%) in a sodium

silicate system with 25 wt % SiO2,30 however, those authors did not study a range of solutions with varying silica concentrations, which limits comparison with the current investigation. In light of the points raised in the previous paragraphs, the sizable fraction of Q4 sites at high-silica concentrations in the model presented here probably overestimates the percentage of these sites. However, as mentioned before, existence of these sites in high-silica solutions has been seen in the past,30 and therefore the DFT/CGMC methodology is again able to replicate (although probably strictly semiquantitatively in this aspect) the structural mechanisms occurring during the initial stages of zeolite formation. 9885

dx.doi.org/10.1021/jp2006217 |J. Phys. Chem. C 2011, 115, 9879–9888

The Journal of Physical Chemistry C

Figure 7. Monomer and cluster size distributions in the sodium silicate systems at various stages during the simulation, for (a) silica concentration of 1.87 m and (b) 9.36 m. Note that both axes for (b) are logarithmic.

Cluster Size Distribution and Ostwald Ripening. To further investigate the types of small clusters formed during the MC simulation (as seen in Figure 3, particularly for the 1.87 m system), 2D slices of the lattice obtained after 400 000 iterations at a silica concentration of 1.87 m have been obtained by sectioning the lattice parallel to the Y
Z plane, and five of these slices are displayed in Figure 6. The slices in this Figure show that there exist both clumps (cyclic species) and chains. The types of silica bonding present in the clusters vary, and only bonds involving lower charge states are in evidence (neutral
neutral, neutral
singly deprotonated, singly
singly deprotonated), since doubly deprotonated species are only present in clusters at very low concentrations due to the relatively high second deprotonation constant of silica (Table 2). As was shown in Figure 3, there are significant differences in cluster size distributions when the silica concentration is altered, especially between silica concentrations of 1.87 and 4.68 m. Further mechanistic information regarding the initial stages of zeolite formation can be gained by studying the change in cluster size distribution as a function of the number of MC events (as a proxy for time) as the various systems move toward equilibrium. Figure 7 displays the changes in cluster size distribution for the systems with silica concentrations of 1.87 m (part a of Figure 7) and 9.36 m (part b of Figure 7). There are evident differences between the two systems, especially with regard to the sizes of the clusters formed, and also in the variations in cluster size distributions as the systems approach equilibrium, with the lowersilica system dominated by small species, and the higher-silica system containing a single large cluster in equilibrium with a lower concentration of monomeric species, as noted above. For the lower-silica system (1.87 m, part a of Figure 7), there are apparent changes in cluster size taking place during the first 40 000 iterations; however, all of these changes involve small

ARTICLE

clusters (20 monomers or less), and equilibrium is reached when 63% of the SiO2 is present as monomers. In contrast, for the highsilica system (9.36 m, part b of Figure 7) practically all of the monomers are consumed once equilibrium is reached, with less than 1% of the SiO2 remaining dissolved as monomers in solution by 6 000 000 iterations. Furthermore, there are obvious changes happening in cluster distribution between 0 and 400 000 iterations, as is evident in part b of Figure 7, where after 40 000 iterations there are ∼67% monomers, 13% of sites made up of small and medium clusters (up to 450 occupied sites in size), and one large cluster consisting of ∼19 000 sites (20% of occupied sites). By 400 000 iterations, all of the small and medium size clusters have disappeared in preference to the growth of the largest particle. Although the simulation does not inherently provide a time scale for the reaction, it is interesting to note that the dissolution of solid silicates into alkaline solution is believed to take place mainly via the release of silicate monomers,42,43 and that the hydrolysis of an alkoxide precursor in a sol
gel silicate process will also provide a source of monomeric silicate units which can then participate in reactions. This means that the initial state of the MC simulations here (taken as a box containing only monomeric silicate units) may in some way represent the early stages of these processes and that the evolution toward equilibrium as measured in terms of MC events may to some extent parallel the time evolution of such a system. Thus, the observed behavior as the extent of reaction progresses for the high-silica system (9.36 m) can be identified as displaying an Ostwald ripening-type process. In this system, as the simulation progresses, the smaller clusters disappear, and there is a preference for growth of the larger precipitates at their expense. This is seen in part b of Figure 7 by the disappearance of small and medium clusters (less than ∼450 sites), together with the growth of the largest particle. Ostwald ripening is a common mechanistic observation in many aqueous silica-based and aluminosilicate-based systems, including silica nanoparticle evolution15,17 and zeolite formation.44 The results presented in this investigation also show a clear tendency for Ostwald ripening to occur in the high-silica concentration system. The phenomenon of Ostwald ripening includes the dissolution of intermediate clusters, and therefore requires the simulation methodology to be able to replicate such events. Hence, the methodology must not only be able to model dissolution but also diffusion. The methodology used in this current investigation artificially models the start and end points of diffusion through implementation of the swap function. This simple implementation of diffusion enables the phenomenon of Ostwald ripening to evolve, as has been previously proposed to occur in zeolite synthesis systems. Primary and Secondary Amorphous Phases. Part b of Figure 7 displays another mechanistic detail regarding the structural changes occurring during zeolite formation. Even though it takes 6 000 000 iterations for the 9.36 m system to reach equilibrium, the single large cluster (>90 000 sites) forms before 400 000 iterations. This indicates that the initial large cluster undergoes extensive structural rearrangement (without significant changes in cluster size) before the system reaches equilibrium, which agrees with the energy profile in part d of Figure 1, where by 400 000 iterations the system is approximately 75% of the way to reaching equilibrium. These results agree with the discussion in the review of Cundy and Cox,3 where it is highlighted that numerous previous investigations have reported 9886

dx.doi.org/10.1021/jp2006217 |J. Phys. Chem. C 2011, 115, 9879–9888

The Journal of Physical Chemistry C the presence of two amorphous phases in zeolite synthesis that occur consecutively, the primary amorphous phase and the secondary amorphous phase. These phases exist prior to crystal nucleation, and therefore are difficult to characterize experimentally due to (i) their amorphous nature and (ii) their temporary existence. Hence, this investigation provides direct evidence of these two intermediate products. Analysis of the silicate bonding environment for the 400 000 and 6 000 000 iteration structures reveals that 3% of the Q3 sites at 400 000 iterations transition to Q4 sites by 6 000 000 iterations, which is indicative of structure densification during the transition from the primary amorphous phase to the secondary amorphous phase. All other Qn site populations remain relatively unchanged. It should be noted that this investigation does not replicate nucleation of the zeolitic crystals, because our DFT/CGMC simulation methodology (i) does not model an increase in temperature after the formation of the secondary amorphous phase, and (ii) does not include intermolecular interactions. As mentioned by Cundy and Cox,3 in the final stage of zeolite synthesis (without the inclusion of templating), an increase in temperature is usually required for a prolonged period in order for the intermediate phase (secondary amorphous phase) to transition into the crystalline zeolite product. It may be possible to incorporate such functions in future simulations and therefore be able to mechanistically model the entire reaction process (from precursor hydrolysis or dissolution to crystal nucleation) for zeolite formation. Implications and Outlook. As was mentioned in the Model Description section, one of the main differences between this coarse-grained MC methodology and the model presented by Jorge et al.17,18 is that the interaction energy values used in this investigation are taken directly from quantum chemical calculations, whereas in the case of Jorge et al.17,18 the values were relative and required calibration. Hence, the DFT/CGMC methodology used in this investigation is truly a multiscale simulation technique. As is evident from the results in this study, the agreement with previous literature is quantitative, not just qualitative, and therefore this investigation provides significant further advances toward accurately replicating the mechanistic changes occurring during complex sol
gel reactions such as the initial stages of zeolite synthesis. The simulations presented in this investigation form the foundation for a model, which can be developed further to quantitatively describe zeolite and aluminosilicate gel formation. Much can be expanded, including description of the impact of silicate/aluminosilicate precursor dissolution at the beginning of the process, and replicating the effects of elevated temperature and subsequent crystal nucleation. Furthermore, the effect of varying the charge-balancing cation (including organic structuredirecting cations) still remains unexplored using the multiscale DFT/CGMC methodology. It should also be noted that application of this methodology is by no means limited to investigating the initial stages of zeolite synthesis; many related silicate and aluminosilicate processes could potentially benefit from such methodology, including and not limited to mineral dissolution, geopolymerization, calcium silicate hydrate structure formation, and biomineralization.

V. CONCLUSIONS Coarse-grained Monte Carlo simulations, based on nearestneighbor interaction energies derived directly from DFT

ARTICLE

computation, have been used to elucidate important details of the initial stages of silica polymerization. Sodium silicate solutions of varying silica concentrations were investigated, with silica concentrations ranging from 0.019 to 9.36 m. Analysis of the intermediate and equilibrium configurations revealed that at high-silica concentrations (above 4 m) the systems tend to form a single large cluster in preference to the many small oligomeric clusters observed at lower silica concentrations. The distribution of connectivities among the dissolved silicate species replicates the connectivity distribution observed by 29Si NMR, showing that the use of fundamentally derived interaction energies in a Monte Carlo model can provide quantitatively accurate descriptions of complex chemistry. Analysis of the structural changes occurring during the simulation of the highest-silica system (9.36 m silica concentration) reveals that the initial solid phase that forms is not yet equilibrated. This phase, known as the primary amorphous phase, is colloidal in nature, and undergoes structural rearrangements and densification to form the secondary amorphous phase, which is the phase present prior to zeolite crystallization. Structural rearrangement of the initial colloidal phase has been previously reported in the literature and is also able to be reproduced by the model. Hence, this multiscale DFT/CGMC methodology is, uniquely among the models reported to date, capable of quantitatively replicating chemical and microstructural aspects of the initial stages of silicate polymerization prior to crystallization from a fundamental standpoint.

’ AUTHOR INFORMATION Corresponding Author

*E-mail: [email protected], phone: þ61 3 8344 8755, fax: þ61 3 8344 4153.

’ ACKNOWLEDGMENT The authors would like to thank Dr. Daniel Riley, University of Melbourne/Australian Nuclear Science and Technology Organisation, for computational hardware support and insightful discussions. This work was funded in part by the Australian Research Council (ARC) (including some funding via the Particulate Fluids Processing Centre, a Special Research Centre of the ARC), and in part by a studentship paid to Claire White by the Centre for Sustainable Resource Processing via the Geopolymer Alliance. ’ REFERENCES (1) Auerbach, S. M.; Carrado, K. A.; Dutta, P. K. Handbook of Zeolite Science and Technology; Dekker Inc.: New York, 2003. (2) Cundy, C. S.; Cox, P. A. Chem. Rev. 2003, 103, 663–702. (3) Cundy, C. S.; Cox, P. A. Microporous Mesoporous Mater. 2005, 82, 1–78. (4) Caullet, P.; Guth, J. L. In Zeolite Synthesis; Occelli, M. L., Robson, H. E., Eds.; American Chemical Society: Washington DC, 1989; p 83
97. (5) Gomes, J. R. B.; Cordeiro, M. N. D. S.; Jorge, M. Geochim. Cosmochim. Acta 2008, 72, 4421–4439. (6) Kinrade, S. D.; Pole, D. L. Inorg. Chem. 1992, 31, 4558–4563. (7) Knight, C. T. G.; Balec, R. J.; Kinrade, S. D. Angew. Chem., Int. Ed. 2007, 46, 8148–8152. (8) Mora-Fonz, M. J.; Catlow, C. R. A.; Lewis, D. W. Angew. Chem., Int. Ed. 2005, 44, 3082–3086. 9887

dx.doi.org/10.1021/jp2006217 |J. Phys. Chem. C 2011, 115, 9879–9888

The Journal of Physical Chemistry C

ARTICLE

(9) Mora-Fonz, M. J.; Catlow, C. R. A.; Lewis, D. W. J. Phys. Chem. C 2007, 111, 18155–18158. (10) Nangia, S.; Garrison, B. J. J. Phys. Chem. A 2008, 112, 2027–2033. (11) Pereira, J. C. G.; Catlow, C. R. A.; Price, G. D. J. Phys. Chem. A 1999, 103, 3252–3267. (12) Rimer, J. D.; Lobo, R. F.; Vlachos, D. G. Langmuir 2005, 21, 8960–8971. (13) Sefcík, J.; McCormick, A. V. AIChE J. 1997, 43, 2773–2784. (14) Sefcík, J.; Rankin, S. E. J. Phys. Chem. B 2003, 107, 52–60. (15) Provis, J. L.; Vlachos, D. G. J. Phys. Chem. B 2006, 110, 3098–3108. (16) Rankin, S. E.; Kasehagen, L. J.; McCormick, A. V.; Macosko, C. W. Macromolecules 2000, 33, 7639–7648. (17) Jorge, M.; Auerbach, S. M.; Monson, P. A. J. Am. Chem. Soc. 2005, 127, 14388–14400. (18) Jorge, M.; Auerbach, S. M.; Monson, P. A. Mol. Phys. 2006, 104, 3513–3522. (19) Jin, L.; Auerbach, S. M.; Monson, P. A. J. Phys. Chem. C 2010, 114, 14393–14401. (20) White, C. E.; Provis, J. L.; Kearley, G. J.; Riley, D. P.; van Deventer, J. S. J. Dalton Trans. 2011, 40, 1348–1355. (21) Chatterjee, A.; Vlachos, D. G.; Katsoulakis, M. A. J. Chem. Phys. 2004, 121, 11420–11431. (22) Katsoulakis, M. A.; Vlachos, D. G. J. Chem. Phys. 2003, 119, 9412–9427. (23) Somvarsky, J.; Dusek, K. Polym. Bull. 1984, 33, 369–376. (24) Somvarsky, J.; Dusek, K. Polym. Bull. 1984, 33, 377–384. (25) Herrmann, H. J.; Stauffer, D.; Landau, D. P. J. Phys. A: Math. Gen. 1983, 16, 1221–1239. (26) Park, S. K.; Miller, K. W. Commun. ACM 1988, 31, 1192–1201. (27) Frenkel, D; Smit, B. Understanding Molecular Simulation: From Algorithms to Applications; Academic Press: San Diego, 1996. (28) Provis, J. L.; Duxson, P.; Lukey, G. C.; Separovic, F.; Kriven, W. M.; van Deventer, J. S. J. Ind. Eng. Chem. Res. 2005, 44, 8899–8908. (29) Pitzer, K. S. J. Phys. Chem. 1973, 77, 268–277. (30) Bahlmann, E. K. F.; Harris, R. K.; Metcalfe, K.; Rockliffe, J. W.; Smith, E. G. J. Chem. Soc. Faraday Trans. 1997, 93, 93–98. (31) Walker, A. J.; Whitehead, N. J. Appl. Chem. 1966, 16, 230–238. (32) Haouas, M.; Petry, D. P.; Anderson, M. W.; Taulelle, F. J. Phys. Chem. C 2009, 113, 10838–10841. (33) Provis, J. L.; Gehman, J. D.; White, C. E.; Vlachos, D. G. J. Phys. Chem. C 2008, 112, 14769–14775. (34) Mintova, S.; Olson, N. H.; Valtchev, V.; Bein, T. Science 1999, 283, 958–960. (35) Tosheva, L.; Valtchev, V. P. Chem. Mater. 2005, 17, 2494–2513. (36) de Moor, P. E. A.; Beelen, T. P. M.; van Santen, R. A. Micropor. Mater. 1997, 9, 117–130. (37) Vail, J. G. Soluble Silicates - Their Properties and Uses; Reinhold: New York, 1952. (38) Provis, J. L. In Geopolymers: Structures, Processing, Properties and Industrial Applications; Provis, J. L., van Deventer, J. S. J., Eds.; Woodhead: Cambridge, UK, 2009; p 50
71. (39) Kinrade, S. D.; Swaddle, T. W. Inorg. Chem. 1988, 27, 4253–4259. (40) Provis, J. L.; van Deventer, J. S. J. J. Mater. Sci. 2007, 42, 2974–2981. (41) Harris, R. K.; Bahlmann, E. K. F.; Metcalfe, K.; Smith, E. G. Magn. Reson. Chem. 1993, 31, 743–747.  izmek, A.; Subotic, B. J. Chem. Soc. Faraday Trans. (42) Antonic, T.; C 1994, 90, 1973–1977. (43) Meza, L. I.; Anderson, M. W.; Slater, B.; Agger, J. R. Phys. Chem. Chem. Phys. 2008, 10, 5066–5076. (44) Navrotsky, A. Proc. Natl. Acad. Sci. U.S.A. 2004, 101, 12096– 12101.

9888

dx.doi.org/10.1021/jp2006217 |J. Phys. Chem. C 2011, 115, 9879–9888