THE JOURNAL OF CHEMICAL PHYSICS 142, 054503 (2015)

The wet solidus of silica: Predictions from the scaled particle theory and polarized continuum model G. Ottonello,1,a) P. Richet,2 and M. Vetuschi Zuccolini1

1 2

DIPTERIS, Università di Genova, Corso Europa 26, 16132 Genoa, Italy Institut de Physique du Globe, Rue Jussieu 2, 75005 Paris, France

(Received 30 November 2014; accepted 12 January 2015; published online 3 February 2015) We present an application of the Scaling Particle Theory (SPT) coupled with an ab initio assessment of the electronic, dispersive, and repulsive energy terms based on the Polarized Continuum Model (PCM) aimed at reproducing the observed solubility behavior of OH2 over the entire compositional range from pure molten silica to pure water and wide pressure and temperature regimes. It is shown that the solution energy is dominated by cavitation terms, mainly entropic in nature, which cause a large negative solution entropy and a consequent marked increase of gas phase fugacity with increasing temperatures. Besides, the solution enthalpy is negative and dominated by electrostatic terms which depict a pseudopotential well whose minimum occurs at a low water fraction (XH2O) of about 6 mol. %. The fine tuning of the solute-solvent interaction is achieved through very limited adjustments of the electrostatic scaling factor γel which, in pure water, is slightly higher than the nominal value (i.e., γel = 1.224 against 1.2), it attains its minimum at low H2O content (γel = 0.9958) and then rises again at infinite dilution (γel = 1.0945). The complex solution behavior is interpreted as due to the formation of energetically efficient hydrogen bonding when OH functionals are in appropriate amount and relative positioning with respect to the discrete OH2 molecules, reinforcing in this way the nominal solute-solvent inductive interaction. The interaction energy derived from the SPT-PCM calculations is then recast in terms of a sub-regular Redlich-Kister expansion of appropriate order whereas the thermodynamic properties of the H2O component at its standard state (1-molal solution referred to infinite dilution) are calculated from partial differentiation of the solution energy over the intensive variables. C 2015 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4906745] I. INTRODUCTION

We have recently shown that the solubility of gaseous species in silicate melts may be precisely assessed through an application of the Scaling Particle Theory (SPT) coupled with a Polarized Continuum Model (PCM) assessment of solute-solvent interactions.1 The investigation, based on the experimentally observed solubility of inert gases, showed that two kinds of solvents as different as water and silicate melts may be treated within the framework of a unifying first principles theory. Within the same theoretical framework, it thus appears now possible to investigate the mixing properties of H2O and SiO2 under various Pressure (P), Temperature (T) conditions in the liquid aggregation state and to depict the nature of solute-solvent interactions in this system.

II. METHODS A. The solubility of molecular H2O

The chemical potential of a solute “i” in a liquid solvent can be expressed as2 ( ) Λ3 ¯ i + PV ¯ i + RT ln * i + + RT ln ni , µi,soln = −U V , ji -

(1)

a)Author to whom correspondence should be addressed. Electronic mail:

[email protected] 0021-9606/2015/142(5)/054503/9/$30.00

¯ i and −U ¯ i are the partial molar volume and the where V molar potential energy in solution of the ith solute relative to infinite separation, respectively, Λ3i is the momentum partition function, ji is the internal degree of freedom of the solute, and ni is the number of solute molecules in the solvent volume. As shown by Pierotti,3 the first two terms on the right-hand side of Eq. (1) may be replaced by the summation of the cavitation ¯ cav + G ¯ int, energy and of the solute solvent interaction G ( ) Λ3 ¯ cav + G ¯ int + RT ln * i + + RT ln ni . µi,soln = G (2) V , ji The chemical potential of the solute in a gas phase in equilibrium with the solution is given by ( ) Λ3 fi µi,gas = RT ln * i + + RT ln , (3) RT , ji where fi is the fugacity of the solute. Assuming the momentum and translational partition function of the ith component to be identical in the two aggregation states and equating the chemical potentials, one gets ( ) ( ) ¯ cav G ¯ int fi G RT ln = + + ln . (4) ni RT RT Vsln Considering one mole of liquid composed of solute i and solvent s (i.e., ni + ns = 1), we can substitute molar fraction xi

142, 054503-1

© 2015 AIP Publishing LLC

054503-2

Ottonello, Richet, and Vetuschi Zuccolini

J. Chem. Phys. 142, 054503 (2015)

for the number of moles. For a simple binary system, we have ( )   ¯ cav G ¯ int fi G RT ln = + + ln . (5) ¯ i + (1 − xi) V ¯s xi RT RT xiV ¯ sln  V ¯s When the ith component is sufficiently diluted V and the right-hand side term reduces to( the ) Henry’s law constant of the ith solute in solution, ln xfi = lnKH,i. Note, i however, that our equation is more generally valid and applies to any concentration of “i” in solution. Also, purposely, it does not contain any configurational contribution arising from mixing in solution terms. Accounting for such additional contributions requires the adoption of model-dependent inverse methods. Assuming for instance a nonreactive solution (i.e., an amount of dissolved H2O equal to the nominal concentration) and completely random mixing in liquid, it is found (see below) that the magnitude of the second term on the right-hand side of the solution equation is marginally modified without affecting the model accuracy. B. The cavitation, electronic, dispersive, and repulsive energies

¯ cav required to The partial molar Gibbs free energy G obtain a cavity in a fluid of hard spheres was obtained by Reiss and coworkers4–7 who relied on a statistical mechanical theory of fluids based upon the properties of exact radial distribution functions. The interaction energy in all electron ab initio calculations based on the PCM10,11 is conceived as the summation of electrostatic dispersive and ¯ int = G ¯ el + G ¯ disp + G ¯ rep. repulsive solute/solvent interactions G 12,13 In the classical approximation, the sum of solute-solvent dispersive and repulsive interactions Gdis + Grep is computed as double summations of terms extending to the solute atoms placed within the cavity and solvent atoms outside a cavity limited by a tessellated surface, Gdis + Grep =

solute  solvent  a

s

ρs

tesserae 

aiΓas,

(6)

i

where ρs is the number density of the solvent in Å−3, ai is the area of the ith tessera composing the surface, and Γas is a set of dispersive and repulsive parameters obtained from crystallographic data.14 III. RESULTS

We sampled a few points along the wet solidus of the H2OSiO2 system delineated by experiments15,16 and previously adopted by us to illustrate the bulk solubility behavior in silica liquid.17 The selected compositions span almost the entire solubility curve from the nominally anhydrous low-P melting limit to the critical point located at a pressure of ∼1 GPa. Two additional points were added to represent the compositional limits of the system (i.e., pure H2O and H2O at infinite dilution in SiO2) at 25 ◦C and 1 bar. We adopted for the cavitation energy the formulation of Pierotti8,9 and performed calculations for the interaction energy with the PCM model. Ab-initio calculations were carried out with the GAUSSIAN computer package18 at the B3LYP/6-31+G(d,p) theory level. Diffuse and polarized functions were adopted to

ensure accurate sampling of the electron density in the outer region. Results do not change appreciably when the diffuse and polarization components of the basis set (6-31G) are omitted, however, or when calculations are performed at a simple Hartree-Fock level with Slater-type Orbitals (HF/STO-3G). The various solvent parameters are summarized in Table I. The dielectric constant was obtained through application of the Clausius-Mossotti equation1 and the solvation diameter σs of the mixed solvents was derived from a simple expression of the form σS = 2 × (1.1383rcoll√− 0.3836), where rcoll is the collisional radius1,19 (i.e., rcoll = 3 V/8). For the solute radius, we used the United Atom Topological Model applied on radii optimized for the HF/6-31G(d) theory level (UAHF). These are the recommended radii for the calculation of the solvation Gibbs free energy (∆Gsolv) when the calculation of the solvation terms is preceded by a gas-phase optimization of the cluster. For comparative purposes, we also performed the calculations with United Atom Topological Model (UA0) applied on atomic radii of the Universal Force Field (UFF). In this case, hydrogens are enclosed in the sphere of the heavy atom to which they are bonded. The number density of the mixed solvent under room conditions and its thermal expansion at T of interest was based on molar volumes calculated with the procedure of Bottinga and coworkers.20 The isothermal compressibility was based on the experimental findings of Rivers and Carmichael.21 In the ab initio calculation of the electrostatic interaction, one must finely tune the electrostatic radius of the solute, which may appreciably differ from the nominal radius adopted for conforming the cavity (σi/2). The default electrostatic scaling factor γel adopted for UAHF is 1.2, whereas a default value of 1.0 is adopted for UA0.18 Because not all dissolved water remains in molecular form OH2 but partly reacts with the silica polymers to form hydroxyl terminations, tuning was done to ensure agreement of the calculated molar fraction of molecular H2O with experimental observations (or to the experimental gas-phase fugacity for pure water at room conditions). To obtain the amount of molecular OH2 in the H2O-SiO2 liquid at the selected compositions, we adopted a model22 based on the simple assumption that the reactive solubility of H2O in silica is well represented by an homogeneous equilibrium of type [OH2]melt + [O]0melt ⇔ 2[OH]0melt

(7)

with a T,P independent constant Kreactive = 0.11. The model returns correct bulk solubilities up to a maximum of ∼6% in weight of H2O. The gas phase fugacity of water was computed with the SUPERFLUID package23 which spans all the P,T conditions of interest with the appropriate accuracy. In Table II we summarize the SPT calculations in terms of Gibbs free energy (Gcav) and enthalpic, entropic, and volumetric components (Hcav, Scav, Vcav). The values of Gcav are listed for not only the reference temperature Tr = 298.15 K but also for the temperature of interest T to emphasize the dominant role of T on the solubility of OH2. The Henry’s law constant is obviously related to the Gibbs free energy of solution (∆Gsln), ∆Gsln = lnKH.

(8)

Theory

UA0 UAHF UAHF UAHF UAHF UAHF UAHF UA0 UAHF UAHF UAHF UAHF UA0 UAHF UAHF UAHF UAHF UA0 UAHF UAHF UAHF UAHF UA0 UAHF UAHF UAHF UAHF UA0 UAHF UAHF UAHF UAHF UA0 UAHF UAHF UAHF UAHF UA0 UAHF UAHF UAHF

DFT DFT DFT DFT HF DFT DFT DFT DFT DFT DFT DFT DFT DFT DFT DFT DFT DFT DFT DFT DFT DFT DFT DFT DFT DFT DFT DFT DFT DFT DFT DFT DFT DFT DFT DFT DFT DFT DFT DFT DFT

Basis set 6-31+G(d,p) 6-31+G(d,p) 6-31G(d,p) 6-31G STO-3G 6-31+G(d,p) 6-31+G(d,p) 6-31+G(d,p) 6-31+G(d,p) 6-31+G(d,p) 6-31+G(d,p) 6-31+G(d,p) 6-31+G(d,p) 6-31+G(d,p) 6-31+G(d,p) 6-31+G(d,p) 6-31+G(d,p) 6-31+G(d,p) 6-31+G(d,p) 6-31+G(d,p) 6-31+G(d,p) 6-31+G(d,p) 6-31+G(d,p) 6-31+G(d,p) 6-31+G(d,p) 6-31+G(d,p) 6-31+G(d,p) 6-31+G(d,p) 6-31+G(d,p) 6-31+G(d,p) 6-31+G(d,p) 6-31+G(d,p) 6-31+G(d,p) 6-31+G(d,p) 6-31+G(d,p) 6-31+G(d,p) 6-31+G(d,p) 6-31+G(d,p) 6-31+G(d,p) 6-31+G(d,p) 6-31+G(d,p)

Hybrid ... sp3 sp3 sp3 sp3 sp3 sp3 ... Sp3 sp3 sp3 sp3 ... Sp3 sp3 sp3 sp3 ... sp3 sp3 sp3 sp3 ... sp3 sp3 sp3 sp3 ... Sp3 sp3 sp3 sp3 ... Sp3 sp3 sp3 sp3 ... sp3 sp3 sp3

γel

σs

ν × 102

α0 × 105

α1 × 109

β

ζ

ε

Gint

Gel

Gdisp

Grep

1.0000 1.2000 1.2000 1.2000 1.2000 1.0000 1.2240 1.0000 1.2000 1.0000 1.0556 1.0925 1.0000 1.2000 1.0000 1.0557 1.0903 1.0000 1.2000 1.0000 1.0487 1.0770 1.0000 1.2000 1.0000 1.0290 1.0527 1.0000 1.2000 1.0000 1.0040 1.0218 1.0000 1.2000 1.0000 0.9864 0.9958 1.0000 1.2000 1.0000 1.0945

2.77 2.77 2.77 2.77 2.77 2.77 2.77 3.091 3.091 3.091 3.091 3.091 3.136 3.136 3.136 3.136 3.136 3.167 3.167 3.167 3.167 3.167 3.193 3.193 3.193 3.193 3.193 3.211 3.211 3.211 3.211 3.211 3.232 3.232 3.232 3.232 3.232 3.262 3.262 3.262 3.262

3.332 3.332 3.332 3.332 3.332 3.332 3.332 2.569 2.569 2.569 2.569 2.569 2.48 2.48 2.48 2.48 2.48 2.422 2.422 2.422 2.422 2.422 2.375 2.375 2.375 2.375 2.375 2.343 2.343 2.343 2.343 2.343 2.306 2.306 2.306 2.306 2.306 2.254 2.254 2.254 2.254

27.67 27.67 27.67 27.67 27.67 27.67 27.67 9.301 9.301 9.301 9.301 9.301 6.692 6.692 6.692 6.692 6.692 4.999 4.999 4.999 4.999 4.999 3.611 3.611 3.611 3.611 3.611 2.682 2.682 2.682 2.682 2.682 1.616 1.616 1.616 1.616 1.616 0.1 0.1 0.1 0.1

−66.049 −66.049 −66.049 −66.049 −66.049 −66.049 −66.049 −7.519 −7.519 −7.519 −7.519 −7.519 −4.042 −4.042 −4.042 −4.042 −4.042 −2.312 −2.312 −2.312 −2.312 −2.312 −1.23 −1.23 −1.23 −1.23 −1.23 −0.687 −0.687 −0.687 −0.687 −0.687 −0.254 −0.254 −0.254 −0.254 −0.254 0 0 0 0

0.556 0.556 0.556 0.556 0.556 0.556 0.556 4.391 4.391 4.391 4.391 4.391 4.911 4.911 4.911 4.911 4.911 5.26 5.26 5.26 5.26 5.26 5.553 5.553 5.553 5.553 5.553 5.752 5.752 5.752 5.752 5.752 5.984 5.984 5.984 5.984 5.984 6.32 6.32 6.32 6.32

1.408 1.213 1.213 1.213 1.213 1.213 1.213 1.262 1.087 1.087 1.087 1.087 1.244 1.071 1.071 1.071 1.071 1.231 1.061 1.061 1.061 1.061 1.221 1.052 1.052 1.052 1.052 1.215 1.046 1.046 1.046 1.046 1.207 1.040 1.040 1.040 1.040 1.196 1.030 1.030 1.030

78.3 78.3 78.3 78.3 78.3 78.3 78.3 5.4 5.4 5.4 5.4 5.4 4.7 4.7 4.7 4.7 4.7 4.4 4.4 4.4 4.4 4.4 4.1 4.1 4.1 4.1 4.1 3.9 3.9 3.9 3.9 3.9 3.8 3.8 3.8 3.8 3.8 3.5 3.5 3.5 3.5

−52.6 −48.2 −42.5 −49.1 −31.2 −99.2 −45.1 −38.7 −35.7 −71.4 −56.1 −48.9 −36.9 −34.2 −66.9 −52.8 −46.6 −35.8 −33.1 −64.7 −52.6 −47.2 −34.9 −32.3 −62.9 −55.3 −50.8 −34.2 −31.7 −61.6 −60.5 −56.7 −33.5 −31.4 −61.0 −65.1 −62.2 −32.4 −30.1 −58.9 −40.5

−38.3 −32.7 −26.9 −33.6 −15.8 −83.6 −29.5 −27.7 −23.8 −59.5 −44.1 −36.9 −26.4 −22.6 −55.4 −41.2 −35.0 −25.4 −21.8 −53.4 −41.3 −35.9 −24.7 −21.2 −51.8 −44.2 −39.7 −24.2 −20.8 −50.7 −49.6 −45.7 −23.6 −20.5 −50.3 −54.4 −51.5 −22.8 −19.6 −48.4 −30.0

−16.3 −21.7 −21.7 −21.8 −21.8 −21.7 −21.6 −12.5 −16.7 −16.7 −16.7 −16.7 −12.1 −16.1 −16.0 −16.0 −16.0 −11.8 −15.7 −15.6 −15.6 −15.6 −11.6 −15.4 −15.4 −15.4 −15.4 −11.4 −15.2 −15.1 −15.1 −15.1 −11.3 −15.1 −14.9 −14.9 −14.9 −11.0 −14.6 −14.6 −14.6

2.1 6.2 6.1 6.3 6.4 6.2 6.0 1.5 4.7 4.7 4.7 4.7 1.5 4.5 4.4 4.4 4.4 1.5 4.4 4.4 4.4 4.4 1.4 4.3 4.3 4.3 4.3 1.4 4.3 4.2 4.2 4.2 1.4 4.2 4.2 4.2 4.2 1.4 4.1 4.1 4.1

J. Chem. Phys. 142, 054503 (2015)

Radii

1 1 1 1 1 1 1 0.335 0.335 0.335 0.335 0.335 0.244 0.244 0.244 0.244 0.244 0.184 0.184 0.184 0.184 0.184 0.133 0.133 0.133 0.133 0.133 0.098 0.098 0.098 0.098 0.098 0.0583 0.0583 0.0583 0.0583 0.0583 0 0 0 0

Ottonello, Richet, and Vetuschi Zuccolini

XH2O

054503-3

TABLE I. Polarized continuum model calculations of interaction energies (Gint = Gel + Gdisp + Grep) in the H2O–SiO2 binary join at liquid state and T and P of reference (Tr = 298.15 K; Pr = 105 Pa). All energy terms are in kJ/mol. γel is the electrostatic scaling factor; σs is the solvent diameter (Å); ν is the number density (molecules/Å3); α0 and α1 are the isobaric coefficients of thermal expansion (αT = α0 + α1 × T) in K−1 and K−2, respectively; β is the isothermal compressibility in megabars−1; ζ is the solute/solvent diameter ratio; ε is the dielectric constant. Values in italics account for the presence of ideal configurational entropic effects in the liquid mixture computed assuming non-reactive dissolution of molecular H2O.

054503-4

Ottonello, Richet, and Vetuschi Zuccolini

J. Chem. Phys. 142, 054503 (2015)

TABLE II. Scaling particle theory energy terms. The cavitation enthalpy at P of interest and T of reference (Hcav,P,Tr), the interaction energy (Hint,P,Tr, assumed as purely enthalpic), and the cavitation Gibbs free energies (Gcav,P,Tr; Gcav,P,T) are in kJ/mol. The fugacity was calculated with the SUPERFLUID package23 assuming the gas phase to be pure H2O. The cavitation entropy (Scav,P,Tr) is in J/(mol × K). Volume terms are in cc/mol. KPr,Tr and KP,T are the Henry’s law constants for H2O in pure silica liquid at Tr, Pr of reference (Tr = 298.15 K; Pr = 105 Pa) and T,P of interest, respectively. Values in italics account for the presence of ideal configurational entropic effects in the liquid mixture computed assuming non-reactive dissolution of molecular H2O. XH2O P (bar) f (bar) T (◦C) Hcav,P,Tr Scav,P,Tr Vcav,P,Tr Gcav,P,Tr Gcav,P,T Hint,P,Tr Vint,P,Tr Vsolute,P,Tr lnKPr,Tr lnKP,T

1

0.335

0.244

1 1 25 2.9 −52.6 12.2 18.6 18.6 −45.1 −45.1 −0.250 −0.250 11.934 11.934 −3.45 −3.45

10 000 19 610 1 080 13.0 −57.1 13.7 30.0 94.1 −56.1 −48.9 −2.462 −2.147 11.300 11.615 −3.53 11.81

7 500 11 380 1 065 9.7 −57.8 13.9 27.0 89.8 −52.8 −46.6 −2.593 −2.290 11.386 11.689 −3.49 11.73

The equations governing enthalpy of solution, entropy of solution, and partial molar volume of the OH2 solute in the H2O-SiO2 mixed solvent have been discussed elsewhere.1,9 As already outlined, in the present formulation, these relations are valid not only at the infinite dilution condition but also over the entire compositional field.

0.184

0.133

0.098

0.0583

0

5000 2500 1250 500 1 5940 2520 1180 577 1 1070 1100 1170 1300 1723 6.6 3.4 1.8 0.8 1.3 −58.3 −58.7 −59.0 −59.3 −59.7 14.0 14.1 14.2 14.3 14.4 23.9 20.9 19.4 18.5 17.8 86.9 85.5 88.1 94.8 119.2 −52.6 −55.3 −60.5 −65.1 −40.6 −47.2 −50.8 −56.7 −62.2 −40.6 −2.764 −3.069 −3.480 −3.898 −2.565 −2.484 −2.820 −3.259 −3.724 −2.565 11.361 11.179 10.852 10.532 12.007 11.642 11.428 11.073 10.760 12.006 −4.64 −6.98 −9.71 −11.98 −2.35 11.46 11.04 10.73 10.77 13.47

In contrast to the cavitation energy, the interaction energy is parameterizable, namely, in terms of its electrostatic component, which depends upon the scaling factor γel assigned to the solute radius (Table I). Such a scaling has been justified by the fact that the first solvation layer cannot be assigned the same dielectric properties as the bulk solvent.24 In our case, moreover, hydrogen bonding (HB) likely operates between the solute molecule H2O and the hydroxyl functionals of proton

IV. DISCUSSION

Figure 1 shows schematically the cavity depicted for a spherical solute molecule i in a liquid solvent s with a solvation diameter σs. The dotted line is the Solvent Accessible Surface (SAS) depicted by a probe sphere rolling on the surface generated by the interlocking spheres, which represent the atoms with van der Waals atomic radii (or “van der Waals surface,” here coinciding for simplicity with a single sphere) and the dashed line is the envelope of the volume formed from van der Waals atomic radii, scaled by the electrostatic factor γel, and excluded to the rolling probe sphere (usually denoted24 “Solvent Excluding Surface,” SES). The cavitation energy and its components are not parameterizable, because they depend on a few well established quantities, namely, the radius of the solute (1.68 Å according to the SCFVAC-UAHF theory or 1.98 Å when using the UFF-UA0 approach) and the solvation radius of the solvent (intrinsically related to the collisional radius based on the molecular volume of the solvent1). In Figure 2 we may appreciate the effect of the σi/σs ratio on the two limiting compositions (i.e., H2O at infinite dilution in pure H2O and in liquid SiO2) at reference Tr = 298.15 K and Pr = 105 Pa. The composition of the solvent has only a limited effect on the bulk magnitude of the cavitation energy which is dominated by a strong negative entropic term. The cavitation enthalpy is virtually negligible.

FIG. 1. Spherical cavity of radius r defining the SAS (dotted line) originated by a solute of diameter σi in a hard spheres solvent of molecules of radius σ s. The electrostatic scaling parameters γel act on the nominal van der Waals radius of the solute, expanding (or contracting) it and generating the SES (dashed line). Electrostatic scaling operates only on the inductive energy terms.

054503-5

Ottonello, Richet, and Vetuschi Zuccolini

FIG. 2. Gibbs free energy of cavitation for OH2 in pure H2O and (fictive) SiO2 liquid and its enthalpic and entropic components at Tr = 298.15 K, Pr = 105 Pa.

terminated polymers in the liquid when these are present in the complex solvent in sufficient fraction to surround the solute molecule. In Figure 3 we see that the electrostatic radius optimized on the basis of the bulk interaction energy necessary to reproduce the solubility behavior of H2O in the investigated media depicts a pseudopotential well with a minimum localized at rel ∼ 1.65 (configurational effects absent) or rel ∼ 1.67 (ideal random mixing and non-reactive behavior). The minimum energy is attained when the bulk molar amount in the system is XH2O ∼ 0.06 (Figure 4). There

J. Chem. Phys. 142, 054503 (2015)

FIG. 4. Mixing enthalpy of OH2 in H2O-SiO2 liquid. Solid and open circles: discrete values calculated with and without ideal Gibbs free energy of mixing contributions, respectively. In this second case, there are no apparent deviations from ideal configuration (filled triangles). The solid lines encompassing the enthalpy terms are generated by the adopted Redlich-Kister 6th order expansions (Table IV). Although extended to XH2O = 0.4 for graphical purposes, their practical utilization is recommended only up to the consolute limit (XH2O = 0.335).

is no doubt that the electrostatic part of the solute-solvent interaction is well resolved by the perturbed Hamiltonian. However, the description of the dispersive and repulsive part of the interaction energy (Eq. (6)) is still subject to refinement.25 Moreover, though the fine tuning of the scaling factor γel assigned to the solute radius is quite appropriate, the fact that this is nonlinear with respect to composition in the binary system needs some justification. It is moreover almost compulsory to admit that the very large stabilizing effect observed at low H2O contents in the mixed solvent originates in additional bonding terms not easily recast in terms of the current formulation of the dispersion. HB could form whenever OH terminated polymers contour in an efficient way the discrete OH2 molecules without the need of severe distortion. To evaluate the soundness of our hypothesis, we carried out ab initio computations adopting the Swart-Solà-Bickelhaupt functional with Grimme’s dispertion correction (SSB-D) functional as implemented in NEWCHEM.26 The dispersion term is added empirically through long-range contribution Density Functional Theory with Dispersion (DFT-D), i.e., EDFT-D = EDFT − KS + Edisp,

(9)

where Edisp = −s6

Natom atom −1 N i=1

FIG. 3. Electrostatic radius of OH2 in the H2O-SiO2 mixed solvent. The shrinking of the electrostatic radius with the increase of H2O component in the solvent can be ascribed to the formation of hydrogen bonding between H2O and hydroxyl functionals of polymer clusters in solution when the latter are in sufficient amount to surround water molecule in an appropriate position.

ij

C6

  −α(Rij/Rvdw−1) × 1 + e . 6 j=i+1 Rij

(10)

In this equation, the s6 term (s6 = 0.847 455 in our case) ij depends upon the functional and basis set used, C6 is the dispersion coefficient between pairs of atoms, and Rvdw and Rij are related with van der Waals atom radii (hence intimately connected, though in nonexplicit way to γel) and the nuclear

054503-6

Ottonello, Richet, and Vetuschi Zuccolini

J. Chem. Phys. 142, 054503 (2015) FIG. 5. Hydrogen bonding between discrete OH2 molecules and hydroxyl terminations of silica polymers. In position I, molecular water approaches the [Si5O6(OH)8] neutral cluster at an OH functional termination pertaining to a singly hydroxylated monomeric group. In position II, molecular water approaches a doubly hydroxylated unit. In both cases (and also in the case of HB with the orthosilicic acid unit), the HB distance is ∼1.98 Å.

distance, respectively. The α value contributes to control the corrections at intermediate distances.27–31 A good proxy of the tridimensional fully polymerized network in the presence of H2O in the system is the neutral cluster [Si5O6(OH)8]0 (a structure of intermediate Lewis acidity corresponding to the OH terminated analogue of the structure 4 of Civalleri and coworkers32 with three 4membered rings; Figure 5). This structure which most likely forms in the SiO2-rich portion of the H2O-SiO2 binary has 8 hydroxyl terminations, 2 of them pertaining to the two singly hydroxylated monomers. In the intermediate to low acidity range, the system is most likely completely depolymerized and orthosilicic acid units [H4SiO4]0 are present. Our SSB-D calculations return ∆HHB energies of −45.7 kJ/mol when one OH2 molecule approaches the OH functional of the singly hydroxylated monomer and of −33.7 kJ/mol when it approaches the OH functional of the doubly hydroxylated monomer. The ∆HHB energy is −39 kJ/mol when OH2 approaches the orthosilicic acid unit (Table III). These energies are quite consistent with the estimates of the ∆HHB energy in BLYP-D geometrically optimized structures of water33 which returned energies spanning from −20.6 kJ/mol (dimer) to −42.0 kJ/mol (hexamer). If we now compute the electrostatic energy difference for calculations carried out with the refined electrostatic scaling factor with respect to the nominal one (γel = 1.200) (Table I), we obtain at low H2O content elthalpy differences consistent with the formation of an additional HB with respect to the hydrogen bonding state of H2O in pure water (i.e., ∆H = −39 kJ/mol at XH2O  0.06). Moving in the two compositional directions the nonlinear contribution to solvation energy is still evident

but less intense (∆H = −20.3 at XH2O  0.34; ∆H = −18.6 at XH2O  0.24; ∆H = −19.5 at XH2O  0.18; ∆H = −23 at XH2O  0.13; ∆H = −28.8 at XH2O  0.10; ∆H = −10.4 at XH2O  0.00). It should be noted that the [OH2] to [OH] proportion experimentally observed at XH2O  0.06 corresponds to the proportion depicted as “position I” in Figure 5 (i.e., molecular water linked via HB to a singly hydroxylated silica group) and also for hydrogen bonding with the silicic acid group (Fig. 5). The calculated ∆HHB energies (−45.7 and −39 kJ/mol) are nearly similar to the magnitude of the potential well at its bottom (−39 kJ/mol). These calculations, though confirming the likely formation of HB between molecular water and the hydroxyl terminations of the polymeric chains (and/or hydroxylated monomers), are too approximate for a precise assessment of the statistical number of bonds eventually formed at any composition in the system. It must be stressed also that the HB energy contributions, though being able to explain the nature of the pseudopotential well, have no practical effect on the liquid/gas partitioning of molecular water, due to the huge restoring entropy effects involved in complexation. As a matter of fact, it may be readily seen by summing up the zero point and thermal correction to G (ZPC, Therm-G, respectively, in Table III) to the electron energy (Eel in the same Table) that HB in most cases does not survive the thermal stress already at room condition (Table III). V. APPLICATION

Obviously, the complex nature of the reactive solubility of H2O in silicate melts we analyzed for the simple H2O-SiO2

TABLE III. SSB-D functional calculations for hydroxyl-terminated clusters bonded to a water molecule. Eel = electronic energy; ZPC = zero point correction; Therm-U,H,G = thermal corrections to internal energy, enthalpy, and Gibbs free energy at T = 298.15 K; S = entropy. All magnitudes are expressed in hartree if not otherwise specified. Single and double asterisks identify positions I and II, respectively, in Fig. 5. Cluster [Si5O6(OH)8] [OH2] [Si5O6(OH)8(OH2)]∗ [Si5O6(OH)8(OH2)]∗∗ [Si(OH)4] [Si(OH)4(OH2)] HB reaction (T = 298.15 K, P = 105 Pa) [Si5O6(OH)8] + [OH2] [Si5O6(OH)8(OH2)]∗ [Si5O6(OH)8] + [OH2] [Si5O6(OH)8(OH2)]∗∗ [Si(OH)4] + [OH2] [Si(OH)4(OH2)]

Eel

ZPC

Therm-U

Therm-H

Therm-G

S [J/(mol × K)]

−2519.545 040 48 −76.952 333 56 −2596.522 398 20 −2596.518 160 44 −596.251 567 13 −673.226 329 26

0.150 209 00 0.021 671 00 0.176 327 00 0.176 607 00 0.058 275 14 0.084 342 00

0.174 792 00 0.024 505 00 0.203 428 00 0.203 470 00 0.065 386 68 0.094 013 00

0.175 736 00 0.025 448 00 0.204 372 00 0.204 414 00 0.066 330 87 0.094 957 00

0.097 340 61 0.003 382 39 0.120 995 80 0.121 504 86 0.028 930 90 0.050 381 05

690.3 194.3 734.2 730.1 329.3 392.5

∆U (kJ/mol) −43.2

∆H (kJ/mol) −45.7

∆S (J/(mol × K)) −150.4

∆G (kJ/mol) −0.8

−31.2

−33.7

−154.6

12.4

−36.5

−39.0

−131.1

0.1

054503-7

Ottonello, Richet, and Vetuschi Zuccolini

J. Chem. Phys. 142, 054503 (2015)

TABLE IV. Redlich-Kister parameterization of the SPT-PCM enthalpy of mixing in the H2O-SiO2 liquid. Coefficient A0 A1 A2 A3 A4 A5 A6 Mean error (kJ) Standard error (kJ)

Gmixing = XAXB

a Model

(1)

b Model

(2)

Model 1a

Model 2b

−132.860 −898.290 −1142.000 1235.400 1745.200 −3641.200 −4882.100 0.084 0.120

−19.703 −204.120 −740.470 17.163 3188.700 3124.100 −539.320 0.141 0.224

6 

Ai × (XA − XB)i . i=0 +RT [XA ln (XA) + XB ln (XB)]

Gmixing = XAXB

6 

Ai × (XA − XB)i.

i=0

system prevents any sort of application from being made if not translated into a simple thermodynamic assessment. To achieve this goal, it is sufficient to refer to the standard state of “hypothetical 1-molal solution referred to infinite dilution,” commonly adopted in solution chemistry. Repeating the calculations at various temperatures and at Pr = 105 Pa, one may retrieve the appropriate value assigned to the heat capacity of solution, which turns out to be represented by a constant term [Cp,sln = −8.2768 J/(mol × K)]. The enthalpy of solution at the standard state is −43.050 kJ/mol and the entropy of solution is −124.823 J/(mol × K). The molar volume of the solvent is 12.007 cc/mol, in line with the density and Brillouin scattering observations34,35 (actually the standard state interaction was parameterized in such a way as to reproduce this value through the appropriate γel and the ensuing interaction volume; cf. Table II). The parameters at the other limiting condition represented by H2O in solution in a pure H2O solvent are easily obtained from the experimental values of the vapor pressure of liquid H2O at Pr,Tr. Because the cavitation energy depends linearly upon composition at all P,T conditions, the properties of the mixture may be easily recast in terms of a sub-regular Redlich-Kister expansion accounting for the complex behavior of the electrostatic energy terms. To avoid error propagation near the maximum solubility limit, we deemed necessary to adopt a 6th-order expansion whose coefficients are listed in Table IV. In this respect, we stress that the two sets of listed values refer to the hypothetical conditions of (1) a liquid mixture deprived of configurational terms and (2) a liquid mixture where configurational terms are ideal and all water dissolves as molecular H2O. Accounting for the presence of ideal mixing terms slightly modifies the magnitude of the electrostatic interaction and softens the potential well, as shown in Figure 4, but does not affect significantly the model results. This simple thermodynamic assessment is a powerful tool to determine partitioning of H2O between the gas and liquid phases. For practical purposes, we may use it in three distinct ways.

(1) Calculate the amount of molecular H2O in the liquid knowing the P,T conditions of the system (and hence the fugacity of H2O in the gas phase). (2) Refine the P or T of the wet solidus on the basis of the bulk amount of H2O in the liquid. (3) Determine the fugacity of H2O in the gas phase based on the observed OH2-OH disproportionation in the liquid. The first procedure is the one applied in an inverse fashion to parameterize the electrostatic interaction. It needs not to be further considered. Concerning the second procedure, we note that the reactive disproportioning of the H2O component into molecular OH2 and OH hydroxyls is not appreciably affected by T. Assuming the reaction constant to be known, one readily derives the relative fractions of molecular H2O and OH in the liquid.22 We then see in Figure 6 that the imposed speciation is correctly reproduced by the SPT-PCM calculus. Because the SPT-PCM calculations return the appropriate gas/liquid partitioning based on the equality of equilibrium of the molecular form in the two aggregation states, and because the gas phase fugacity may be readily obtained from the compressibility parameter of gaseous H2O,23 there is, for a given amount of bulk H2O in the liquid, a unique value of T that satisfies the equality. In Figure 7 we see the results of this exercise obtained when adopting the nominal reactive constant22 (Kreactive = 0.11). The maximum bias between input T and model refinement is observed at the critical point where the SPT-PCM procedure predicts a temperature about 20 ◦C higher than the experimental value. The third procedure is perhaps the most interesting from a predictive point of view. The SPT-PCM calculation is a unifying conceptual framework for spanning the entire compositional range under the whole P,T intervals of interest. As such, it is not restricted to the loci of P,T conditions representing the wet solidus of silica since the liquidus region and the supercooled liquid may be explored as well (Figure 8). The possibility of exploring the P-T

FIG. 6. Relative abundances of OH2 and hydroxyls in H2O-SiO2 liquid. The homogeneous partitioning is reproduced by a reactive constant Kreactive = 0.11 not sensibly affected by T and P.22 The SPT-PCM discrete calculations are superimposed.

054503-8

Ottonello, Richet, and Vetuschi Zuccolini

J. Chem. Phys. 142, 054503 (2015)

VI. CONCLUSIONS

FIG. 7. Wet solidus of silica. Model results derived from the Redlich-Kister parameterization of the interaction energy are compared with the input values adopted to conform the electrostatic interaction through electrostatic scaling of the UAHF radii of the solute molecule.

trajectory of a compositionally evolving liquid has obvious advantages in materials science. The increase of H2O in the liquid causes an increased fugacity of H2O in the gas phase, if not counterbalanced by a concomitant decrease in temperature which favors the solubility of the same component in molecular form in the liquid. The response in terms of bulk pressure of the system is highly non-linear and unpredictable through simple assumptions unable to depict the complex solubility illustrated by the SPT-PCM procedure here developed. The generalization to chemically complex melts and industrial slags is straightforward.

The solution of gases in polar solvents like water is characterized by large negative heats and negative solution entropies which determine an inverse relationship between temperature and solubility (i.e., the fugacity of the gaseous species increases with the increase of T). Long ago Frank and Evans36 ascribed this behavior to the formation of an ordered hydrogen-bonded structure around the solute molecule (the so called “iceberg theory” or “frozen layer theory”). More recently34 through density and Brillouin scattering measurements it was found that the partial molar volume of H2O at the infinite dilution limit in silicate glasses (∼12.0 cc/mol) has volume properties approximating those of ice VII, the densest form of ice. Incidentally, as shown here, this is also the volume one would observe when H2O dissolves at molecular level (i.e., cavitation + interaction volume). Investigating by first principles calculations the solubility of inert gases in silicate liquids,1 we showed that the solvent properties of chemically complex silicate melts are intermediate between those of a polar solvent like H2O and a non-polar solvent like C6H6, which may be accounted for through a simple operative modification of the collisional theory.19 More importantly, it was confirmed through detailed SPT-PCM calculations that in silicate liquids as well as in H2O, the fugacity of noble gases decreases by several orders of magnitude with the decrease of T due to an increased solubility in the liquid. Stemming from the above considerations, we analyzed in terms of SPT-PCM theory the wet solidus of the H2O-SiO2 system in the light of the existing experiments15,16 and found that the electrostatic interaction, reproduced through a scaling factor applied to the van der Waals UAHF radius, conforms a potential well at low H2O content in the system, most likely ascribable to the formation of hydrogen bonding when the amount of OH functionals in the melt are in sufficient amount to contour the discrete H2O molecules in solution. We then recast the results of the computation in terms of a Redlich-Kister sub-regular expansion to render them useful practically. The model has wide potential application, because it may be easily generalized to compositionally complex systems. For this purpose, it will be sufficient to assess the standard state properties of molecular H2O in the various pure liquid oxides that makes up the composition of complex melts (i.e., Al2O3, CaO, MgO, Na2O, K2O, etc.) and to find a conceptually appropriate expression of the reactive constant in the complex melt, in the light of the simple model developed for discrete binary systems.22 ACKNOWLEDGMENTS

This work has been partly supported by the MIUR PRIN Project No. 2009B3SAFK (topology of phase diagrams and lines of descent). We gratefully thank IPGP for an invited professorship to G.O. during which part of this work was concluded. FIG. 8. H2O isopleths generated by Model (1) in Table IV, coupled with the speciation calculations of Silver and Stolper22 (their Eqs. (5.1) and (5.2)). Fugacity coefficients in the gas phase computed with the SUPERFLUID code.23

1G.

Ottonello and P. Richet, J. Chem. Phys. 140, 044506 (2014). H. Fowler and E. A. Guggenheim, Statistical Thermodynamics (Cambridge University Press, Cambridge, 1939).

2R.

054503-9

Ottonello, Richet, and Vetuschi Zuccolini

J. Chem. Phys. 142, 054503 (2015)

3R.

21M.

4H.

22L.

A. Pierotti, J. Phys. Chem. 67, 1840 (1963). Reiss, Adv. Chem. Phys. 9, 1 (1965). 5H. Reiss and S. W. Mayer, J. Chem. Phys. 34, 2001 (1961). 6H. Reiss, H. L. Frisch, and J. L. Lebowitz, J. Chem. Phys. 31, 369 (1959). 7H. Reiss, H. L. Frisch, E. Helfand, and J. L. Lebowitz, J. Chem. Phys. 32, 119 (1960). 8R. A. Pierotti, J. Phys. Chem. 69, 281 (1965). 9R. A. Pierotti, Chem. Rev. 76, 717 (1976). 10J. Tomasi and M. Persico, Chem. Rev. 94, 2027 (1994). 11J. Tomasi, B. Mennucci, and E. Cancès, J. Mol. Struct.: THEOCHEM 464, 211 (1999). 12F. M. Floris and J. Tomasi, J. Comput. Chem. 10, 616 (1986). 13F. M. Floris, J. Tomasi, and J. L. Pascual-Ahuir, J. Comput. Chem. 12, 784 (1991). 14J. Caillet and P. Claverie, Acta Crystallogr., Sect. B: Struct. Crystallogr. Cryst. Chem. 34, 3266 (1978). 15G. C. Kennedy, G. J. Wasserburg, H. C. Heard, and R. C. Newton, Am. J. Sci. 260, 501 (1962). 16A. L. Boettcher, Am. Mineral. 69, 823 (1984). 17P. Richet and G. Ottonello, Elements 6, 315 (2010). 18M. J. Frisch, G. W. Trucks, H. B. Schlegel et al., Gaussian 03 Revision B.05, Gaussian, Inc., Wallingford, CT, 2004. 19A. E. Stearn and H. Eyring, J. Chem. Phys. 5, 113 (1937). 20Y. Bottinga, P. Richet, and D. F. Weill, Bull. Mineral. 106, 129 (1983).

L. Rivers and I. S. E. Carmichael, J. Geophys. Res. 92, 9247 (1987). Silver and E. Stolper, J. Geol. 93, 161 (1985). 23A. B. Belonoshko, P. Shi, and S. K. Saxena, Comput. Geosci. 18, 1267 (1992). 24M. Cossi, V. Barone, R. Cammi, and J. Tomasi, Chem. Phys. Lett. 255, 327 (1996). 25M. Cossi and V. Barone, J. Chem. Phys. 112, 2427 (2000). 26M. Valiev, E. J. Bylaska, N. Govind, K. Kowalski, T. P. Straatsma, H. J. J. van Dam, D. Wang, J. Nieplocha, E. Apra, T. L. Windus, and W. A. de Jong, Comput. Phys. Commun. 181, 1477 (2010). 27M. Swart, M. Solà, and F. M. Bickelhaupt, J. Chem. Phys. 131, 094103 (2009). 28S. Grimme, J. Comput. Chem. 25, 1463 (2004). 29S. Grimme, J. Comput. Chem. 27, 1787 (2006). 30Q. Wu and W. Yang, J. Chem. Phys. 116, 515 (2002). 31U. Zimmerli, M. Parrinello, and P. Koumoutsakos, J. Chem. Phys. 120, 2693 (2004). 32B. Civalleri, E. Garrone, and P. Ugliengo, Chem. Phys. Lett. 294, 103 (1998). 33K. Wendker, J. Thar, S. Zahn, and B. Kirchner, J. Phys. Chem. A 114, 9529 (2010). 34P. Richet and A. Polian, Science 281, 396 (1998). 35P. Richet, A. Whittington, F. Holtz, H. Behrens, S. Ohlhorst, and M. Wilke, Contrib. Mineral. Petrol. 138, 337 (2000). 36H. S. Frank and M. W. Evans, J. Chem. Phys. 13, 507 (1945).

Journal of Chemical Physics is copyrighted by AIP Publishing LLC (AIP). Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. For more information, see http://publishing.aip.org/authors/rights-and-permissions.