Article pubs.acs.org/JPCC
Intrinsic Acidity of Surface Sites in Calcium Silicate Hydrates and Its Implication to Their Electrokinetic Properties Sergey V. Churakov,*,† Christophe Labbez,§ Luis Pegado,† and Marialore Sulpizi‡ †
Laboratory for Waste Management, Paul Scherrer Institute, CH-5232 Villigen-PSI, Switzerland Johannes Gutenberg University Mainz, Staudinger Weg 7, 55099 Mainz, Germany § Laboratoire Interdisciplinaire Carnot de Bourgogne, UMR 6303 CNRS, Université de Bourgogne, FR-21078 Dijon, France ‡
S Supporting Information *
ABSTRACT: Calcium Silicate Hydrates (C−S−H) are the major hydration products of portland cement paste. The accurate description of acid−base reactions at the surface of C−S−H particles is essential for both understanding the ion sorption equilibrium in cement and prediction of mechanical properties of the hardened cement paste. Ab initio molecular dynamics simulations at the density functional level of theory were applied to calculate intrinsic acidity constants (pKa’s) of the relevant SiOH and CaOH2 groups on the C−S−H surfaces using a thermodynamic integration technique. Ion sorption equilibrium in C−S−H was modeled applying ab initio calculated pKa’s in titrating Grand Canonical Monte Carlo simulations using a coarse-grained model for C−S−H/solution interface in the framework of the Primitive Model for electrolytes. The modeling results were compared with available data from electrophoretic measurements. The model predictions were found to satisfactorily reproduce available experimental data.
■
INTRODUCTION The macroscopic behavior of nanoparticle dispersions, from diluted gels and liquids to concentrated pastes, and of nanomaterials is controlled by their interfacial properties.1 The surface charge density, among all other parameters, plays a key role in controlling the behavior of nanoparticle dispersions. For example, the sol−gel transition, the colloidal stability, and self-assembly of nanoparticles as well as the reactivity (dissolution and growth rate), morphology, and size of particles are known to be strongly dependent on their surface charge density.2−8 In other words, a large variety of macroscopic properties such as color and rheology of particle dispersions as well as sequestration capacity of ionic pollutants and cohesion strongly depend on the surface charge density.9−15 Among the most widely used manmade materials, hydrated cement is a great example where electrostatics and thus surface charge density play a fundamental role. As much as 1 ton of concrete (mix of cement, gravels and water) per person is produced every year, causing a great impact on the environment due to the high energy consumption and large CO2 emission during clinker production (i.e., burning of limestones and clays to produce hard nodules of tri- and dicalcium silicates). Thus, there is an environmental and economical need to improve cement manufacturing and use efficiency. Hydrated cement can be described best as a nanomaterial.16 It is constituted by a cohesive and rather dense three-dimensional network of nanoparticles, acting as the glue in concrete, sticking the gravels and rocks together to form a consolidated solid material used in most common housing and civil infrastructure © 2014 American Chemical Society
constructions. This nanoparticle network is formed in a heterogeneous dissolution−nucleation−growth process driven by dissolution of anhydrous cement clinker grains and precipitation of nanoparticle hydrates from a highly supersaturated pore solution.17,18 Among the hydrate phases formed during cement hydration, nanoplatelets of calcium−silicate− hydrates (C−S−H) are dominant (>60%). It is generally accepted that C−S−H is the main source of cohesion in cement and thus responsible for its final hardening and setting.19 Strong short-ranged attractive interactions between single C−S−H platelets have been shown to be the main reason for cohesion in cement paste.9,11 This strong attraction happens to be mainly driven by electrostatics, as a result of the strong surface charge density carried by the C−S−H particles and the high concentration of divalent calcium counterions present in the interstitial pore solution.12 Under such conditions, the fluctuations of the calcium ions between the negatively charged C−S−H particles have been shown to give rise to so-called ion−ion correlation attraction.9,10 One way to better control the setting of cement and to devise new environmentally friendly cementitious materials is thus to obtain a better understanding of the surface charge formation at the C−S−H surface and to describe the physicochemical properties of the C−S−H/water interface at a microscopic scale by means of molecular modeling. Received: March 13, 2014 Revised: April 14, 2014 Published: May 7, 2014 11752
dx.doi.org/10.1021/jp502514a | J. Phys. Chem. C 2014, 118, 11752−11762
The Journal of Physical Chemistry C
Article
C−S−H is made up of platelets of typical dimensions 60 × 30 × 5 nm, which at higher level are further assembled in stacks of up to 100 nm.16 The chemistry and structure of C−S−H phases depend on the synthesis route and the composition of the equilibrium pore solution.11,20 C−S−H platelets are built by calcium layers flanked on each side by linear silicate “dreierketten” chains, as revealed by MAS NMR measurements.21−23 These building blocks have apparent similarities with the structure of the natural minerals tobermorite24−26 and jennite,27,28 based on which models for C−S−H have been developed.29 C−S−H builds up a negative charge at its surface due to the ionization of surface silanol groups SiOH + H 2Oaq ↔ SiO− + H3O+aq 10−pKa = K =
consistently at the same level of the electronic structure theory. This allows to accurately describe the heterogeneous environment at the interfaces, explicitly taking into account, e.g., polarization effects. The method has been successfully applied to different solid/liquid interfaces of relevance to geochemistry, including quartz, alumina, and clays.43,44,46 In particular, this approach has permitted, e.g., to explain the anomalous behavior of silica identifying two different groups with different reactivity at the silica/water interface.43 Applying this simulation technique to the structurally different surface groups the intrinsic acidity constants can be obtained and directly used in the GCMC model for C−S−H. In this paper, we apply ab initio simulations to study the structural and acid/base properties of the C−S−H/water interface at an atomistic level. The obtained structural and thermodynamic properties of the C−S−H surface, e.g., site density distribution and intrinsic pKa values of SiOH and Ca(OH)2 sites were confronted with electrophoretic mobility measurements on C−S−H dispersions in a wide range of experimental conditions using GCMC simulations at the level of the Primitive Model.
(1)
aSiO−a H3O+ a H2Oa≡ SiOH
(2)
where ai are activities. The ionization degree of the C−S−H surface depends on the intrinsic dissociation constants of the silanol groups pKa, surface site density, proton activity aH3O+, and composition of the pore solution. The strong negative surface charge density of C−S−H and the accumulation of Ca2+ ions at the C−S−H/solution interface has been evident from electrophoresis and potentiometric titration experiments.30 Grand Canonical Monte Carlo (GCMC) simulations in the framework of the Primitive Model of electrolyte solutions have been successfully applied to model the development of surface charge density and formation of the electric double layer at the C−S−H/electrolyte interface in a wide range of experimental conditions.31 The C−S−H surface was modeled employing a regular distribution of titration sites and assuming the protolisys constant of all surface groups to be equal to the first dissociation constant of silicic acid. Despite the simplicity, this model was found to reproduce experimental data over a wide range of electrostatic coupling, from a weakly charged surface in contact with a reservoir containing monovalent ions to a highly charged one in contact with a reservoir with divalent ions. In particular, these simulations successfully predict cohesion between C−S−H particles9,10 as well as the charge reversal of C−S−H observed through electrophoresis, surface charge formation, and adsorption of sodium cations at the surface of C−S−H.31,32 Although very successful, the model can still be further refined by including molecular information about the C−S−H surface and by evaluating the acidity of titration sites. Such information is generally difficult to obtain directly in experiments but can be, instead, evaluated from appropriate molecular simulations. A number of approaches have been suggested to correlate the structural properties of surface OH groups with their acidity based on data for small aqueous molecules.33−36 Such approaches fail to take into account proper reorganization of solvent molecules at the interface and are thus less reliable for complex heterogeneous surfaces. Already for a quite simple system, such as the silica surface in contact with water, the surface pKa is quite different from the pKa of silicic acid in solution.37−39 Recently, a direct method for calculating acidity constants of molecules and interfaces based on molecular dynamics simulations at the density functional level of theory has been developed. In this approach, the free energy of titration reactions (for sites at surfaces or in molecules) is obtained by thermodynamic integration.40−45 Remarkably, solute (interface) and solvent are considered
■
MOLECULAR SIMULATIONS System Setup. The structure of C−S−H platelets is consistent with the crystal structure of 11 Å tobermorite.47 The variation of Ca:Si ratio in C−S−H is explained by depolymerization of Si chains and by sorption of Ca ions in the interlayer space and onto the external surfaces of C−S−H platelets. In our simulation setup, the surface of C−S−H was built based on a single structural layer of tobermorite (Figure 1). We consider two limiting cases for the state of the Si chains, namely fully polymerized and fully depolymerized chains. Fully polymerized chains are formed by alternating two pairing and one bridging Si tetrahedra. The protons were attached to apical oxygen sites of the bridging tetrahedra (SibO1H) and to the
Figure 1. (Left) Snapshot from the ab initio molecular dynamics simulations of pKa constants of OH groups on the C−S−H surface. Oxygen atoms are red. Hydrogen atoms are indicated by empty circles. Calcium atoms are gray. Silicon atoms are light brown. (Right) Schematic view on the C−S−H surface with the investigated surface OH sites shown in blue. The hydrogen atoms participating in the deprotonation reactions are colored in yellow. Arcs with arrows indicate angles restrained during the simulations (see method and Table 1 for details). 11753
dx.doi.org/10.1021/jp502514a | J. Phys. Chem. C 2014, 118, 11752−11762
The Journal of Physical Chemistry C
Article
Table 1. Parameters for the Restraint Potential (eq 7) Used in the Simulations r0i (Å) kri (eV Å−2) θ0j (deg) kθj (eV rad−2)
SibO1H
SibO2H
SipO1H
SipO2H
CaOH2
H3O+
H2O
1.01 27.21 120 2.721
1.01 27.21 114 2.721
1.01 27.21 116 2.721
1.01 27.21 116 2.721
0.99 27.21 105 2.721
0.99 27.21 111 2.721
0.99 27.21
oxygen sites at Si−(O)−Ca linkages (SibO2H). Fully depolymerized chains were represented by pairs of Sitetrahedra in which the terminating oxygen atoms were the titrating sites (SipO1H and SipO2H). Use of one and the same simulation setup is essential for consistent calculations and comparison of the pKa values for different titration sites. Both polymerized and depolymerized Si chains were maintained in a single simulation setup as follows. On one side of the Ca layer there are fully polymerized Si chains of tobermorite, and on the opposite side of the layer only depolymerized pairs of tetrahedra are present. This allows determination of deprotonation constants for distinct structural sites using a single setup. The composition of the half layer with polymerized chains was 4 × (Ca2[Si3O7](OH)2 × H2O). The composition of the corresponding half-layer with depolymerized chains was 4 × (Ca2[Si2O5](OH)2 × 2H2O). The interlayer region was filled with 158 water molecules. The dimensions of the orthorhombic simulation supercell were 11.26 × 14.56 × 31.02 Å, where a and b vectors were taken equal to the lattice parameters of natural tobermorite,48 whereas the c spacing was obtained from an NPT preequilibration run using a classical force field.49 Ab Initio Calculations. All ab initio calculations were performed on the basis of the density functional theory using the Gaussian and Plane Waves method (GPW) as implemented in the CP2K simulation package.50,51 The Kohn−Sham orbitals were expanded using a linear combination of atom-centered Gaussian-type orbital functions. In this study, a “short range” double-ξ valence polarized basis set for each atom kind was used.52 This basis set was developed for simulations of condensed matter systems and was proved to give accurate structural and thermochemical results.52 The electron charge density was expanded using an auxiliary basis set of plane waves up to a 320 Ry cutoff. The BLYP exchange and correlation functional53,54 used in this work is known to reproduce the structural properties of water accurately.55 Dual space normconserving pseudopotentials56 were applied to avoid explicit consideration of core electrons. Dispersion interaction was taken into account using the DFT+D2 method.57 Ab initio MD simulations based on the Born−Oppenheimer approximation were performed with a time step of 0.5 fs at 330 K using the Nose−Hoover thermostat.58,59 Such a slightly elevated temperature was used to prevent the commonly reported nonergodic behavior of BLYP water at ambient conditions.60 The applied simulation protocol (330 K and BLYP+D2) is the recommended approach to obtain the best agreement between structural and dynamics properties of DFT and real water at ambient conditions (300 K, 1 bar).61 Before each force evaluation step, the energy was converged to within a value62 of 6 × 10−10 au/atom using a single k-point in the origin of the Brillouin zone (Γ-point sampling). This set of potentials and basis sets has been used successfully in a number of similar systems.43,44,46,63,64 Intrinsic pKa Calculations. Free energies ΔdpA of proton transfer reactions like eq 1 were derived on the basis of vertical
energy gaps in ab initio molecular dynamics simulations at the DFT level of theory40−42 by the thermodynamic integration ΔdpA =
∫0
1
dλ⟨ΔdpE⟩H(λ)
(3)
where ΔdpE is the potential energy of the proton-transfer reaction. The simulations were performed in an ensemble defined by a λ-dependent Hamiltonian Hλ = (1 − λ)H0 + λH1
(4)
where H0 and H1 are the Hamiltonians describing the system of adduct and product in eq 1, respectively. The equilibrium constant of the deprotonation reaction pKa (eq 1, 2) was calculated at the low temperature limit neglecting the difference in the zero-point vibration energy of the proton in SiOH and H3O+aq complexes as 2.3kBT pK a = ΔdpA + kBT ln[c 0 Λ 3H+]
(5)
where kB is the Boltzmann constant, T is temperature, c0 = 1 mol L−1 is the standard state concentration, and ΛH+ is the thermal wavelength of a proton. The last term in eq 5 is a thermochemical correction constant equal to −3.2 pK units. This term accounts for the configurational entropy of the proton in water.40 The integral in eq 3 was approximated by the 5-point Simpson formula: 7 ΔdpA ≈ ΔdpA̅ = [⟨ΔdpE⟩0 + ⟨ΔdpE⟩1] 90 32 12 [⟨ΔdpE⟩0.25 + ⟨ΔdpE⟩0.75 ] + + ⟨ΔdpE⟩0.5 90 90 (6)
Thus, to calculate a single pKa constant five independent trajectories described by a Hamiltonian for adducts, products and a weighted mixture of both have to be obtained. The vertical deprotonation gap for the adducts (λ = 0) was obtained by transferring the proton at the titrating site to a virtual site at the selected H2O molecule in the middle of the interlayer. The system size is chosen big enough to ensure a bulk like behavior of water in the middle of the interlayer. Similarly, for the simulations of reaction products (λ = 1) a proton of the H3O+ in the middle of the interlayer is transferred to a virtual site at deprotonated surface. The dynamics of the virtual sites is defined by an additional harmonic potential (eq 7). This harmonic potential restrains the virtual site to a proper OH distance to ensure an efficient sampling Vrstr
1 = 2
nr
∑ i=1
kir(ri
−
ri0)2
1 + 2
nj
∑ kjθ(θj − θj0)2 (7)
j=1
r0i ,θ0j
where the equilibrium bond distances and angles and corresponding force constants kri ,kθj are summarized in Table 1. The equilibrium values were obtained from simulations of an unrestrained surface. The harmonic potential has a negligible influence on the dynamics of the reference ensemble due to the 11754
dx.doi.org/10.1021/jp502514a | J. Phys. Chem. C 2014, 118, 11752−11762
The Journal of Physical Chemistry C
■
Article
COARSE-GRAINED SIMULATIONS Atomistic structural information about the distribution of titration sites on the C−S−H surface and their pKa values obtained from AIMD simulations were incorporated in a coarse-grained model for C−S−H31 to simulate the charging process of the C−S−H surface and the electrokinetic potential of the C−S−H/water interface. The obtained results were then compared with experimental data to test the performance of our multiscale modeling approach. Model. The coarse-grained model for the C−S−H surface is illustrated in Figure 2. It consists of an infinite planar wall with
small mass of the virtual H-site compared to the mass of the oxygen to which the site is attached. Restraints on solvating water molecules were also necessary to prevent the spontaneous reprotonation of titrating sites if the pKa of the studied sites was larger than the pKa of water. Because these restraints are present in both the initial (protonated) and final (deprotonated) states, the contribution of the imposed potential to the calculated vertical deprotonation gap cancels.43 Preliminary analysis of the surface dynamics has shown that the dihedral angles of the OH sites can have a polymodal distribution which cannot be easily described with simple harmonic restraints. We therefore did not use restraints on dihedral angles as is often done44 and instead performed weighted sampling of the trajectories for λ = 1 as described in the Supporting Information. The length of the sampling trajectory varied from 10 to 18 ps depending on the convergence of the target quantities, preceded by at least 2 ps of equilibration. The initial configuration for the ab initio runs was taken from a pre-equilibration in the NPT ensemble using a classical force field.49 Uncertainties. A number of technical issues related to the thermodynamic meaning of the Helmholtz free energy in eq 3 and the practical thermodynamic definition of the protolysis constant of an OH group were painstakingly discussed in ref 40. The correct calculation of pKa values is by far determined by the accurate treatment of the proton’s solvation energy implicit in eq 5.41 The thermochemical correction term in eq 5 includes the leading contribution of the proton’s translational entropy. Other potential sources of uncertainties which were not corrected for are differences in the zero-point vibration energies of a proton attached to the surface site and in a hydronium ion, system size effects, and a possible mismatch in frequencies of the virtual and real protons. Equation 5 also neglects the difference between Gibbs and Helmholtz free energies of deprotonation, which is expected to be small, however. Another source of uncertainty which often remains unmentioned is the difference in thermodynamic properties between BLYP and real water. Comparative studies of different GGA−water models revealed a glassy behavior of BLYP water at ambient conditions, suggesting that its phase diagram is shifted compared to that of real water.60,65 It has been argued that the best agreement between simulated and experimentally measured properties of bulk water (e.g., density, structure, diffusion) is obtained applying the BLYP exchange-correlation functional augmented with dispersion correction57 and performing the simulation at slightly elevated temperature (330 K).61 Using such a simulation protocol, we thus expect to obtain a good estimate of ΔdpA298K by eq 6. In contrast, T = 298 K has to be used in eq 5 to obtain acidity of the surface OH groups at ambient conditions. The maximum error associated with the use of T = 330 K rather than T = 298 K in eq 5 does not exceed 11%, which is less than the error associated with the sampling. Since one and the same simulation setup was used to calculate the acidity of the different surface sites, it can therefore be expected that the above-mentioned uncertainties, if at all significant, result in one and the same shift for all calculated pKa values. Thus, we can expect the simulations to provide reliable relative acidity constants of the titration sites. Experimental data can be used to correct for such a potential systematic error in the calculated pKa values. The accuracy of the later was analyzed by comparing the results for ⟨ΔdpE⟩λ from the first and second halves of the AIMD trajectories.
Figure 2. Snapshot of a simulation box used in the coarse-grained model of the C−S−H surface. The pairs of titrating sites, modeling the silanols groups, are explicitly treated. Simulation box used in the GCMC simulations showing a model C−S−H surface in contact with a Ca(OH)2-CaCl2 salt solution. Ca2+ ions are shown in blue, Cl− ions are green, and OH− are orange. The dimensions of the box are Lx*Ly*Lz. Periodic boundary conditions are only applied along the x and y directions.
explicit titrating sites. These sites are distributed by pairs on lines, representing silicate chains. The chains are separated by a distance of 5.65 Å. The site pairs represent either silanol groups of bridging tetrahedra, i.e., SibO1H and SibO2H, or two silanol groups belonging to neighboring pairing tetrahedra facing each other, i.e., SipO1H and SipO2H (Figure 1). In tobermorite, the distance between oxygen sites of SibO1H and SibO2H groups, and of SipO1H and SipO2H groups, is about 2.6 Å. However, the height of the oxygen sites in SibO1H and SibO2H groups is not the same. In the coarsegrained model all sites were distributed in plane and the distances between sites in a pair were all set to 2.5 Å, for the sake of simplicity. The distance between centers of pairs was 4.9 Å. Such site distribution reproduces the surface site density of tobermorite, i.e., 4.82 sites nm−2. The C−S−H surface is in equilibrium with an electrolyte solution confined in an orthorhombic box with dimensions Lx × Ly × Lz. Periodic boundary conditions were applied in the x and y directions. In the z-direction the cell was confined by the C−S−H surface on one side and a neutral, impenetrable wall on the opposite side. In addition, the simulations were also performed for a system with equidistant distribution of the surface sites with the same surface site density (4.82 sites nm−2) to investigate the effects of the sites distribution on the titration process. In the primitive model used in this work, each ion i was assigned a valence charge qi and a hard sphere radius Ri. The solvent was approximated by a dielectric continuum with the relative permittivity εs. Given the elementary charge e and the absolute permittivity of vacuum ε0, the interaction between two 11755
dx.doi.org/10.1021/jp502514a | J. Phys. Chem. C 2014, 118, 11752−11762
The Journal of Physical Chemistry C
Article
activity variable, pH = −log10aH+ (Note that in the coarsegrained model there is no distinction between H+ and H3O+ (aq) in eq 1). The model further employs the usual thermodynamic convention for dilute solutions, aH2O = 1. The surface ionization is modeled in two steps. First, a proton (H+) is transferred from a surface site to the simulation box, followed by a transfer of a salt pair HB, where B− is an anion, from the simulation box to the bulk. The corresponding Boltzmann factor for the trial energy is given by
charged species i and j separated by a distance rij is calculated according to ⎧ ⎪ u el(r ) = ij ⎪ ⎪ ⎨ hs ⎪ u (rij) = ⎪ ⎪ u hs(rij) = ⎩
qiqje 2 4πε0εsrij ∞ , if rij ≤ R i + R j 0, if rij > R i + R j
(8)
exp( − βΔU ) =
The hard sphere radii of all ions in the system were set to 2 Å. A one-body external field was used to correct for the longrange electrostatic interactions and to account for the cell boundary constraints along the z-direction acting as hard walls, following the charged sheets method66,67 ext ⎧ ⎪ v (zi) = ∞ , when zi ≥ Lz ||zi ≤ 0 ⎨ ext ⎪ ext ⎩ v (zi) = qieφ (zi), when 0 < zi < Lz
exp[ln 10(pH − pK ai )]
U
No + Ns
=
∑ i
ex
v (zi ) +
No + Ns
∑ i
