- © 2003, The Clay Minerals Society
The atomic structure of dioctahedral 2:1 phyllosilicate edge surfaces was calculated using pseudopotential planewave density functional theory. Bulk structures of pyrophyllite and ferripyrophyllite were optimized using periodic boundary conditions, after which crystal chemical methods were used to obtain initial terminations for ideal (110)- and (010)-type edge surfaces. The edge surfaces were protonated using various schemes to neutralize the surface charge, and total minimized energies were compared to identify which schemes are the most energetically favorable. The calculations show that significant surface relaxation should occur on the (110)-type faces, as well as in response to different protonation schemes on both surface types. This result is consistent with atomic force microscopy observations of phyllosilicate dissolution behavior. Bond-valence methods incorporating bond lengths from calculated structures can be used to predict intrinsic acidity constants for surface functional groups on (110)- and (010)-type edge surfaces. However, the occurrence of surface relaxation poses problems for applying current bond-valence methods. An alternative method is proposed that considers bond relaxation, and accounts for the energetics of various protonation schemes on phyllosilicate edges.
- Ab Initio
- Clay Edge Surfaces
- Density Functional Theory
- Dissolution Kinetics
- Surface Structure
Phyllosilicate surface chemistry is an important control on ion mobility and pH in natural waters, permeability reduction in reservoir rocks, groundwater flow through engineered barriers, and clay rheology. In addition, phyllosilicate clays have a number of industrial applications due to their distinctive surface chemical properties. Therefore, a detailed atomic-scale understanding of phyllosilicate surface chemistry is desirable. While the basal surfaces of phyllosilicates have been the subject of numerous studies, the reactivity of edge surfaces is less well known because edges are difficult to isolate experimentally, and more difficult to model atomistically.
Edge and basal surfaces of phyllosilicates vary greatly in reactivity due to the extreme anisotropy of phyllosilicate crystal structures. For example, edge surfaces exhibit pH-dependent surface charging behavior, whereas siloxane-terminated basal surfaces do not (White and Zelazny, 1988; Bleam, 1993; Bleam et al., 1993). Furthermore, acid dissolution of phyllosilicates proceeds almost exclusively at edge surfaces (Kaviratna and Pinnavaia, 1994; Turpault and Trotignon, 1994; Rufe and Hochella, 1999; Bosbach et al., 2000; Bickmore et al., 2001). These differences in acid-base reactivity arise because edge surfaces are terminated by valence unsaturated oxygen atoms, whereas oxygen atoms on siloxane-terminated basal surfaces are valence saturated. Underbonded oxygen atoms electrostatically attract protons, and participate in pH-dependent charging behavior or bond dissociation (Xiao and Lasaga, 1994; Ganor et al., 1995; Bickmore et al., 2001). Solution-chemical investigations of dissolution and pH-dependent charging cannot be interpreted on a molecular level without assuming some sort of surface structural model.
The most useful models of phyllosilicate edge surfaces take into account site type, distribution and density. First attempts included schematic representations of surface site types (Schofield and Samson, 1953; Muljadi et al., 1966), but did not address site density or distribution. White and Zelazny (1988), Bleam (1993), Bleam et al.(1993), Brady et al.(1996) and Bickmore et al.(2001) used crystal chemical methods to predict site type, distribution and density. The precise termination of each edge surface is determined by cutting the crystal structure parallel to the plane of interest in such a way as to break the weakest bonds and preserve stoichiometry (Bleam et al., 1993). The same results are achieved by terminating each surface at the edges of periodic bond chains, as defined by Hartman-Perdock periodic bond chain (PBC) theory (Hartman and Perdock, 1955a, 1955b, 1955c; White and Zelazny, 1988; Bickmore et al., 2001). Crystal chemical methods predict that only two types of edge surfaces should exist, (010) and (110) (Figure 1⇓) (White and Zelazny, 1988; Bickmore et al., 2001). (This assignment of Miller indices assumes a 1M polytype. The (110) face is not symmetrically equivalent to the (11̅0) face in triclinic structures, but bond lengths and single-layer unit-cell parameters are nearly identical for 1M and 1Tc pyrophyllite structures (Bleam et al., 1993), so the two faces are treated as equivalent here.) Edge surfaces are stabilized in the presence of water by chemisorption of water species (Bleam et al., 1993). However, crystal chemical rules cannot always be used to distinguish between different surface protonation schemes, and there are differences in this regard among the edge surface models mentioned above (Figure 2⇓). Another complication is surface relaxation (i.e. shifts in electron density, and hence bond lengths, with respect to the bulk structure), which Bleam et al.(1993) did not consider. White and Zelazny (1988) proposed a protonation scheme for the (110) surface based on relaxation patterns observed for related minerals (Figure 2e⇓), but this model remains hypothetical.
Molecular modeling methods can be used to calculate simultaneously the lowest-energy configuration of attached protons and surface structural relaxation. Yet few molecular modeling calculations are available for edge structures of phyllosilicates, and none at the ab initio level of theory to our knowledge (Bleam, 1993; Bleam et al., 1993). In this paper we report the results of calculations using pseudopotential planewave density functional theory (PPW-DFT) for neutral, protonated, edge-surface structures of dioctahedral 2:1 phyllosilicates. The calculations test and improve the models discussed above by treating edge structures at a more fundamental level than has been attempted before. Among other things, we specifically investigated the structural relaxation and energetic stability of proton configurations postulated by White and Zelazny (1988). Implications of these calculations were explored with respect to the pH-dependent charging behavior and acid dissolution of dioctahedral phyllosilicates.
All structures and total energies were calculated using the PPW-DFT method as implemented in CASTEP (Payne et al., 1992). The generalized gradient approximation (GGA) and generalized gradient local spin approximation (GGS) were applied using the Perdew-Wang (Perdew and Wang, 1992) parameterization of the exchange-correlation functional, modified to work with planewave calculations (White and Bird, 1994). We used the CASTEP parameterization of ultrasoft pseudopotentials (Vanderbilt, 1990) without core corrections. Pseudopotentials were generated using the local density approximation (LDA/LSDA), meaning that the screening effect of the core electrons was modeled using LDA/ LSDA, whereas the screening effect of the valence electrons was modeled using GGA/GGS. This approach has been validated previously (Garcia et al., 1992) and successfully applied to structure optimizations of 2:1 phyllosilicates (Rosso et al., 2001).
Pyrophyllite edge structure models with Al3+ or Fe3+ occupying the octahedral sites were generated from optimized bulk structures ((Al,Fe)Si2O5(OH)). The atomic coordinates and cell parameters of bulk 1Tc pyrophyllite were optimized simultaneously without symmetry constraint (i.e. in the P1 space group). The GGA was used for Al2Si2O5 and GGS for Fe2Si2O5. For the latter, knowledge of the ordering of unpaired spins in ferripyrophyllite is incomplete, but anti-ferromagnetic coupling of high-spin Fe3+ within ideal dioctahedral layers has been observed (Coey, 1988). Thus, in our calculations, a net spin of zero was assigned to the unit-cell so as to impose an anti-ferromagnetic spin distribution for each 2:1 layer. One k-point was used (gamma point), which gave satisfactory results for closely related mineral structures (Bridgeman et al., 1996; Chatterjee et al., 2000; Rosso et al., 2001). Optimization was performed using a cutoff energy of 300 eV and a conjugate gradient electronic minimizer using a density mixing scheme (Kresse and Furthmuller, 1996). Finite basis set error estimations using methods described in Payne et al.(1992) were <0.1 eV/atom, indicating that the bulk structure calculations were converged with respect to basis set size. Although these computational conditions are chosen in part due to the need for computational efficiency for our large structures, similar conditions were found to provide accurate structures and energetics in a previous study on phyllosilicates (Rosso et al., 2001).
Edge structures were then excised from optimized bulk configurations. Because the PPW calculations in this study always have three-dimensional periodic boundary conditions, surfaces were constructed by building in a vacuum layer (Rosso, 2001). In our case, two mutually perpendicular vacuum layers were present for each edge model. Because we used only one 2:1 layer for the edge structure calculations, the first vacuum layer is oriented parallel to the basal plane. The other vacuum layer in each model separates the edge surfaces of interest. Thus, the edge models are best described as polymer-like, because they are structurally contiguous and are treated as infinitely periodic units along only one dimension. Edge terminations, generated by PBC theory, were identical to those described in White and Zelazny (1988) and other previous studies cited above. The vacuum layers separating the edge models between adjacent cells were always >10 Å; and models were designed with an inversion center to ensure that interactions across the vacuum layers were negligible. The infinitely repeating chains were chosen to run along the  and  vectors, generating edge surfaces of the (110) and (010) types, respectively (Figure 3⇓). Protons were added to edge sites to correspond with neutral protonation schemes proposed by White and Zelazny (1988) and Bleam (1993), and a new scheme that has never been suggested before (see Figure 2⇑), while keeping the basal surfaces unchanged (vacuum terminated). The coordinates of all atoms in each of the protonated edge models were energy optimized using the same treatment as described for the bulk structures above, except with fixed lattice parameters, where the unit-cell dimension along the edge vector was deter mined from the corresponding bulk optimization.
RESULTS AND DISCUSSION
Atomic positions and unit-cell parameters derived from the bulk structure optimizations are reported in Table 1⇓. Although these calculations were performed without symmetry constraints, the bulk structures converge to their proper space group symmetries (C1̅), as determined using a structural search tolerance of ~0.1 Å. The accuracy of the calculated structures can be evaluated further in two ways.
First, calculated unit-cell parameters can be compared with experimental values for 1Tc pyrophyllite (Lee and Guggenheim, 1981) and a 2M ferripyrophyllite (Chukhrov et al., 1979) (Table 2⇓). Unit-cell parameters for the calculated structures are given in the C1̅ space group for direct comparison to the 1Tc pyrophyllite structure. However, only the calculated a and b cell dimensions are directly comparable to those measured for ferripyrophyllite (Chukhrov et al., 1979), which was a 2M polytype. Calculated ferripyrophyllite values agree well with those measured. Calculated unit-cell parameters for pyrophyllite are close to measured values, except that the calculated c dimension is larger by 0.49 Å. Overestimated c-axis lengths may be due to the particular computational conditions applied. Our results for pyrophyllite and those from Teppen et al. (2002) were obtained using slightly different methods and produced slightly different c-axis dimensions. Their c-axis estimate is also overestimated, but is closer to the experimental value (larger by 0.16 Å). The overestimated lengths could also be related to the incomplete description of van der Waals attractive forces in DFT (Wu et al., 2001). Pyrophyllite layers are held together predominantly by van der Waals forces, and the c-axis dimension is sensitive to the description of interactions across the interlayer (Giese, 1975). However, the accuracy of the c-axis dimension becomes insignificant in the context of the edge surface calculations because the calculations are based on single isolated 2:1 layers making interlayer interactions irrelevant.
Second, we suggest that the set of predicted bond lengths in each structure should show systematic correspondence with the bond valence model of Brown and Altermatt (1985). The total Brown (1981) bond strength reaching each cation in the structure, computed using the relevant bond lengths predicted ab initio, should equal the cation’s formal charge. This evaluation method has at least one advantage – systematic deviations in the computed bond lengths can be identified for structures (such as that of pyrophyllite) comprising asymmetrically distorted polyhedra. For example, it is not clear that simply tabulating and averaging all computed Al–O bond lengths for comparison with a ‘grand universal’ mean observed Al–O distance would be sensitive to subtle computational errors. In fact, Brown’s (2002) ‘distortion theorem’ states, "For any ion, lengthening some of its bonds and shortening others, keeping the bond valence sum the same, will always increase the average bond length". Thus, comparing average Al–O bond length with some idealized value could lead to the conclusion that the computations overestimate bond lengths when, in fact, the opposite might be the case.
Brown and Altermatt (1985) fitted the following bond valence expression to a large number of crystal structures in the Inorganic Crystal Structure Database:
where s is bond valence in valence units (v.u.), r is bond length, r0 an arbitrary bond length fitted for each pair of atoms, and B is a fitted parameter usually equal to 0.37 Å. The summed strength of bonds reaching a cation should be equal to its formal charge. Based on our models, calculated total bond strengths reaching the Si atoms in pyrophyllite were 4.08–4.11 v.u. and 4.09 v.u. in ferripyrophyllite. Total bond strengths reaching Al atoms were 3.09–3.10 v.u. in pyrophyllite and total bond strengths reaching Fe atoms in ferripyrophyllite were 3.19–3.20 v.u. This means that the calculated bond lengths were, on average, ~0.01 Åtoo short for Si–O and Al–O bonds, and ~0.02 Åtoo short for Fe–O bonds. Estimated standard deviations for these bond lengths, reported by Brown and Altermatt (1985), range from 0.001 to 0.003 Å. Thus, the bond length underestimation is probably a real artifact of the PPW-DFT calculations. Nonetheless, the agreement between calculated and expected bond lengths is still quite good.
Using the PPW-DFT optimizations, the lowest energy configurations of protons were sought for each of the edge models, and the extent of surface relaxation was explored. For both (010)- and (110)-type faces, the protonation schemes suggested by White and Zelazny (1988) were clearly energetically favorable. Furthermore, the surface relaxation predicted for the (110)-type faces by White and Zelazny (1988) was qualitatively confirmed.
Edge structure calculations were based on models that were infinitely periodic in one direction. To avoid excessive computational costs, the models were designed to be no more than ~ 10 Åwide perpendicular to the edges. In cases where models are ‘thin’, it is important to verify that relaxation effects from the opposing surfaces do not penetrate the full width of the model. To test for such an effect, distortion of bond lengths in the relaxed structures was checked. Table 3⇓ shows the range of distortion for each type of bond in the interior polyhedra of each edge structure (Figure 2⇑), with respect to the same bonds in unrelaxed structures. These bond shifts are relatively small, so we conclude that the edge structure models are of sufficient width to adequately mimic both the bulk and edge structural environments.
Our surface structure calculations did not include an overlying bulk water phase, which might affect the relative stabilities of surface proton configurations. However, we will show that the valences of O–H bonds in the surface hydroxyl groups in our calculations were similar to those of O–H bonds in bulk water. That is, the O–H bonds were about 0.8 v.u., leaving 0.2 v.u. on the H to donate to H bonds. Therefore, the vacuum-terminated surface structures probably would not be affected much by the addition of bulk water.
The structures of (010)-type faces were calculated with the protonation schemes shown in Figure 2a,b⇑, and their total energies compared. The (010)-type faces are terminated by silanol and aluminol or ferrinol groups, and the proton configuration to neutralize this surface typically has been assumed to be like that depicted in Figure 2a⇑. The proton configuration depicted in Figure 2b⇑ should also be considered as a possibility, based on the results of this study which show that deprotonation of a silanol group occurs at a lower pH than deprotonation of one of the aluminols. However, when the starting structure in Figure 2b⇑ was optimized, it was so energetically unfavorable that a proton on one of the aluminol (pyrophyllite)/ferrinol (ferripyrophyllite) groups ‘hopped’ over to the adjacent deprotonated silanol group. This result is discussed further below.
The (110)-type edge structures with starting proton configurations as shown in Figure 2c,d⇑ were also optimized. The minimal edge unit of the vacuum-terminated (110)-type edge exposes four underbonded oxygen atoms for each 2:1 layer, consisting of two silanol groups, one Si–O–Al bridging oxygen, and one aluminol/ferrinol group. Each edge has a net charge of −4 that must be neutralized with protons. Adding one proton to each underbonded O atom at an edge is a commonly assumed scheme (Bleam, 1993) (Figure 2c⇑). Assuming that surface relaxation tends toward the neutralization of unsaturated valence on individual surface groups, White and Zelazny (1988) proposed an alternative protonation scheme that differs simply by doubly protonating the aluminol/ferrinol and leaving the bridging O atom unprotonated (Figure 2d,e⇑). Note that this difference requires only a small translocation of a proton from its optimized location on the bridging oxygen atom to the aluminol/ferrinol oxygen atom (a distance of ~2.2 Å). We optimized both configurations to assess their relative stabilities. The total energies of the two configurations indicate that White and Zelazny’s (1988) scheme is lower in energy by −0.44 eV for pyrophyllite and −0.25 eV for ferripyrophyllite per 2:1 edge. This indicates that the singly protonated aluminol/ferrinol groups have a higher proton affinity, i.e. they are more basic than the bare bridging O atom. Such a configuration appears to be partly stabilized by the participation of one of the Al–OH2/Fe–OH2 protons in hydrogen bonding to the bridging O atom. This result agrees in principle with the electrostatic calculations of Bleam et al.(1993) and the predictions of White and Zelazny (1988). The hydrogen-bonding interaction between the proton on the aluminol/ ferrinol groups and the bridging O atom on the (110)-type surface is weak (0.026–0.032 v.u.). Thus, although it is probable that this interaction would be screened by intervening water molecules, we would not expect the screening to cause a reversal in the relative stabilities of the two protonation schemes. Therefore, it is concluded that the protonation scheme in Figure 2c⇑ is less likely to be important than the scheme proposed by White and Zelazny (1988) (Figure 2d,e⇑).
White and Zelazny (1988) also predicted a certain amount of structural relaxation at the (110)-type surfaces, resulting in the configuration illustrated in Figure 2e⇑. Table 4⇓ shows the total strength of bonds reaching the oxygen atoms of the surface functional groups on the (010)- and (110)-type surfaces shown in Figures 2a and ⇑2d,e⇑, excluding O–H bonds. The values in Table 4⇓ are corrected for the tendency of our calculations to underpredict bond lengths by subtracting the average overestimation of bond valence for each type of bond in our calculated bulk structures. That is, 0.017 v.u. was subtracted from the computed valence of Al–O bonds, 0.033 v.u. from that of Fe–O bonds, and 0.025 v.u. from Si–O bonds. The O–H bonds on the surface functional groups were computed to have bond strengths of 0.75–0.80 v.u., in good agreement with the assumption of Hiemstra et al.(1996) that surface hydroxyl bond strengths are close to O–H bond strengths for bulk water (0.80 v.u.). White and Zelazny (1988) predicted that Si–O, Al–O and Fe–O bond lengths reaching the bridging oxygens should contract significantly, and the Al–O and Fe–O bond lengths in the aluminol and ferrinol groups should expand significantly. Therefore, due to this relaxation of the edge structure, the total strength of bonds reaching the bridging oxygens on (110)-type surfaces was predicted to be 2.00 v.u., and the strength of the Al–O and Fe–O bonds in the aluminol and ferrinol groups was predicted to be 0.00 v.u. On the other hand, if no structural relaxation occurs (Bleam et al., 1993), the total strength of bonds reaching the bridging oxygens should be ~1.50 v.u., and the strength of the Al–O and Fe–O bonds in the aluminol/ferrinol groups would be ~0.50 v.u.
The PPW-DFT calculations confirm the structural relaxation predicted by White and Zelazny (1988), although the degree of relaxation is not so extreme as they envisaged. The Si–O, Al–O and Fe–O bonds reaching the bridging oxygens do contract by 0.04–0.14 Å, and the Al–O and Fe–O bonds in the aluminol and ferrinol groups expand by 0.10–0.13 Å. The total corrected valence of bonds reaching the bridging oxygen is computed to be 1.861 v.u. for pyrophyllite and 1.757 v.u. for ferripyrophyllite. The Al–O bond strength in the aluminol (pyrophyllite) is 0.327 v.u., and the Fe–O bond strength on the ferrinol (ferripyrophyllite) is 0.316 v.u. Therefore, most of the structural relaxation predicted by White and Zelazny (1988) for the bonds reaching the bridging oxygens occurs, whereas predicted relaxation for the Al–O and Fe–O bonds in the aluminol/ferrinol groups is less.
The small relaxation of the aluminol/ferrinol groups on (110)-type surfaces takes on added meaning when compared to similar functional groups on the (010)-type surfaces. In the unrelaxed bulk structures, the Al–O bond on the (110)-type aluminol has a corrected bond valence of 0.431 v.u., and the Fe–O bond on the equivalent ferripyrophyllite site has a bond valence of 0.460 v.u. Unrelaxed Al–O and Fe–O bonds on the (010)-type surfaces have bond strengths of 0.569 and 0.551–0.569 v.u., respectively. Table 4⇑ shows that the Al–O and Fe–O bonds of singly-protonated (010)-type aluminol/ferrinol groups have bond strengths of 0.674 and 0.698 v.u., respectively. On the other hand, the Al–O and Fe–O bonds of doubly-protonated (010)-type aluminol/ferrinol groups have bond strengths of 0.290 and 0.332 v.u. These values are close to the bond strengths of doubly-protonated Al–O and Fe–O bonds in (110)-type aluminol/ferrinol groups, 0.327 and 0.316 v.u., respectively. Since the Al–O and Fe–O bond valences for aluminol/ferrinol groups are different between unrelaxed (010)- and (110)-type surfaces, but nearly identical on relaxed surfaces when the protonation is the same, it can be argued that the valence of Al–O and Fe–O bonds in the aluminol/ferrinol groups is controlled largely by the degree of protonation.
Finally, the corrected bond valences of the Si–O bonds of the silanol groups on the surfaces studied ranged from 0.99–1.05 v.u. These values bracket the ideal bond strength of 1.00 for Si–O bonds in silica tetrahedra.
Acidity of phyllosilicate edges
In the past, treatments of phyllosilicate surface-charging behavior have employed various surface complexation models (e.g. constant capacitance, triple-layer, etc.) to rationalize potentiometric titration data (e.g. Brady et al., 1996; Kraepiel et al., 1998). However, most previous modeling has assumed the existence of two or three generic functional groups with integer charges, and the possibility of fractional charges has been ignored. Furthermore, with few exceptions (Chang and Sposito, 1994, 1996; Brady et al., 1996) surface complexation models of clay minerals have not separated permanent structural charge located at the siloxane basal surfaces from pH-dependent charge on the edge surfaces (Kraepiel et al., 1998). Neither have crystallographically distinct edge surfaces been treated separately.
In this section, we attempt to provide the basis for a more realistic treatment of phyllosilicate acid-base reactivity by using calculated surface structures to estimate acidity constants for specific functional groups. However, the current model (Hiemstra et al., 1996) designed for predicting site acidity may be oversimplified with respect to surface relaxation. Therefore, site acidities will be predicted using the model of Hiemstra et al. (1996), and also by combining ab initio surface structure calculations with a revised bond-valence approach.
The current model (MUSIC).
Hiemstra and coworkers (e.g. Hiemstra et al., 1989, 1996, 1999) developed a multisite complexation (MUSIC) modeling approach capable of considering the contribution of specific sites on distinct crystalline surfaces. Site types are defined by crystallographic methods and intrinsic acidity constants are predicted based on the valence saturation of the oxo-and hydroxo-surface groups (Hiemstra et al., 1996). Valence saturation is calculated using the model of Brown and Altermatt (1985).
It should be possible to create MUSIC models for the acid-base behavior of euhedral clay samples using the information gained from the ab initio calculations reported here, including site coordination and bond lengths. Although we will not fit a MUSIC model to potentiometric titration data for phyllosilicates in this paper, the location, density, coordination, unsaturated valence and predicted intrinsic acidity constants of reactive sites on dioctahedral phyllosilicate edge surfaces are calculated.
Hiemstra et al.(1996) predicted site acidities using a modified form of the following equation:
where Ka is the intrinsic acidity constant, A equals 19.8, V is the valence of the surface oxygen (−2), and sj is the valence saturation of the surface oxygen, defined by:
Here sMe is the valence of the Me–O bond, and sH is the valence of the O–H bond to the surface oxygen if the base is a hydroxo- group. (This value was assumed by Hiemstra et al. (1996) to be 0.80 v.u., close to that for H–O bonds in bulk water.) Some of the surface oxygen valence is assumed to be neutralized by hydrogen bonds from surrounding water molecules. Oxygens have four empty orbitals with which to form Me–O, O–H, or hydrogen bonds, and n is the number free to form hydrogen bonds. However, for oxygens coordinated to one metal cation, it is usually assumed that one of the free orbitals is sterically hindered from hydrogen bonding at surfaces (m + n = 2 in equation 3), unless the surface structure is relatively open, in which case m + n = 3. One of the free orbitals may or may not be sterically hindered on surface oxygens coordinated by two metal cations. (For solution monomers, it is assumed that none of the orbitals is sterically hindered from bonding.)
Even though the bond valence method is an attempt to create physically realistic surface complexation models, it still requires several simplifying assumptions. For instance, Hiemstra et al.(1996) assumed the bond strength of H–O bonds on surface functional groups to be 0.80 v.u. Similarly, acidity constants for crystallographically distinct surface groups can sometimes be treated as equivalent, even though slight variations in bond valence exist. Also, surface reactions are written assuming Pauling bond valences for simplicity in charge bookkeeping (Hiemstra et al., 1996).
We calculated pKa values for reactions taking place on the various surface functional groups on (010)- and (110)-type surfaces using the MUSIC model (Table 5⇓). Bulk bond lengths were used to calculate unsaturated valence.
Generally, the calculated pKa values are similar to those assigned to analogous functional groups on simple oxides. For instance, we calculate pKa values of 7.6–8.7 for silanol groups, and Hiemstra et al.(1996) assigned an average pKa of 7.9 for these groups on silica surfaces. Similarly, calculated pKa values for the loss of a proton from a doubly-protonated silanol group ranged from −4.3 to −3.1, whereas Hiemstra and coworkers used a value of −4.0. Aluminol groups were calculated to have pKa values of 8.5–10.6, whereas Hiemstra et al.(1999) assigned an average pKa of 9.9 for aluminols on gibbsite surfaces. The calculated pKa values for ferrinol groups (8.5–10.7) are somewhat higher than the value of 7.7 assigned to these groups on goethite surfaces (Hiemstra et al., 1996). However, this is easily explained by the highly distorted nature of Fe3+ octahedra in the 2:1 structure.
Problems with the current model.
The above assignment of pKa values is somewhat unsatisfying for the following reasons. First, the bond-valence method of Hiemstra et al.(1996) explicitly assumes that no surface relaxation occurs, and therefore, bulk Me–O bond lengths can be used to calculate the unsaturated valence of surface oxygens. On the contrary, our calculations show that the Me–O bond lengths in the aluminol/ferrinol groups relax rather dramatically, the singly-protonated forms having shorter, and the doubly-protonated forms having longer, bond lengths than in the bulk structure. This difficulty is especially acute with respect to the estimation of pKa values for bridging oxygens on (110)-type faces. We have shown that the bonds to these oxygens relax considerably on neutral surfaces, so it is not apparent why bulk bond lengths should be adequate for predicting the acidity of these groups. Alternatively, we cannot simply substitute bond valences from calculated surface structures into equation 3 and expect to obtain realistic pKa values. For example, if the relaxed Al–O bond valence for singly-protonated aluminol groups on (010)-type surfaces were used (0.674 v.u.) in equation 3, the calculated pKa for the deprotonation of doubly-protonated aluminol groups would be 6.5, instead of 8.5–10.6.
Second, the treatment of solvation in this method is inconsistent. For instance, singly-coordinated surface oxygens are defined to have only two orbitals available for proton binding or H-bonding interactions (m + n = 2 in equation 3) because of steric constraints, but in order to predict properly the acidity constants for silanol groups, it has to be assumed that three orbitals are available (m + n = 3). Similarly, surface oxygens bound to two metal atoms can have either one or two orbitals available (m + n = 1 or 2) for proton binding or H-bonding interactions (Hiemstra et al., 1996). Thus, the number of free oxygen orbitals is effectively another fitting parameter.
Third, the idea that acidity constants would be correlated with total unsaturated valence on surface oxygens of the base is also surprising. Since O–H bonds involve only one O orbital, it seems more likely that acidity would correlate with the unsaturated valence available in a single orbital of the basic form.
Finally, application of the current bond-valence model predicts that silanol groups on (010)-type faces should be slightly more acidic than adjacent aluminol groups, but our calculations show that the protonation scheme in Figure 2b⇑ is quite unstable.
An older model revised.
It appears that bond relaxation may actually be a key to understanding surface acidity constants. Experimentation with different methods of predicting acidity constants for functional groups on our relaxed surfaces has led us to propose a revision of the bond-valence model first proposed by Brown (1981, 2002). This model is capable of incorporating surface relaxation, does not require inconsistent treatment of solvation, and properly divides bonding power among orbitals. While full development of this model awaits more computational and experimental data, early results for the calculated phyllosilicate surface structures are promising.
Brown’s (1981) model is based on a correlation between the Lewis base strength of the base in an acid-base reaction and measured pKa values of solution monomers. The correlation equation is
where Sb is the base strength, which is the total unsaturated valence divided by the number of free possible bonds. Consider, for example, the following reaction.
The valence of sulfur (+6) is divided among four S–O bonds, so that each S–O bond has a valence of 1.5. This leaves a total unsaturated valence of 2 v.u. Since oxygen is four-coordinated on average, we assume that there are four bonding ‘orbitals’ per oxygen, one of which is taken up by an S–O bond. Since there are three remaining possible external bonds per oxygen, and 4 oxygens, we divide the total unsaturated valence by 12 to obtain a Lewis base strength of 0.17 (2 v.u./12 possible external bonds). Substituting this value into equation 4, we obtain a pKa of 3.0. This is fairly close to the value reported by Perrin (1982) of 1.99.
It is somewhat more complicated to calculate the Lewis base strength of a protonated oxyanion. Consider the following reaction.
Assuming Si–O bonds of 1.00 v.u. and O–H bonds of 0.80 v.u., each of the oxygen atoms not connected to a hydrogen would have 1.00 v.u. of unsaturated valence. In the previous example, it was assumed that each oxygen atom could form four bonds, but Brown (2002) maintained that the protonated oxygen atoms only have enough unsaturated valence for one weak bond, so we assume that each of these is three coordinated, while the deprotonated oxygen is four coordinated. This leaves six possible external bonds among which to divide the unsaturated valence. Since there is 0.2 v.u. of unsaturated valence on each protonated oxygen, and 1.0 v.u. on the deprotonated oxygen, we divide the total 1.6 v.u. unsaturated valence between six possible external bonds to obtain a Lewis base strength of 0.27. Substituted into equation 4, this yields a pKa of 9.9, which exactly matches the experimentally determined value of 9.9 reported by Perrin (1982).
There are, however, potential problems with this method of calculating Lewis base strength. First, it is not clear why the unsaturated valence would be evenly distributed about the possible external bonds. Would not the base strength of possible bonds to the bare oxygen necessarily be greater than that of the possible bonds to protonated oxygens? Second, our calculations show that Me–O bond lengths tend to contract or expand in response to varying protonation of the oxygen. Therefore, it is not clear that all the Si–O bond valences in H3SiO−4 should be 1.00 v.u.
It is probable that similar base strengths can be estimated if bond relaxation is taken into account. For example, results from PPW-DFT optimizations of charged pyrophyllite and ferripyrophyllite edge surfaces (K.M. Rosso, unpublished results) produce Si–O bond valences of 1.1–1.3 v.u. for deprotonated silanol groups, whereas the calculations of neutral surface structures yield Si–O bond valences of ~1.00 for singly-protonated silanol groups. If we substitute 1.2 v.u. for the Si–O bond valence on the deprotonated oxygen, then divide the unsaturated valence of that specific oxygen between the three non-bonded orbitals, we obtain a Lewis base strength of 0.27 v.u. (0.8 v.u./3 possible bonds). This is exactly the same base strength as calculated above, and can be used directly in equation 4 to predict a pKa of 9.7. Sefcik and Goddard (2001) recently used DFT methods to calculate silicate anion geometries. Calculated Si–O bond lengths for the deprotonated oxygen in gas and aqueous phases yielded bond strengths of 1.17 and 1.18 v.u., respectively. Following the procedure outlined above, these values yield respective base strength values of 0.28 and 0.27, and pKa values of 10.4 and 9.9.
Will this method of Lewis base strength and pKa estimation hold true for oxyanions in general? Structural calculations on oxyanions with well-defined acidities need to be carried out to test this hypothesis. However, at this point the method can be used with our calculated phyllosilicate surfaces to obtain reasonable pKa values (Table 6⇓). Furthermore, the approach can explain aspects of our data that the model of Hiemstra et al.(1996) cannot.
Acid-base reactions on neutral surfaces are supplemented by those on charged surfaces containing fewer protons (K. Rosso, unpublished results) in Table 6⇑. Comparison of pKa values with those in Table 5⇑ shows that the new method predicts values near those predicted using the method of Hiemstra et al.(1996). Exceptions include silanol groups on (110)-type surfaces and (110)-type bridging oxygens. (We note that calculated pKa values for the deprotonation of doubly-protonated silanol groups varied widely in the negative range. This is probably due to the sensitivity of equation 4 in this range, but even still several of the calculated values agree well with those in Table 5⇑.)
The values of silanol pKas on (110)-type surfaces can be rationalized by observing that significant relaxation takes place on these faces, corresponding to shifts in bond valence. For example, since the Si atom bonded to the bridging oxygen shifts some of its bond valence to that oxygen, it cannot transfer as much valence to the adjacent silanol if that group is deprotonated. Therefore, the silanol becomes more basic.
If the relaxation of bonds between a metal cation and one functional group can affect the acidity of another functional group bonded to the same metal cation, this new method of pKa calculation may help explain our inability to optimize the surface structure in Figure 2b⇑. It is possible that doubly protonating both aluminols causes adjacent bonds in the structure to relax, with the net effect of making the silanol more basic. Thus, it would be energetically favorable for a proton to hop from one of the doubly-protonated aluminols to the silanol.
Implications of the proposed model.
More calculations and experiments must be performed to test and refine this model of acidity constant prediction, but if the model premise is correct then there are profound implications for multisite complexation modeling.
First, the relaxation effects described above imply that one cannot simply assign a single pKa to each type of surface group. Acidity of one group may be strongly affected by the protonation state of adjacent groups. In retrospect, this conclusion should hardly be surprising, and we are not the first to put forward this idea (Stumm, 1992). Existing surface complexation models are not equipped to deal with shifting acidity constants, but detailed ab initio studies may correct this deficiency by providing information on the proportion of various site types that can be expected to exhibit a given pKa.
Second, our proposed model for pKa prediction treats the solvation of surface functional groups consistently. In fact, no differences in solvation are postulated, whether between surface functional groups with different coordination, or between functional groups on surfaces and solution monomers. This conclusion contradicts recent work postulating that differences in solvation behavior contribute significantly to differences in the acidity of surface functional groups and analogous groups on solution monomers. Sverjensky and Sahai (1996) and Sahai (2002) showed that points of zero charge for simple oxide surfaces can be correlated to equations that take into account both average bond valence and average dielectric constant for the mineral structures. Although the correlation is completely empirical, the contribution of the dielectric constant was qualitatively related to Born solvation theory. We suggest that the correlation of the dielectric constant with surface acidity is due to the fact that the dielectric constant is a measure of the ‘polarizability’ of the crystal structure, and might thus be related to the ability of bonds at the solid surface to relax. It may be differences in bond relaxation on surfaces and analogous solution monomers that control differences in acidity.
Acid dissolution of phyllosilicates
It is generally agreed that acid dissolution of phyllosilicates proceeds almost exclusively at edge surfaces, and recent in situ atomic force microscopy (AFM) observations of the process confirm this (Rufe and Hochella, 1999; Bosbach et al., 2000; Bickmore et al., 2001). However, different mechanisms of proton attack have been proposed. After fitting a surface complexation model to kaolinite, Wieland and Stumm (1992) noted that the dissolution rate could be modeled as a first-order function of the concentration of protonated aluminol sites at the edge faces. It was proposed that the rate-controlling activated complex in the dissolution reaction has a stoichiometry of 1 H to 1 Al, and protonated aluminol groups are precursors to the activated complex. On the other hand, results of ab initio cluster calculations (Xiao and Lasaga, 1994) showed that the key step in the dissolution of aluminosilicate groups should be proton attack at the bridging oxygen. Protons sorbed at neighboring locations did not significantly affect the dissolution reaction. Ganor et al.(1995) proposed that the rate-controlling dissolution step is the breaking of bonds to Si–O–Al bridging oxygen atoms. They also showed that kaolinite dissolution rates were approximately proportional to the total surface proton concentration, explaining this by hypothesizing equilibrium between the various types of surface sites. Zysset and Schindler (1996) showed that the dissolution rate of montmorillonite is proportional to the total surface proton concentration, including protons sorbed in charged interlayers. As did Ganor et al.(1995), they proposed equilibrium among all surface protons, but noted that they could not use their data to distinguish whether the rate-determining step in acid dissolution is hydrolysis of Si–O–Al or Al–O–Al bonds at edge surfaces.
Bickmore et al.(2001) attempted to settle this controversy by observing phyllosilicate dissolution with in situ AFM. Lath-shaped particles of nontronite (a dioctahedral smectite with mainly Fe3+ in the octahedral sites) and hectorite (a trioctahedral smectite with mainly Mg2+ in the octahedral sites) were fixed on a polyelectrolyte-coated mica surface (Bickmore et al., 1999) and dissolved pH 2 HCl in the AFM fluid cell. The hectorite dissolution data were previously reported by Bosbach et al.(2000). It was observed that both minerals dissolved inward from the edges, but the reaction fronts on nontronite quickly became pinned at the pseudohexagonal face angles and then dissolution slowed to an imperceptible rate. In contrast, the hectorite reaction fronts showed no preferential orientation. Why were the pseudohexagonal nontronite faces so much more stable than the randomly oriented edges of nontronite and all the hectorite edges? Bickmore et al.(2001) used periodic bond chain theory to generate models of trioctahedral and dioctahedral edge faces. It was shown that if the surface relaxation proposed by White and Zelazny (1988) for the dioctahedral (110)-type surface occurs (Figure 2e⇑), all the stable faces have valence-saturated bridging oxygen atoms, whereas the bridging oxygens on unstable faces are valence unsaturated. No scheme of bond relaxation could be conceptualized at the unstable faces that would neutralize the entire unsaturated valence of bridging oxygens. Underbonded bridging oxygens would electrostatically attract protons, and if the protons bonded to bridging oxygens, they would become overbonded, destabilizing the structure. If the predicted surface relaxation does not occur (Figure 2d⇑), no such differences between stable and unstable faces could be detected. The acid dissolution behavior described seems to indicate that the rate-controlling step is the breaking of bonds to bridging oxygens at the edge surfaces. Furthermore, it seems to support White and Zelazny’s (1988) surface structure model, but the evidence is indirect and Bickmore et al.(2001) concluded that "molecular modeling calculations are needed to theoretically justify the predicted surface relaxation".
As discussed above, the PPW-DFT calculations reported here confirm the structural relaxation predicted by White and Zelazny (1988), but the calculated relaxation is not so severe that the formal charge on the (110)-type bridging oxygens is fully neutralized. In fact, the bridging oxygens are predicted to be slightly underbonded, but less so than bridging oxygens on trioctahedral (110)-type surfaces (Bickmore et al., 2001). This conclusion is consistent with a modified form of the explanation of Bickmore et al.(2001) for hectorite and nontronite dissolution behavior. That is, the rate-determining step of the dissolution process is the breaking of bonds to bridging oxygens at the edge surfaces, and hence the most stable edge faces during acid dissolution are those with little or no unsaturated valence on the bridging oxygens.
The structure of ideal pyrophyllite and ferripyrophyl-lite neutral edge surfaces has been calculated ab initio using PPW-DFT methods. The results of the calculations favor the surface protonation scheme predicted by White and Zelazny (1988) and Bleam et al.(1993), and predict that much of the surface relaxation envisaged by White and Zelazny (1988) for (110)-type surfaces should occur. Relaxed surface structures were also used to predict intrinsic acidity constants for acid-base reactions involving the surface functional groups present on (110)- and (010)-type edge surfaces. It was demonstrated that surface relaxation must be taken into account in the prediction of surface acidity constants, and a revision of the bond-valence method for pKa prediction was proposed for this purpose. Collectively, these results are consistent with the AFM observations of phyllosilicate dissolution behavior made by Bickmore et al.(2001).
This work was made possible by a user grant from the Environmental Molecular Sciences Laboratory, Pacific Northwest National Laboratory. Portions of this work were funded by a grant to KLN from the US Department of Energy Environmental Management Science Program (DE-FG07-99ER15009), and internal support funds to BRB from the BYU College of Physical and Mathematical Sciences. Also, work by RTC was made possible by funding provided by the US Department of Energy, Office of Basic Energy Sciences, Geosciences Research Program. Sandia is a multiprogram laboratory operated by Sandia Corporation, a Lockheed Martin company, for the US Department of Energy under contract DE-AC04-94AL85000. We also thank B. Teppen and J. Kubicki for helpful reviews of the initial manuscript. In particular, B. Teppen suggested we try the protonation scheme illustrated in Figure 2b⇑. I.D. Brown also offered valuable advice about calculating Lewis basicity.
- Received 8 August 2002.
- Revised 3 February 2003.