Ab initio calculations for equilibrium fractionations in multiple sulfur isotope systems

Available online at www.sciencedirect.com

Chemical Geology 249 (2008) 357 – 376 www.elsevier.com/locate/chemgeo

Ab initio calculations for equilibrium fractionations in multiple sulfur isotope systems Tsubasa Otake a,⁎, Antonio C. Lasaga b , Hiroshi Ohmoto a a

NASA Astrobiology Institute and Department of Geosciences, The Pennsylvania State University, University Park, PA 16802, USA b Geokinetics, State College, PA 16803, USA Received 18 May 2007; received in revised form 6 December 2007; accepted 20 January 2008 Editor: S. L. Goldstein

Abstract Ab initio calculations of vibrational frequencies (ν), reduced partition function ratios (β) of various gaseous and aqueous sulfur-bearing compounds, and isotopic fractionation factors (α) between them were performed on the multiple isotope systems using the Hartree-Fock (HF) and the Becke three-parameter Lee-Yang-Parr (B3LYP) density functional methods. Based on the comparison of the geometry, vibrational frequencies, and fractionation factors with those determined experimentally, the HF/6-31G(d) was chosen as the optimal method for the sulfur isotope systems. The results of computation for multiple sulfur fractionation factors indicate that the mass-dependent relationships, i.e., (33α − 1) / (34α − 1) and (36α − 1) / (34α − 1) values, converge only at high temperature (e.g., T N 500 °C) to 0.515 and 1.89, respectively, which are the values previously derived from the Bigeleisen-Mayer approximations. However, they significantly deviate from these values with decreasing temperature, and depending on the sulfur species pairs that are involved in equilibrium isotope exchange reactions. The (33α − 1) / (34α − 1) values and the (36α − 1) / (34α − 1) values range from 0.505 to 0.517 and from 1.88 to 1.96, respectively, at 0 °C. An isotopic crossover relationship, in which the fractionation factor changes from N 1.0 to b 1.0 (or from b1.0 to N 1.0) about a specific temperature (i.e., crossover temperature), was observed in several species couples (e.g., S2–H2S, S8–H2S, CS2–H2S, SO3–SO2− 4 (aq)). Since the crossover temperatures vary among different isotope molecules, large changes in the (33α − 1) / (34α − 1) and (36α − 1) / (34α − 1) values (e.g., N 1 to b 0) also occur near the crossover temperatures. The effects of solvation on the fractionation factors were investigated using the PCM models. The dielectric constant (ε) of solvent, which is the dominant solvent characteristic, affects the fractionation factors between gaseous and aqueous species (e.g., SO2(g) and SO2(aq)), but the effect levels out for ε N 20. Among the solvent models, the IEF-PCM model yields the greatest fractionation factors between gaseous and aqueous species. The effects of anharmonicities were also investigated for H2S and HS−. When anharmonicity is corrected without scaling, the DFT methods, rather than the HF methods, yield frequency values closer to those determined spectroscopically. Our results show that anharmonicity can significantly affect the reduced partition function ratios (0.5‰ for the 1000(34βH2S − 1) value and 0.1‰ for the 1000(34βHS− − 1) value) and fractionation factors; however, it has little effect on the multiple isotope fractionation factor ratios. Overall, temperature, harmonic frequency, and mass of atoms bonded to sulfur (M) are the three important parameters affecting the reduced partition function ratios, fractionation factors, and mass-dependent relationships. By decreasing temperature and increasing frequency and mass, the 1000(34β − 1) values increase and the (33β − 1) / (34β − 1) values decrease (i.e., deviate from “expected” 0.515). The (33β − 1) / (34β − 1) values converge to 0.516 only at higher temperature (e.g., T N 600 °C), lower frequency (ωe b 500 cm− 1), and smaller mass (M b 7 amu) and become as low as 0.503 at T = 50 °C, ωe = 2500 cm− 1, and M = 56 amu. © 2008 Elsevier B.V. All rights reserved. Keywords: Ab initio calculation; Sulfur isotopes; Equilibrium; Mass-dependent fractionation; Solvent model; Anharmonicity

⁎ Corresponding author. Department of Geosciences, The Pennsylvania State University, 437 Deike Bldg. University Park, PA 16802, USA. Tel.: +1 814 360 1225; fax: +1 814 863 2001. E-mail address: [email protected] (T. Otake). 0009-2541/$ - see front matter © 2008 Elsevier B.V. All rights reserved. doi:10.1016/j.chemgeo.2008.01.020

358

T. Otake et al. / Chemical Geology 249 (2008) 357–376

where 34S is substituted with 33S and 36S for δ33Si and δ36Si, respectively. In the early stages of the development of sulfur isotope geochemistry (1960s and 1970s), the δ33S and δ36S values of some natural samples were measured. However, they were mostly limited to extraterrestrial samples, such as meteorites and lunar rocks (e.g., Hulston and Thode, 1965; Thode and Rees, 1971). Very few studies were carried out to determine δ33S and δ36S values of terrestrial samples before the 1990s (e.g., Hoering, 1989). This was because precise determination of these values was difficult due to the low abundance of 33S and 36S, and because it was generally believed that measurements of minor isotopes of terrestrial samples would not give additional information since the theory of massdependent isotope fractionation (Bigeleisen and Mayer, 1947) predicted the following linear relationships among the δ values:

Farquhar et al. (2001) showed that SO2 photolysis with wavelength of 193 nm radiation produced S0 and SO42− with strong MIF-S signatures, and suggested a linkage between the presence/absence of MIF-S signatures in sedimentary rocks and the evolution of atmospheric oxygen. All the discussions presented in recent papers dealing with multiple sulfur isotope data of natural samples have been based on the acceptance of the Bigeleisen-Mayer approximations (i.e., Eqs. (2) and (3)) as precise expressions for all equilibrium isotope exchange reactions. This has led some recent researchers (e.g., Johnston et al., 2005; Jamieson et al., 2006; Ono et al., 2006) to use very small deviations (b 0.1‰) from the approximations to distinguish kinetic processes (e.g., biological reactions) from equilibrium processes. However, the validity of these approximations has not been thoroughly evaluated either experimentally or theoretically. Previous experimental studies have investigated 32S–34S equilibrium isotope exchange reactions between various sulfurbearing species (see Ohmoto and Goldhaber, 1997 for a compilation), but no experimental study has focused on other sulfur isotope systems (32S–33S, 32S–36S). The main purpose of this study is, therefore, to evaluate the mass-dependent approximations for various sulfur-bearing species using ab initio calculations to estimate the vibrational frequencies of four different sulfur isotope molecules. These frequencies are required for the calculation of isotope fractionation factors and described in detail in the following section. Recently, ab initio calculations have been applied to predict equilibrium fractionation factors in other isotope systems, including H, O (Driesner et al., 2000), Li (Yamaji et al., 2001), B (Oi, 2000; Oi and Yanase, 2001; Liu and Tossell, 2005; Zeebe, 2005), Cr (Schauble et al., 2004; Ottonello and Zuccolini, 2005), Fe (Anbar et al., 2005), Cu (Seo et al., 2007), and Mo (Tossell, 2005).

d33 Si ¼ 0:515  d34 Si

2. Computational methods

1. Introduction Sulfur has four stable isotopes, 32S, 33S, 34S, and 36S with an average abundance of 95.04, 0.75, 4.20, and 0.01%, respectively (Ding et al., 2001). Until recently, applications of sulfur isotopes in studies of geological processes (e.g., Ohmoto and Goldhaber, 1997; Rye, 2005; Riccardi et al., 2006), biological processes (e.g., Canfield and Thamdrup, 1994; Detmers et al., 2001), and environmental problems (e.g., Edraki et al., 2005) were primarily based on interpretations of the ratio of the two most abundant isotopes, 34S/32S, of a sample. The sulfur isotopic composition of a compound, i, is usually expressed as a δ34Si value, which is defined as a per mil deviation of the 34 32 S/ S ratio of the compound relative to that of an international standard (CDT, Canyon Diablo Troilite, or V-CDT, a virtual standard adopted by the IAEA (Vienna)): 
34 32
  ð1Þ d34 Si ¼ S= SÞi = 34 S=32 SÞstd  1  1000

ð2Þ

2.1. Calculation of fractionation factors

and d36 Si ¼ 1:89  d34 Si

ð3Þ

The recent surge of interest in determining the δ33S and δ36S values of terrestrial samples has stemmed from the finding by Farquhar et al. (2000) that some sulfide and sulfate minerals in sedimentary rocks older than ~ 2.0 Ga do not fall on the massdependent fractionation lines (i.e., Eqs. (2) and (3)); this abnormal fractionation is termed “mass-independently-fractionated sulfur isotopes (MIF-S)”. Deviations in δ33S and δ36S values from the mass-dependent lines are expressed respectively by Δ33S and Δ36S values, which are defined as: D33 Si ¼ d33 Si  0:515  d34 Si

ð4Þ

and D36 Si ¼ d36 Si  1:89  d34 Si

ð5Þ

The isotope exchange fractionation factor between two sulfur-bearing compounds, A and B, iαA–B is defined as: i



 aA−B ¼ ði S=32 SÞA =ði S =32 SÞB ¼ di SA þ 1000 = di SB þ 1000 ð6Þ

where (iS/32S)A and (i S/32S)B are the abundance ratio of sulfur isotopes iS and 32S in compounds A and B, respectively. Eqs. (2) and (3) were derived from an approximation (see also Section 6) by Bigeleisen and Mayer (1947): i



 aA−B i1 þ ðCA  CB Þ i m  32 mÞ= i m 32 m

ð7Þ

where im and 32m are masses of isotopes iS and 32S; and CA and CB are specific constants for A and B, respectively, but are independent of the isotope mass. Therefore, the following

T. Otake et al. / Chemical Geology 249 (2008) 357–376

359

mass-dependent relationships are expected among 33α, 34α and 36 α:

harmonic vibrational frequency of the jth vibrational mode in s− 1:

ð33 a  1Þ=ð34 a  1Þið1=ð32  33ÞÞ=ð2=ð32  34ÞÞ ¼ 0:515

mj ¼ x j c

ð8Þ and ð36 a  1Þ=ð34 a  1Þið4=ð32  36ÞÞ=ð2=ð32  34ÞÞ ¼ 1:89 ð9Þ Since α values are generally very close to 1, the following approximations are often used: ða  1Þilna and


i

ð10Þ


a  1 = 32 a  1Þilni a=ln32 a

ð100 Þ

The theoretical calculation of isotope fractionation factors in equilibrium systems was developed by Urey (1947) and Bigeleisen and Mayer (1947). Subsequently, it was extensively treated by Richet et al. (1977) and more recently reviewed by Chacko et al. (2001) and Schauble (2004). Equilibrium isotope fractionation is a quantum mechanical phenomenon usually arising from a difference in the vibrational energies between compounds that contain different isotopes. Suppose the following isotope exchange reaction between H2S and SO42− is considered:

ð14Þ

where ωj is the harmonic frequency in cm− 1 and c is the speed of light. One of the earliest theoretical estimations of sulfur isotope fractionation was performed by Thode's group (Tudge and Thode, 1950; Thode et al., 1971). They calculated 34α values between simple gaseous sulfur compounds using available frequency data for 32S compounds and an empirical force field model to compute frequencies of the 34S-substituted compounds. Their works were later reviewed by Richet et al. (1977), particularly for the SO2–H2S isotope exchange reaction. Tudge and Thode (1950) also calculated an 36α value for the SO42−–H2S isotope exchange reaction and demonstrated that the (36α − 1) value is approximately twice the (34α − 1) value (1.96 at 25 °C). Sakai (1968) and Oi et al. (1985) expanded the calculation for 34 α values to higher temperatures and more varied sulfur-bearing gaseous compounds and minerals. More recently, Farquhar and his colleagues (Farquhar et al., 2003; Farquhar and Wing, 2003; Ono et al., 2007) calculated both 33α and 34α values for some sulfur compound pairs (e.g., SO42−−H2S, SO32−−H2S) using an empirical force field model, and confirmed the mass-dependent relationships during sulfur isotope exchange reactions between those species. 2.2. Ab initio methods

H2 i S þ

32

32 i 2 SO2 4 ¼ H2 S þ SO4

where i is 33, 34 or 36, the equilibrium constant K, also defined as α here, can be obtained from the reduced partition function ratios (RPFRs), iβ: K ¼ i aSO2− ¼ 4 −H2 S

i i S S = ¼ i bSO2− =i bH2 S 32 32 4 S SO2− S H2 S

ð11Þ

4

The RPFRs of a compound for sulfur masses of 32 and i is the ratio of the quantum partition function and the classical vibrational partition function, which, under the Born-Oppenheimer and harmonic oscillator approximations (Bigeleisen and Mayer, 1947), can be written as:


Uj exp i U j =2 1  exp 32 U j


b ¼ j 32 j¼1 Uj exp 32 U j =2 1  exp i Uj

ð12Þ

where Uj ¼ hmj =kT

ð13Þ

i

f

i

f is the number of degrees of freedom for the vibrational mode (f = 3N − 5 for linear molecules and 3N − 6 for other molecules, where N is number of atoms), h is Plank's constant, k is Boltzmann's constant, T is temperature in Kelvin, and νj is the

We calculated vibrational frequencies of sulfur isotopic molecules with the Gaussian 03 program package (Frisch et al., 2004) using both the Hartree-Fock (HF) method and the Becke three-parameter Lee-Yang-Parr (B3LYP) method, which is a gradient-corrected hybrid HF-density functional theory (DFT) (Lee et al., 1988; Becke, 1993). An HF calculation considers the electron–electron repulsion as an average effect in the Hamiltonian; each electron interacts with an averaged electron density from the other electrons. In contrast, DFT calculations, while not strictly first-principles, include some of the electron correlation accounting for the instantaneous interactions of pairs of electrons (Foresman and Frisch, 1996; Jensen, 1999). We chose two types of extended basis sets for both HF and DFT calculations: 6-31G(d) and 6-311G(d,p). Although both basis sets use 6 Gaussian functions for inner core atomic orbitals and split valence and polarization functions, the 6-311G(d,p) basis has additional split valence and polarization functions of d- and p-type (Foresman and Frisch, 1996). Vibrational harmonic frequencies calculated by first-principle methods tend to be systematically higher than experimentally-derived fundamental frequencies (see Section 3.2). This is, partly, because the fundamental frequencies obtained from vibrational spectra are affected by anharmonicities, which tend to decrease the frequencies from the true harmonic frequencies. In addition, there are intrinsic errors from the limitations of the basis set and electron correlation corrections. In order to correct for anharmonic effects and to be in overall good agreement with

360

T. Otake et al. / Chemical Geology 249 (2008) 357–376

Table 1 Comparison of geometry (bond distance (Å) and angle (°)) of some gaseous and aqueous sulfur molecules between the values determined experimentally and those calculated in this study using various ab initio methods Experimental (gaseous) HS− H2S S2 S8 CS2 SO2 SO3 SO2− 4

dSH dSH aHSH dSS dSS aSSS dSC dSO aOSO dSO aOSO dSO aOSO

1.328 c 92.2 c 1.889 d 2.059 e 107.9 e 1.554 c 1.432 c 119.5 c 1.418 f 120 f 1.486 g

1σ h a b c d e f g h

HF/6-31G(d)

HF/6-311G(d,p) PCM

1.342 1.326 b 94.37 1.878 2.059 107.62 1.545 1.414 118.82 1.405 120.00 1.487 109.47 0.011

a

1.340 1.330 95.11 1.877

1.545 1.419 116.76 1.405 120.00 1.479 109.47

B3LYP/6-31G(d)

PCM 1.340 1.332 99.20 1.877 2.069 107.43 1.543 1.407 118.73 1.396 120.00 1.481 109.47 0.016

1.344 1.337 94.93 1.875

1.544 1.412 116.66 1.397 120.00 1.473 109.47

B3LYP/6-311G(d,p)

PCM 1.363 1.350 92.76 1.929 2.100 109.20 1.563 1.464 119.12 1.453 120.00 1.525 109.47 0.032

1.362 1.354 93.48 1.928

1.564 1.467 117.43 1.452 120.00 1.516 109.47

PCM 1.355 1.348 92.56 1.928 2.110 109.06 1.560 1.457 118.81 1.446 120.00 1.521 109.47 0.031

1.361 1.355 93.22 1.926

1.561 1.460 117.06 1.512 120.00 1.512 109.47

MP2/631G(d) 1.353 1.340 93.36 1.939 2.067 108.48 1.563 1.478 119.82 1.460 120.00 1.521 109.47 0.033

The IEF-PCM model was used for aqueous species. Bold numbers indicate values closest to the experimental values among those calculated using different methods for gaseous species. Herzberg (1966). Huber et al. (1979). Cotton and Wilkinson (1972). Meyer et al. (1991). McGinnety (1972). Standard deviation of the bond distance compared to the experimental data for gaseous species was calculated using the following equation:

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
2
2
2
2
2
2 dSHðH2 SÞ  dTSHðH2 SÞ þ dSSðS2 Þ  dTSSðS2 Þ þ dSSðS8 Þ  dTSSðS8 Þ þ dSCðCS2 Þ  dTSCðCS2 Þ þ dSOðSO2 Þ  dTSOðSO2 Þ þ dSOðSO3 Þ  dTSOðSO3 Þ 1r ¼ 6 where ⁎ indicates its experimental value.

experimentally-derived frequencies, the calculated HF frequencies are often multiplied by an appropriate scaling factor for each basis set. The following scaling factors were used in our calculations: 0.8929 for HF/6-31G(d) (Foresman and Frisch, 1996), 0.9085 for HF/6-311G(d,p) (CCCBDB, Computational Chemistry Comparison and Benchmark Database — http:// srdata.nist.gov/cccbdb), 0.9613 for B3LYP/6-31G(d) (Foresman and Frisch, 1996), and 0.9668 for B3LYP/6-311G(d,p) (CCCBDB). With these scaling factors, the calculations produce high precision values for the fundamental frequencies. For each combination of method and basis set, we calculated the optimum geometric configurations (bond distances and angles), vibrational frequencies, and RPFRs for all four stable isotopes of sulfur in simple compounds (HS−, H2S, S2, S8, CS2, SO2, SO3, and SO42−) in the temperature range from − 70° to 650 °C. (For molecules containing multiple numbers of S atoms (i.e., S2, S8 and CS2), only one 32S atom was substituted by a 33 S, 34S, or 36S atom.) We determined the optimal ab initio method to estimate 33α, 34α and 36α values for sulfur isotope exchange reactions by comparing our calculated geometry, vibrational frequencies and fractionation factors (34α) with the experimental data. Those 33α, 34α and 36α values calculated using the optimal method were used to evaluate the validity of mass-dependent approximations. Furthermore, we evaluated solvent and anharmonic effects on the RPFRs, fractionation factors, and mass-dependent approximations. The detailed methods are described in each

section (Section 4 for solvent effects and Section 5 for anharmonic effects). 3. Results of ab initio calculations based on harmonic oscillation model 3.1. Geometry of sulfur-bearing compounds Table 1 compares the predicted geometric configurations of sulfur-bearing compounds, obtained using different ab initio methods, with those determined experimentally (Herzberg, 1966; Cotton and Wilkinson, 1972; McGinnety, 1972; Huber and Herzberg, 1979; Meyer et al., 1991). It shows that the B3LYP methods tend to overestimate the bond lengths. The HF methods tend to underestimate them, but the calculated values generally deviate less from the experimental data than calculated values with the B3LYP methods. However, the B3LYP methods give better predictions for the bond angles for H2S and SO2. Table 1 includes results from HF methods as well as methods that also add some electron correlation such as MP2 or the density functional B3LYP. One would think that adding electron correlation would improve the accuracy of the prediction. However, as has been found in numerous other cases in the quantum chemical literature, cancellation of errors (including anharmonic corrections, in this case — see Section 5) can often lead a “lower-level” method to perform better than a “higher-level” method. For the sulfur molecules studied here,

T. Otake et al. / Chemical Geology 249 (2008) 357–376

361

Table 2 Comparison of vibrational frequencies of some gaseous and aqueous sulfur (32S) molecules calculated in this study using various ab initio methods with those determined spectroscopically Method

Frequencies (cm− 1)

1σe

H2S Experimental HF/6-31G(d) HF/6-311G(d,p) B3LYP/6-31G(d) B3LYP/6-31G(d,p)

S2 a

1183 1221 1199 1202 1168

a

2615 2607d 2602 2588 2606

a

2626 2618 2610 2608 2591

CS2 b

720 726 721 670 657

a

397 398 388 389 379

SO2 a

658 647 660 648 652

a

1535 1414 1439 1497 1503

a

518 529 546 483 492

SO2− 4 (aq)

SO3 a

1151 1213 1236 1096 1100

a

1362 1401 1421 1285 1280

a

498 502 513 437 445

a

530 521 532 478 481

a

1065 1083 1106 983 989

a

1391 1388 1406 1299 1298

452c 422 435 399 406

617c 582 595 550 553

982c 950 968 900 900

1105c 1084 1084 1053 1029

34 33 56 59

⁎Each molecule with N atoms has 3N-6 vibrational mode (3N-5 for linear molecules). Shimanouchi (1972). b Huber and Herzberg (1979). c Hester and Plane (1963). d Bold numbers indicate values closest to the experimental values among those calculated using different methods. e Standard deviation of the frequencies to the experimental data for 5 gaseous species was calculated using the following equation: a

1r ¼

vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi uP 1
4
3
6
9 2 P 2 P 2 P 2 P 2 P u3
u m  mTjðH2 SÞ þ mjðS2 Þ  mTjðS2 Þ þ mjðCS2 Þ  mTjðCS2 Þ þ mjðSO2 Þ  mTjðSO2 Þ þ mjðSO3 Þ  mTjðSO3 Þ þ ðmjðSO2− Þ  mTjðSO2− Þ Þ2 4ðaqÞ 4ðaqÞ t j jðH2 SÞ j j j j j 26

where ⁎ indicate its experimental values.

indeed, the HF methods performed slightly better than the electron correlation methods for predicting molecular structure and vibrational properties. Note, however, that electron correlation methods would be essential if one were interested in calculating accurate changes in energies involving bond formation or breaking (e.g., calculating dissociation energies, ΔEor, or activation energies, ΔE‡). In addition, increasing the basis set beyond the 6-31G(d) level does not appreciably improve the results. This observation is also in accord with other studies (Hehre et al., 1986; Scott and Radom, 1996). The geometric configurations of the sulfur-bearing compounds in water, which are computed from the IEF-PCM models, differ only slightly from those computed for simple gaseous species (Table 1). For example, the lengths of S–S bonding in a S2 molecule and S–O bonding in a SO42− molecule became shorter by 0.002 Å and 0.01 Å, respectively, while those of other molecules became longer when the IEF-PCM model was applied. However, solvent effects are much more important in the vibrational frequency calculations.

popular in theoretical studies for vibrational frequencies and isotope fractionation (e.g., Anbar et al., 2005; Ottonello and Zuccolini, 2005), it has been noted that higher-levels of theory do not always guarantee better results in frequency calculations (Liu and Tossell, 2005). Note, however, that it is the scaled Hartree-Fock frequencies which outperform the scaled DFT frequencies. The scaling factors for the DFT methods employed here (0.9613 and 0.9668) are closer to unity than the HF methods (0.8929 and 0.9085), in accord with the improvement one would expect from the addition of electron correlation. It is in the combination of harmonic calculations and scaling that the HF/6-31G(d) method comes up effectively better. Furthermore, the agreement is between the scaled ab initio results and the fundamental (i.e., experimental) frequencies. If anharmonic effects are small, these scaling factors are precisely what are needed to also accurately calculate partition functions. However, if anharmonic effects are not negligible, they will affect the values for the fundamental frequencies as well as the partition functions (and hence the fractionation factors; see Section 5).

3.2. Vibrational frequencies of sulfur-bearing compounds Table 2 compares the scaled harmonic vibrational frequencies calculated for gaseous 32S compounds with those obtained experimentally from spectroscopic measurements (Hester and Plane, 1963; Shimanouchi, 1972; Huber and Herzberg, 1979). Unscaled vibrational frequencies of all four sulfur isotope compounds calculated using the HF/6-31G(d) method are given in the Supplementary material (Table S1). Table 2 shows that both HF and B3LYP methods are in reasonable agreement with the experimental data (within 3–6% of uncertainty as 1σ). As in the geometry calculations, frequencies calculated using HF methods are generally in better agreement with the experimental data than those calculated using B3LYP methods (~ 35 versus ~ 50 cm− 1 as 1σ). Although B3LYP methods are becoming

3.3. Reduced partition function ratios and equilibrium fractionation factors for 34S–32S exchange reactions The RPFRs between 32S- and 34S-species (i.e., 34β in Eq. (11)) of H2S, HS−, S2, S8, CS2, SO2, SO3, SO42−, and their aqueous species (i.e., the IEF-PCM modeled species) have been computed over the temperature range of −70 °C to 650 °C using Eqs. (12) and (13). The 34β values calculated with each ab initio method are presented in the Supplementary material (Tables S2, S3, S4, and S5). Equilibrium fractionation factors for 34S–32S exchange reactions between a sulfur species (i) and H2S(g) (i.e., 34 aiH2 S ) were, then, computed using Eq. (11); they are compared in Fig. 1 with those estimated by other theoretical studies (Sakai, 1968; Ono et al., 2007), as well as those based on available

362

T. Otake et al. / Chemical Geology 249 (2008) 357–376

− 34 Fig. 1. Comparison of the fractionation factors for 34S–32S isotope exchange reactions between H2S and (a) SO2− 4 (aq), (b) SO2, (c) HS (aq), and (d) S8 and H2S ( αi−H2S values) among: the values calculated in this study using the ab initio methods (HF/6-31G(d), HF/6-311G(d,p), B3LYP/6-31G(d), and B3LY/6-311G(d,p)); those estimated from an empirical force field model (Sakai, 1968; Ono et al., 2007); and those estimated from experimental data (Ohmoto and Goldhaber, 1997).

experimental data (Ohmoto and Lasaga, 1982; Ohmoto and Goldhaber, 1997) The experimentally-determined fractionation factors are typically limited to high temperature ranges (N 200 °C) because of the difficulty in attaining equilibrium between sulfur species with different valence states at low temperatures (Ohmoto and Lasaga, 1982). The assumption of no isotope fractionation between H2S(g) and H2S(aq), which was adopted by Ohmoto and Goldhaber (1997) during the evaluation of experimental data on SO2–H2S and SO42−−H2S systems, is justifiable, because our calculation predicts a ~0.2‰ fractionation between H2S(g) and H2S(aq) at 25 °C (see Section 4). As expected from the comparison of the geometry and vibrational frequencies, 34αi − H2S values calculated using the B3LYP methods, compared to those calculated by the HF methods, generally show larger deviations from the experimentallydetermined 34 aiH2 S values, by as much as ~ 5‰ in the 1000  (α − 1) values. The 1000 34 aSO H S  1 values calculated using the two HF methods (HF/6-31G(d) and HF/6-311G(d,p)) with the IEF-PCM model are essentially identical with each other (b ± 0.1‰ at N 25 °C), and they agree very closely (± 0.5‰) with the experimentally-determined values (Fig. 1a).
 The two HF methods also gave the 1000 34 aS8 H2 S  1 values that were within 0.7‰ of the experimentally determined values 2 4ðaqÞ

2






(Fig.  1d). However, the calculated 1000 34 aSO H S  1 and 34 1000 aHS H S  1 values are smaller than the experimentally-determined values, by as much as ~ 2‰ at T = 250 °C. These discrepancies will be further discussed in Sections 4 and 5. Based on the results of calculations on the geometry and vibrational frequencies of various sulfur-bearing compounds, and fractionation factors (34α) among them, we recommend the scaled HF/6-31G(d) as the optimal method for studies on chemical and isotopic exchange reactions among sulfur-bearing compounds. Several other recent theoretical investigations have also settled on the HF/6-31G(d) method as optimal. East and Radom (1997), in the formulation of the E1 procedure for accurate calculations of thermodynamic properties of gaseous molecules, recommend the MP2/6-31G(d) method for geometries and the scaled HF/6-31G(d) method for vibrational frequencies. Similarly, accurate methods for calculations of reaction energies, ΔE0, such as the successful G3 method, use the scaled HF/6-31G(d) for vibrational frequencies (to get the zero-point vibrational contribution to the ΔE) (Curtiss et al., 1998). Even the calculations of molecular electrostatic potentials, MEP's, are often carried out with the HF/6-31G(d) (Levine, 2000). 2

 ðaqÞ

2

2

T. Otake et al. / Chemical Geology 249 (2008) 357–376

363

Table 3 Calculated regression coefficients for reduced partition function ratios (33β, 34β, and 36β) of sulfur-bearing compounds as a function of temperature (−70 to 650 °C) calculated using the HF/6-31G(d) method 1000ðb  1Þ ¼

33

H2S HS− S2 S8 CS2 SO2 SO3 SO2− 4 H2S(aq) HS−(aq) S2(aq) CS2(aq) SO2(aq) SO3(aq) SO2− 4 (aq)

β

A B C D  1011 þ 3  108 þ 2  106 þ  103 ðT in KÞ T4 T T T β

β

34

36

A

B

C

D

A

B

C

D

A

B

C

D

0.151 0.065 0.012 − 0.013 0.061 0.422 0.635 0.425 0.265 0.121 − 0.011 0.031 0.409 0.662 0.303

− 1.934 − 0.817 − 0.632 − 0.175 − 1.293 − 6.502 − 9.768 − 8.028 − 2.795 − 1.246 − 0.458 − 0.985 − 6.193 − 9.698 − 6.788

0.903 0.366 0.801 0.747 1.186 3.951 6.145 6.161 1.106 0.471 0.759 1.073 3.766 5.951 5.638

0.309 0.129 − 0.067 − 0.022 − 0.073 − 0.272 − 0.375 − 0.606 0.136 0.055 − 0.037 − 0.062 − 0.337 − 0.501 − 0.457

0.292 0.126 0.021 − 0.025 0.115 0.798 1.187 0.762 0.513 0.234 − 0.022 0.058 0.766 1.227 0.531

−3.745 −1.584 −1.192 −0.316 −2.456 −12.329 −18.253 −14.563 −5.409 −2.415 −0.861 −1.866 − 11.683 −18.041 −12.24

1.753 0.710 1.547 1.444 2.292 7.645 11.869 11.820 2.146 0.914 1.467 2.073 7.268 11.462 10.812

0.597 0.250 − 0.127 − 0.039 − 0.137 − 0.529 − 0.727 − 1.122 0.262 0.107 − 0.070 − 0.118 − 0.643 − 0.944 − 0.832

0.550 0.239 0.034 − 0.043 0.209 1.441 2.107 1.255 0.967 0.443 − 0.045 0.104 1.358 2.139 0.842

− 7.062 − 2.993 − 2.138 − 0.520 − 4.468 − 22.318 − 32.076 − 24.018 − 10.198 − 4.562 − 1.533 − 3.380 7.904 − 31.457 − 19.999

3.323 1.345 2.905 2.716 4.310 14.395 22.267 21.899 4.061 1.730 2.759 3.900 13.626 21.400 20.033

1.121 0.472 − 0.230 − 0.064 − 0.248 − 1.005 − 1.359 − 1.919 0.490 0.201 − 0.126 − 0.213 − 1.173 − 1.687 − 1.383

Typical misfit of the regression lines is b0.05‰ and the maximum misfit is 0.13‰ (36βSO3 at 223 K).

3.4. Reduced partition function ratios and fractionation factors for multiple sulfur isotope systems We have also computed the RPFRs (33β and 36β) and isotope fractionation factors (33α and 36α) for 33S–32S and 36S–32S exchange reactions. Here we will discuss only the results calculated using the HF/6-31G(d) method for the reasons presented in the previous sections. The calculated 33β and 36β values are presented in the Supplementary material (Tables S6 and S7). Table 3 presents the regression lines for the temperature dependence of 33β, 34β, and 36β for the 15 sulfurbearing species. Note that, if the frequencies and the normal modes of vibration for 32S molecules are in good agreement with the experimental data, then those for other sulfur isotope molecules are guaranteed to be accurate. This is because any shift in frequencies due to the mass occurs only from changes in the mass term in front of the potential matrix term: 1 A2 V pffiffiffiffiffiffiffiffiffiffi mi mj Axi Ayj

ð15Þ

where V is the potential energy, and mi and mj are the masses of the atom with coordinate xi and yj, respectively. If we know the ∂2V / ∂xi∂yj terms accurately, the mass variation will follow from Eq. (15) just as accurately for any mass, m. The fact that the ab initio fractionation factors are better than those estimated by previous studies that used the experimental vibrational frequencies and empirical force field models, indicates the importance of knowing directly ∂2V / ∂xi∂yj, which is not the same as using an assumed empirical force field. The calculated fractionation factors between H2S and other S species for the multiple isotope systems are compared in Fig. 2.

Most of the fractionation factors are not linear functions of 1/T2 because of the high vibrational frequencies of H-bearing species within the harmonic approximation (see Section 6). Results of our calculations agree with the general rules of equilibrium isotope fractionation describing that: (i) molecules involving stiffer bonding concentrate the heavier isotopes; and (ii) the magnitudes of fractionations for 33S–32S exchange are roughly one half of those for 34S–32S exchange, and one quarter of those for 36S–32S exchange at a given temperature, as suggested by the mass-dependent approximations of Eqs. (8) and (9). However, an important finding from our study is that these approximations are not always valid. The (33α − 1)/(34α − 1) and (36α − 1) / (34α − 1) values converge to 0.515 and 1.89, respectively, only at high temperatures (e.g., T N 500 °C), but they significantly deviate from these values with decreasing temperature: at 0 °C the ranges are 0.505–0.517 for (33α − 1)/ (34α − 1) and 1.88–1.96 for (36α − 1) / (34α − 1) (Fig. 3). Note that the ranges for lniα/ln34α values are much smaller than those for (iα − 1)/ (34α − 1) values, simply because the logarithm function masks the actual variations in the ratios through the approximation (10). Since iαA−B is related to the δiSA and δiSB values by Eq. (6), the (33αA−B − 1) / (34αA−B − 1) (or (36αA−B − 1) / (34αA−B − 1)) value corresponds to the slope of a line connecting the δS values of A and B on a δ33S–δ34S (or δ36S–δ34S) plot. In other words, the mass-dependent line may have a slope ranging from 0.505 to 0.517 (rather than the unique 0.515) on a δ33S–δ34S diagram and from 1.88 to 1.96 (rather than the unique 1.89) on a δ36S–δ34S diagram. Thus, the mass-dependent “line” is really a massdependent “band”. Sulfur isotope exchange reaction between SO42− and H2S is not likely to attain equilibrium at temperatures below ~100 °C at pH N ~2 (Ohmoto and Lasaga, 1982) and variations in sulfur isotopic composition at near surface conditions are mostly dictated

364

T. Otake et al. / Chemical Geology 249 (2008) 357–376

1.9334 ± 0.0172, respectively. For comparison, the   36 aSO2− H2 SðaqÞ  1 = 34 aSO2− H2 SðaqÞ  1 value at 25 °C is 4ðaqÞ 4ðaqÞ 1.955 in our calculation, while Tudge and Thode (1950) yielded 1.96 for the same value. Farquhar et al. (2003) calculated   ln 33 aSO2−4 H2 S =ln 34 aSO2−4 H2 S in the temperature range of 0– 100 °C, and yielded 0.5148 ± 0.0001, which is almost identical with the value we obtained in the same temperature range (0.5149 ± 0.0001). These agreements suggest that the simple methods used by the previous investigators, at least for the SO42−H2S system, resulted in reasonable values based on rigorous ab initio calculations.

3.5. Crossover of isotope fractionation

Fig. 2. Sulfur isotope fractionation factors between H2S(g) and various sulfurbearing species (HS−(aq), H2S(aq), S2(g), S8(g), SO2(g), and SO2− 4 (aq)), calculated at a temperature range of 0–650 °C using the HF/6-31G(d) method: (a) 1000 (33α − 1), (b) 1000(34α − 1), and (c) 1000(36α − 1).

by kinetic isotope effects (Ohmoto and Goldhaber, 1997). However, the equilibrium fractionation factors between SO42− and H2S provide usefulreferences in interpretation of natural data.  The mean values for 33 aSO2−4ðaqÞ H2 SðaqÞ  1 = 34 aSO2−4ðaqÞ H2 SðaqÞ  1 





and 36 aSO2−4ðaqÞ H2 SðaqÞ  1 = 34 aSO2−4ðaqÞ H2 SðaqÞ  1 values over the temperature range of 0–300 °C are 0.5097 ± 0.0024 and

Anomalous fractionation factor ratios due to differences for the different isotopes in the crossover temperature, at which the (α − 1) value switches signs from positive to negative (or vice versa), were suggested by Deines (2003) for the S2–H2S system, as well as Oi et al. (1985) for other sulfur systems, such as SO2– SOCl2, SO3–SO2Cl2, and CS2–SPCl3. Although Deines' calculation using the Urey-type model gave the crossover temperature of ~ 72 °C in the S2–H2S system, our first-principle calculation resulted in the crossover temperature of ~ 110 °C (Fig. 4a). We have found the crossover relationship in other sulfur-bearing species couples, including S8–H2S, CS2–H2S, and SO3–SO42−(aq) (Fig. 4). The ratios of fractionation factors significantly deviate from the approximated values (0.515 for (34α − 1) / (33α − 1) and 1.89 for (36α − 1) / (33α − 1)) within a significant temperature range around the crossover temperatures, as well as at crossover temperatures where the ratios go from infinitely positive to infinitely negative values. For example, in the S2–H2S system, a few percent deviation from the approximated values may occur within a 10 °C region around 110 °C; the (33α − 1) / (34α − 1) values are 0.4948 (i.e., − 3.9% deviation from 0.515) at 100 °C and 0.5352 (+ 3.9% deviation) at 120 °C. In the SO3–SO42−(aq) system, the (33α − 1) / (34α − 1) values are 0.5568 (+ 8.1% deviation from 0.515) at 115 °C and 0.4792 (− 7% deviation) at 135 °C while the crossover occurs ~ 125 °C. Our study suggests that there are many other species couples with crossover relationships, particularly when the species couples have small fractionation factors and one of them has hydrogen atoms in the molecular structure; this was also suggested by Oi et al. (1985). However, our calculation indicates that Δ33S and Δ36S values arising from crossover relationships are small (b0.01‰) because fractionation factors themselves are very small (b ± 0.05‰ in the 1000(34α − 1) value) in spite of anomalously large deviations from 0.515 and 1.89 in the (iα − 1) / (34α − 1) values. In the case of S2–H2S system, the Δ33S values are almost constant at ~ − 0.002 in a temperature range of ± 10 °C about the crossover temperature. Therefore, a crossover relationship may not yield large Δ33S and Δ36S values unless, as Deines (2003) suggested, a Rayleigh process is involved in the isotope fractionation.

T. Otake et al. / Chemical Geology 249 (2008) 357–376

4. Evaluation of solvent effects on isotope fractionation factors and the mass-dependent relationships A solvent can have a major effect on equilibrium constants and rates of reaction (Lasaga, 1990, 1998). Solvent effects are commonly dealt in quantum mechanics by continuum solvent models, referred as Self-Consistent Reaction Field (SCRF) methods, in which the detailed molecular structure of the solvent is ignored, and the solvent is modeled as a continuous dielectric surrounding a cavity that contains the solute molecule. The method, then, adds a potential energy term, Vsolv, to the molecular electronic Hamiltonian to calculate the effect of the solvent. Vsolv is obtained by placing charges on various parts of the cavity surface and calculating the effect of the electrostatic field on both the polarization of the dielectric continuum and the electric moments of the molecule. The charges themselves depend on the solute electrons and nuclei in a self-consistent manner (in the SCF solvent methods). The self-consistency criteria are added to the usual Hartree-Fock self-consistency relating the average electron density to the electronic wavefunction. The key quantities defining this solvent continuum are the dielectric constant, ε (Robinson and Stokes, 1959), and cavity shape (Tomasi and Persico, 1994). For anisotropic solvents, the dielectric constant is really a second-rank tensor, but we did not pursue fractionation in these systems here. Different methods exist for quantum mechanical calculations of solvent effects, depending on models for the cavity

365

(Foresman and Frisch, 1996). For example, the Onsager method inserts the solute inside a spherical cavity of radius a (Robinson and Stokes, 1959). The PCM (Polarizable Continuum Model) model uses a sphere of radius 1.2 times the Van der Waal's radius around each atom of the molecule, and puts charges on the surface resulting from the intersecting spheres to simulate the external field of the solvent (Tomasi and Persico, 1994). The SCI-PCM model carries out the calculation in a self-consistent fashion, using the electron density of the solute itself (in any one iteration) to determine both the shape of the cavity and the external field potential. An integral equation formulation (IEFPCM) has been used to obtain good results for aqueous ionic solutions (Cancès et al., 1997; Mennucci et al., 1997; Mennucci and Tomasi, 1997; Cossi et al., 1998), including solvation energies for neutral molecules as well as ions. The PCM models have been successfully applied to predict fractionation factors for other isotope systems (e.g., B (Liu and Tossell, 2005) and Fe (Anbar et al., 2005)). We used the IEF-PCM model primarily in this study, but also compared the results with those using the Onsager model and SCI-PCM model and investigated the effects of hydration and H-bond formation by adding 4 innersphere water molecules in conjunction with the Onsager model. We have used the IEF-PCM model as the default solvent model. The key parameter that determines solvent effects is the dielectric constant, ε. Table 4 compares the RPFRs (34β) calculated using the IEF-PCM models for seven sulfur-bearing compounds in the 0–650 °C range. The calculations were

2− Fig. 3. Ratios of fractionation factors between various sulfur-bearing species (HS−(aq)–H2S(g), H2S(aq)–H2S(g), SO2(g)–H2S(g), SO2− 4 (aq)–H2S(g), and SO4 (aq)–H2S(aq)) as a function of temperature (0 to 650 °C) in multiple sulfur isotope systems: (a) (33α − 1) / (34α − 1), (b) ln33α / ln34α, (c) (36α − 1) / (34α − 1), and (d) ln36α / ln34α.

366

T. Otake et al. / Chemical Geology 249 (2008) 357–376

Fig. 4. Large changes in the (33α − 1) / (34α − 1) ratios occurring near the crossover temperatures in some sulfur-bearing species couples: (a) S2–H2S, (b)S8–H2S, (c) CS2–H2S, and (d) SO3–SO2− 4 (aq).

performed for two different conditions where: (a) the dielectric constants were those of water along the liquid–vapor boundary at temperatures below the critical temperature (T = 374 °C), and those at P = 0.5 kbar at higher temperatures (Tanger and Helgeson, 1988); and (b) the dielectric constant of solvent was kept constant throughout the temperature range of 0– 650 °C at 78.39, which is the ε value for water at 25 °C and the default value on the Gaussian 03 package. The calculated RPFRs values were found to be almost identical for each sulfur species under both conditions (a) and (b) at temperatures less than 300 °C; they differ slightly at higher temperatures, but the differences are less than 0.5‰ (Table 4). Fig. 5 shows the calculated fractionation factors (34α) between the seven gaseous species (i.e., ε = 1) and their solvated species (i.e., the IEF-PCM modeled species, ε N 1) at 25 °C, in which the values of ε were varied from 0 to 80. It shows that the magnitudes of fractionation factors do not change considerably in the range of ε = 20 to 80. Sulfate ion gives the largest gassolvent fractionation, ~ − 4.0‰, and it remains constant at ε N 20. However, the fractionation factors rapidly decrease with decreasing ε for ε b 20. This suggests that the calculations of equilibrium fractionation factors involving aqueous species need to incorporate the actual ε values when ε b 20, such as when organic solvents are involved or T N 300 °C for aqueous solutions. However, because the fractionation factors between

gaseous and the IEF-PCM modeled species at T = 300 °C are less than 1.5‰ in the 1000(α − 1) value and become even smaller at higher temperatures (Table 4), we may conclude that, for aqueous solutions, the fractionation factors can be calculated using a constant ε value of 78.39. Fractionations occurring in natural organic solvents (e.g., petroleum), however, involve lower dielectric constants and calculations need to utilize the temperature variation of ε. Note that the data in Fig. 5 and Table 4 enable a calculation of fractionations in organic solvents from the value of the dielectric constant of the solvent at the particular temperature. The Onsager and SCI-PCM model calculations were performed for the comparison with the IEF-PCM model (Table 5). The Onsager model results in the smallest solvent effect on the RPFRs, compared with other SCRF methods. The biggest effect was observed when the IEF-PCM model was applied. Differences in bond lengths between the IEF-PCM and SCI-PCM modeled species were trivial; however, large differences in the vibrational frequencies in some sulfur-bearing compounds (e.g., HS−) were observed. Comparing values of the fractionation factor, the Onsager and SCI-PCM modeled species are less fractionated with respect to the gaseous species (αaq-gas), than IEF-PCM modeled species for all sulfur-bearing compounds and all temperature ranges. Therefore, our IEF-PCM model calculation in the previous sections likely shows the

T. Otake et al. / Chemical Geology 249 (2008) 357–376

maximum solvent effect among the SCRF methods (e.g., ~ 4‰ for SO42− at 25 °C, Fig. 5). While all the SCRF methods take into account long-range electrostatic interactions between solute and solvent, we also need to consider the effect of the short-range interaction (i.e., hydrogen bonding) on the molecular geometry, vibrational frequencies, and RPFRs. Several recent studies (Oi and Yanase, 2001; Liu and Tossell, 2005) have examined the effect of hydration on boron isotope fractionations. While the studies have added from 1 to 34 water molecules, the results have been similar: hydration lowers the RPFR values but the effects are relatively small. We tried to combine the SCRF method and “water-droplet” method by adding 4 water molecules into the Onsager model to include both long-range and short-range water–solute interactions. The results show that some compounds, such as H2S or CS2, were asymmetrically coordinated by 4 water molecules. In particular, SO3 was hydrolyzed and deprotonated, forming HSO4− due to the instability of SO3 as an aqueous species. As a result, RPFRs of SO3 in the model significantly change close to the values for SO42− (Table 5(c)). HS− also shows some increase (~ 0.5‰) in RPFRs, which could be a reason for the inconsistency in the 34 aHS−ðaqÞ H2 S value between our calculation and experimental data. Although this increase is still too small to match the experimental data, the hydration effect may be magnified by adding more water molecules in the system. No large shift in RPFRs was observed in other sulfur-bearing compounds regardless of considerable changes in the vibrational modes and their frequencies. Table 4 The reduced partition function ratios (34β) of sulfur-bearing compounds calculated using the IEF-PCM model under two different assumptions for dielectric constant (ε) of water: (a) ε changes with temperature, and (b) ε is constant at 78.39 (a) T (°C)

0 50 110 170 230 300 400 500 600

ε

1000 (34β − 1) H2 S

HS−

S2

CS2

SO2

SO3

SO2− 4

87.92 69.89 52.96 40.02 29.84 19.85 12.14 3.76 2.26

12.37 10.03 8.06 6.63 5.56 4.61 3.63 2.94 2.42

5.00 4.08 3.30 2.73 2.30 1.92 1.51 1.23 1.01

14.80 11.07 8.18 6.26 4.95 3.86 2.84 2.17 1.71

19.24 14.50 10.78 8.31 6.60 5.20 3.85 3.02 2.43

51.49 40.02 30.64 24.16 19.50 15.57 11.72 9.26 7.49

83.66 64.63 49.26 38.74 31.22 24.89 18.70 14.72 11.87

91.33 69.56 52.23 40.57 32.36 25.57 19.03 14.92 11.97

ε

1000 (34β − 1)

(b) T (°C)

0 50 110 170 230 300 400 500 600

78.39 78.39 78.39 78.39 78.39 78.39 78.39 78.39 78.39

H2 S

HS−

S2

CS2

SO2

SO3

SO2− 4

12.37 10.03 8.05 6.62 5.55 4.60 3.61 2.90 2.38

5.00 4.08 3.30 2.73 2.29 1.91 1.50 1.21 1.00

14.78 11.07 8.18 6.26 4.94 3.86 2.84 2.17 1.71

19.25 14.49 10.77 8.30 6.58 5.17 3.82 2.93 2.32

51.50 40.01 30.62 24.13 19.46 15.50 11.62 9.01 7.18

83.67 64.62 49.23 38.69 31.15 24.79 18.57 14.39 11.46

91.35 69.55 52.19 40.50 32.27 25.44 18.86 14.51 11.49

367

Fig. 5. Fractionation factor (34α) between gaseous species and their solvated species for various sulfur-bearing compounds (S2, HS−, H2S, CS2, SO2, SO3, and SO2− 4 ) at 25 °C. Solvated species are calculated using the IEF-PCM model as a function of dielectric constant (ε) of solvent.

To analyze the solvent effect on the mass-dependent relationships, (33β − 1) / (34β − 1) values for each species are calculated (Table 5(d)). This shows that these solvent models have little effect on the values for most species; however, we observe significant change in the value for H2S when the Onsager + 4 water molecule model was applied (0.502 compared with 0.514 for other models), indicating the possibility that an aqueous phase reaction may vary the mass-dependent relationships. Overall, the solvent effects on fractionation factors, while not very strong, are still significant. For species, such as the sulfate ion, the gas-solvent fractionation can reach values in excess of 4‰. The largest part of the solvent effect occurs as the dielectric constant increases from 1 (gas phase) to values of 20. Once ε exceeds 20, the fractionation factor varies little. Thus, two solvents with dielectric constants both in excess of 20 would have nearly identical fractionation factors at the same temperature for sulfur molecules. Alternatively, for any one solvent, at temperatures below which ε is larger than 20, would not exhibit any additional solvent effects. For example, water at temperatures less than 300 °C would have the same solvent effect on the fractionation factor as water at 25 °C. 5. Evaluation of anharmonic effects on isotope fractionation factors and the mass-dependent relationships Almost all theoretical calculations for isotope fractionation by first-principles methods have been derived from harmonic vibrational frequencies that are multiplied by an appropriate scaling factor to correct for anharmonic effects. Anharmonic effects arise because the potential energy describing variations from equilibrium is not a simple quadratic function of the coordinate changes. For a diatomic molecule with one vibrational mode, the second-order correction to the vibrational energy is given by: "

#   1 1 2 En ¼ hc G0 þ n þ xe  n þ x e xe 2 2

ð16Þ

368

T. Otake et al. / Chemical Geology 249 (2008) 357–376

Table 5 (a) Bond lengths, (b) vibrational frequencies (32S compounds) without scaling, (c) reduced partition function ratios ((1000(34β − 1)) at 25 °C, and (d) (33β − 1) / (34β − 1) ratios for various sulfur-bearing compounds (H2S, HS−, S2, CS2, SO2, SO3, and SO2− 4 ) calculated using the Onsager, Onsager + 4 water molecules, IEF-PCM, and SCIPCM models (a)

Gas Onsager IEF-PCM SCI-PCM

HS−(dSH)

H2S(dSH)

S2(dSS)

CS2(dSC)

SO2(dSO)

SO3(dSO)

SO2− 4 (dSO)

1.342 1.341 1.340 1.337

1.326 1.326 1.330 1.325

1.878 1.878 1.877 1.877

1.545 1.545 1.545 1.545

1.414 1.415 1.419 1.418

1.405 1.405 1.405 1.404

1.487 1.487 1.479 1.480

(b) Frequencies (cm− 1) H2S Gas Onsager IEFPCM SCIPCM

HS−

S2

CS2

SO2

SO2− 4

SO3

1368 1366 1338

2920 2923 2860

2932 2935 2872

2755 2764 2726

813 813 813

445 445 444

725 725 722

1583 1518 1426

592 588 578

1359 1355 1345

1569 1543 1501

562 553 540

583 577 565

1212 1212 1212

1554 1541 1519

469 469 472

668 664 651

1025 1025 1064

1214 1196 1214

1413

2929

2946

2321

815

429

723

1564

602

1360

1543

500

582

1206

1547

457

653

1050

1213

(c)

Gas Onsager Onsager + 4H2O IEF-PCM SCI-PCM

H2S

HS−

S2

CS2

SO2

SO3

SO2− 4

11.27 11.27 11.65 11.10 11.41

4.44 4.46 6.08 4.50 3.60

12.75 12.75 12.78 12.73 12.80

17.52 17.19 17.82 16.62 17.28

47.79 46.76 49.18 45.22 47.40

77.19 75.56 82.76 73.22 74.91

83.91 81.26 81.01 79.36 81.16

H2S

HS−

S2

CS2

SO2

SO3

SO2− 4

(d)

Gas Onsager Onsager + 4H2O IEF-PCM SCI-PCM

0.5144 0.5144 0.5017 0.5144 0.5143

0.5152 0.5152 0.5150 0.5152 0.5154

0.5138 0.5138 0.5140 0.5138 0.5137

where G0 is a constant (cm− 1) and the factor ωexe is the anharmonic correction to the frequency (also in cm− 1). For a polyatomic molecule with f degree of freedom, Eq. (16) generalizes to:  1 En ¼ hcG0 þ hcxi ni þ 2 i   f X 1 1 þ hcxij ni þ nj þ 2 2 jzi f X

ð17Þ

where ωi (i = 1, 2, … f) is the harmonic frequency (xij (in cm− 1) are the coupling factors), and for each quantum energy level, En, n comprises the individual integer quantum numbers, n1, n2, …, nf for each normal mode of vibration. The last term of Eq. (17) arises from the anharmonic coupling among the various normal modes. Note the change in sign between Eqs. (16) and (17) and the omission of ωi in the xij terms in the usual notation (Wolfsberg, 1969; Richet et al., 1977).

0.5131 0.5132 0.5131 0.5132 0.5131

0.5092 0.5094 0.5091 0.5096 0.5093

0.5056 0.5058 0.5050 0.5061 0.5058

0.5048 0.5051 0.5052 0.5054 0.5051

At the simplest level, the significance of anharmonic effects on isotopic fractionation factors can be studied by analyzing the frequency scaling constants obtained within an anharmonic context. A very important contribution to the vibrational partition function comes from the zero-point energy (ZPE) of vibration, E0. E0 for a diatomic is obtained from Eq. (16) by:

1 1 E0 ¼ hcxe  hcxe xe 2 4

ð18Þ

Similarly, the ZPE from Eq. (17) for a polyatomic molecule is given by:

1 X 1 X E0 ¼ hc xi þ hc xij 2 4 i jzi f

f

ð19Þ

T. Otake et al. / Chemical Geology 249 (2008) 357–376

369

Once the ab initio frequencies are scaled, they are usually used in a harmonic calculation of partition functions and fractionation factors. The problem is that, if anharmonic effects are important, the ab initio ZPE being calculated from the scaled frequencies in the harmonic approximation would be (using Eqs. (20) and (21)): 1 1 0 For diatomic Eab ¼ hcxfund ¼ hcxe  hcxe xe 2 2

ð22Þ

1 X fund 0 Eab ¼ hc xi 2 i X 1 X 1 X ¼ hc xi þ hc xii þ hc xij 2 2 i i jNi

For polyatomic Fig. 6. Potential curve of HS− molecule as a function of derivation of the bond length from the optimum value, computed using the B3LYP/6-31G(d) method for anharmonic correction. 1 Hartree = 4.3597482 × 10− 18 J.

The usual scaling of ab initio frequencies strives to obtain the fundamental frequencies (ωfund), which from Eqs. (16) and (17), is given by: For diatomic xfund ¼

E ð 1Þ  E ð 0Þ ¼ xe  2xe xe hc

¼ For polyatomic xfund i 1X xij 2 jpi f

¼ xi þ 2xii þ

ð20Þ

E ð ni ¼ 1Þ  E ð 0Þ hc ð21Þ

Note that xij = xji and the change in summation index in Eq. (21).

ð23Þ

Note the change in summation. Eqs. (22) and (23) underestimate the real ZPE as given P by Eqs. (18) and (19) by  34 hcxe xe and 34 xii hc þ 14 hc jNi xij for diatomic and polyatomic molecules, respectively, as pointed out by Scott and Radom (1996). If one wants to optimize the calculated ZPE as opposed to the fundamental frequencies, then the scaling factors need to be slightly modified. Scott and Radom (1996) used a set of 39 small molecules to recalculate the frequency scaling factors by having the ZPE calculated from the scaled ωifund reproduce the actual ZPE including anharmonic factors. The results are that the scaling factor changes from 0.8929 to 0.9135 for HF/6-31G(d), from 0.9085 to 09248 for HF/6-311G(d,p), and from 0.9613 to 0.9806 for B3LYP/6-31G(d). The data in the Supplementary material (Fig. S1) compare the RPFRs and fractionation factors using these two different (frequency-

Table 6 Results of anharmonic calculations of harmonic frequency (ωe), anharmonic correction (ωexe), fundamental frequency (ωfund), and Go factor for: (a) HS− and (b) H2S, and (c) coupling factors (xij) for H2S (a) HS− H32S− H33S− H34S− H36S−

ωe(cm− 1)

ωexe(cm− 1)

ωfund(a) (cm− 1)

G0(cm− 1)

2576.596 2575.403 2574.282 2572.220

84.306 84.228 84.155 84.020

2407.984 2406.947 2405.973 2404.180

−6.868 −6.861 −6.855 −6.845

(b) H2S ωe(cm− 1) H32 2 S H33 2 S H34 2 S H36 2 S

Δωanh(cm− 1)

2692.464 2691.332 2690.269 2688.316

1249.948 1249.349 1248.786 1247.746

2712.858 2711.543 2710.309 2708.038

149.562 149.429 149.305 149.077

ωfund (cm− 1) 24.549 24.531 24.515 24.485

153.486 153.336 153.196 152.939

2542.902 2541.903 2540.964 2539.239

G0(cm− 1) 1225.399 1224.818 1224.271 1223.261

2559.372 2558.207 2557.113 2555.099

−5.470 −5.468 −5.467 −5.464

(c) xij factors (cm− 1) for H2S (ω1 = 2692, ω2 = 1249, ω3 = 2712) H32 2 S H34 2 S (b) H34 2 S

x11

x12

x13

x22

x23

x33

− 35.467 − 35.403 − 35.409

− 15.513 − 15.492 − 15.486

− 141.744 − 141.504 − 141.495

− 3.439 − 3.433 − 3.433

− 19.828 − 19.808 − 19.791

− 36.350 − 36.27 − 36.282

The calculations were carried out using B3LYP/6-31G(d) method. (a) ωfund = ωe − 2ωexe. (b) Estimated using Darling-Dennsion Approximation.

370

T. Otake et al. / Chemical Geology 249 (2008) 357–376

corrected and ZPE-corrected) scaling factors. In general, the RPFRs using scaling factors for ZPE correction (Scott and Radom, 1996) increase by 0.2–0.6‰ in the 1000(34β − 1) value except for very high values (N 20‰) where the increase could be greater than 1‰. Similar to the β factors, the fractionation factors can vary by several ‰ when fractionation factors are greater than 20‰. In the case of SO42−(aq)–H2S, the ZPEcorrected fractionation factor is in slightly better agreement with the experiment than the frequency-corrected fractionation factor for the HF/6-31G(d) (0.5‰ vs − 0.7‰). The fact that the two different scaling factors yield slightly different RPFRs suggests the significance of anharmonicities in the calculation of RPFRs and fractionation factors. Ultimately, one needs to carry out the ab initio calculations to retrieve the actual G0, ωexe, and xij factors in Eqs. (16) and (17). Scaling can handle most cases, but unusual anharmonic effects will only be understood with the full analyses. Note that, in general, anharmonic effects tend to yield average bond lengths at 0 K, which are slightly higher than (0.01 Å or so) the harmonic values. Interestingly, these anharmonic corrections would make the HF/6-31G(d) agreement better yet than B3LYP methods based on the earlier section regarding the geometry of our sulfur molecules. Vibrational analysis now requires extraction of the G0, ωexe, and xij factors (see Appendix B). Table 7 34 − Reduced partition function ratios
( β) for (a) HS , (b) H2S, (c) fractionation − 34 − S and HS a Þ , and (d) fractionation factor ratios factors between H HS −H S 2 2
33
aHS− −H2 S  1Þ= 34 aHS− −H2 S  1Þ calculated using the B3LYP/6-31G(d) method (a) 1000(34β − 1) for HS− T (°C)

0

50

110

150

190

β β 34 scale β

5.21 4.94 4.97

4.26 4.04 4.06

3.45 3.26 3.28

3.04 2.87 2.89

2.70 2.54 2.56

T (°C)

0

50

110

150

190

β β 34 scale β 34 SR β 34 DD β

12.96 12.37 12.35 12.65 12.35

10.54 10.04 10.03 10.28 10.02

8.48 8.06 8.05 8.27 8.05

7.44 7.06 7.06 7.25 7.05

6.59 6.25 6.24 6.41 6.23

T (°C)

0

50

110

150

190

harmonic anharmonic scaled

− 7.65 − 7.33 − 7.29

− 6.21 − 5.94 − 5.90

− 4.98 − 4.76 − 4.73

−4.37 −4.16 −4.14

− 3.86 − 3.68 − 3.65

34 h

34 anh

(b) 1000(34β − 1) for H2S 34 h

34 anh

In this section, we show results of the higher-level calculation, instead of HF/6-31G(d), because without scaling, the higher-level method performs much better. Results of calculation with other methods are shown in the Supplementary material (Tables S8, S9 and S10). Fig. 6 gives the potential energy of the HS− molecule, V, as a function of the deviation, Δr, of the bond length from the optimum value evaluated using B3LYP/6-31G(d). These V versus Δr data can be used to extract both G0 and ωexe factors using the equations in the Appendix B. While the Gaussian 03 program carries out anharmonic calculations similar to these (Barone, 2004, 2005), it does not do so for linear molecules, including diatomic molecules (or for spherical top molecules). Hence we have carried out the anharmonic calculations for HS− using the ab initio energies directly (see Appendix B). Results of the anharmonic corrections are shown in Table 6 for the B3LYP/6-31G(d) method and the Supplementary material for other methods (Table S8). The ωexe anharmonic correction ranges from ~ 51 to ~ 84 cm− 1. The comparison of the fundamental frequency, ω e − 2ω e x e (Eq. (20)), to the experimental value for HS − and H2S demonstrates that now the higher-level methods without scaling are much closer than the lower-level methods (Table S8). Thus, if scaling is not used, accuracy increases as both higher basis sets and electron correlation as well as anharmonicites are included. Of interest, nonetheless, is that the higher-level methods tend to underestimate the fundamental frequency. For polyatomic molecules, the variation in energy as a function of all f normal modes must be fit numerically to obtain coefficients such as in Eq. (A1) in Appendix B. In this case, we also retrieve coupling coefficients between two different normal modes, which will yield the xij parameters in Eq. (17). The xij coefficients for H2S are given in Table 6 for the 3 normal modes of vibration calculated using the B3LYP/6-31G(d) method. Note that the xij are usually negative in agreement with the sign in Eq. (16) for a diatomic. The biggest coupling factor for H2S is x13, which shows anharmonic coupling between the two H–S bond stretching vibrational modes (the symmetric and antisymmetric stretches). Note also that the Darling-Dennison Approximation (Richet et al. (1977)), which is used to estimate xij for k S-substituted molecule (kxij) from xij for 32S-bearing molecule: k


(c)1000 34 aHS H2 S  1Þ

(d) aHS− −H2 S ratio T (°C)

0

50

110

150

190

α − 1/34αh − 1 33 anh α − 1/34αanh − 1 33 scale α − 1/34αscale − 1

0.5168 0.5167 0.5167

0.5167 0.5165 0.5165

0.5165 0.5163 0.5164

0.5164 0.5163 0.5163

0.5164 0.5162 0.5163

33 h

β : unscaled harmonic RPFRs; 34βanh: full anharmonic RPFRs (no scaling); β : frequency-scaled RPFRs (our default method); 34βSR: ZPE-scaled RPFRs; 34βDD: anharmonic RPFRs calculated using the Dennison-Darling rule to estimate xij.

34 h

xij ¼ xij

k

xi k xj xi xj

ð24Þ

is only roughly followed. For 34xij, the error in the shift 34xij−32xij from the Darling-Dennison Approximation is 9.4, 28.6, 3.8, 0, 8.5, and 15% for x11, x12, x13, x22, x23, and x33, respectively. The effect of anharmonicities on the RPFRs, fractionation factors, and their ratios is evaluated using the B3LYP/6-31G (d) method in Table 7. Here, RPFRs were calculated as follows: i

 i ! X  ! X jj i Uj  En E n b¼ Uj exp exp 32 = ð25Þ kT kT jj 32 n n

34 scale

where iUj = hνj / kT for the iS-substituted molecule as defined earlier, and iEn is the nth energy level as given by Eq. (17) for

T. Otake et al. / Chemical Geology 249 (2008) 357–376

371

unscaled (harmonic) RPFRs, 34βh, are 0.2–0.3‰ bigger in the 1000(34β − 1) value than RPFRs scaled with the ZPE-corrected scaling factor (0.9806) of Scott and Radom (1996), 34βSR. However, RPFRs using the full anharmonic results (i.e., Eq. (25) with Eq. (17)) and no scaling, 34βanh, are still 0.2–0.3‰ lower than the 34βSR values. In fact, the 34βanh values are very close to the 34β values using the standard scaling factor (0.9613), 34βscale, as employed in our earlier section. On the other hand, the difference between the 34βanh values, which is using the actual xij, and the 34β values employing the DennisonDarling rule, 34βDD, is negligible. A universal result of anharmonicity is to decrease the RPFRs (Table 7). The RPFRs for H2S, however, are decreased (~ 0.5‰ for the average in the 0–190 °C range) much more than those for HS− (~ 0.1‰ for the same value). Therefore, the 34 aH2 SHS values are decreased overall, improving agreement with the experimental data, especially for the higher-level calculations. The G0 values affect the 34 aH2 SHS values only slightly by up to 0.05‰ (Table S9). At this point, the remaining persistent discrepancy for the H2S–HS − fractionation between the quantum results and experiments may suggest that the experimental values could be too low. In fact, another experimental study (Tudge and Thode, 1950) indicates that the fractionation between H2S and HS− occurs ~ 6‰ at 25 °C, which is more consistent with our results. Table 7(d) gives the anharmonic effects on the (33α − 1) / (34α − 1) ratios, showing there is almost no difference in the ratios between harmonic and anharmonic calculations. Therefore, anharmonicity for normal molecules can significantly affect RPFRs and fractionation factors, but has little effect on the mass-dependent relationships, at least for the H2S–HS− system. 6. Synthesis: effects of temperature, frequency, and mass of atoms bonded to sulfur on isotopic fractionation

Fig. 7. The harmonic reduced partition function ratios (1000(34βh − 1)) as a function of 32U, calculated using Eq. (26) (solid line), Eq. (31) (dot-dashed line), and Eq. (32) (dashed line) for different masses of an atom bounded to sulfur (a, M = 32; b, M = 16; c, M = 7).

the i S-substituted molecule. Note that the implicit summation in Eq. (25) is really a multiple sum over each of the normal modes of vibration. In the case of H2S with 3 vibrational modes, the sum over n implies a sum over n1, n2, n3 each individually going from 0, 1, 2, … infinity. For each value of n1, n2, n3, the value of En is as given by Eq. (17). For all our molecules, the summation would converge accurately if the ni just ranges from 0 to 6 (i.e., 73 or 343 terms). Eq. (25) still assumes the RedlichTeller rule (Richet et al., 1977) but now uses the full Eq. (17) for En to obtain the quantum vibrational partition function. (The Redlich-Teller rule allows us to eliminate the ratio of translational and rotational partition functions (which are classical) in terms of the harmonic vibrational terms Uj in Eq. (25) as also employed by Bigeleisen and Mayer (1947)). The RPFRs, 34β, for H2S including anharmonic effects and using the B3LYP/6-31G(d) method are given in Table 7. The

At this stage, it is useful to summarize the overall effect of key variables (temperature, vibrational frequency, and mass of atoms bonded to sulfur atom) on the RPFRs and massdependent relationships. A convenient way to display the effects is to use the diatomic case or equivalently focus on the bond between sulfur and another chemical entity. At the outset, to analyze the effect on either 34β or (33β − 1)/(34β − 1) values, it is easiest to use the harmonic approximation: b ¼

i h

U expði U =2Þð1  expð32 U ÞÞ 32 U expð32 U =2Þð1  expði U ÞÞ

i

ð26Þ

Eq. (26) is exact within the harmonic approximation, where pffiffiffiffiffiffiffiffiffi 1 k=i A . k is the mass-independent U = hciωe / kT and i xe ¼ 2pc bond force constant and so i

i

U ¼

qffiffiffiffiffiffiffi 32 A =i A 32 U

ð27Þ

The most accurate approximation (Bigeleisen and Mayer, 1947) sets 32U = iU + Δi U and correctly assumes that in all cases

372 T. Otake et al. / Chemical Geology 249 (2008) 357–376 Fig. 8. The dependence of RPFRs (34β values, a–c) and (33β − 1) / (34β − 1) values (d–f) on the three principal parameters: temperature (T), harmonic frequency (ωe), and mass of atoms bounded to sulfur (M). (a), (d): as a function of M for ωe = 500, 1500, and 2500 cm− 1 at T = 50 °C. (b), (e): as a function of ωe for M = 1, 7, 12., 18, 32, and 58 at T = 50 °C. (c), (f): as a function of ωe for M = 16 and T = 0, 50, 100, 150, and 200 °C.

T. Otake et al. / Chemical Geology 249 (2008) 357–376

ΔiU « 1. To a fairly high degree of accuracy, therefore, Eq. (26) is simplified to ! 1 1 1 i h  U þ i  1 Di U b ¼1þ ð28Þ 2 i eU Using Eq. (27) and defining qffiffiffiffiffiffiffi i g ¼ 32 A =i A

ð29Þ

then
 32 1  i g 32 1
i b ¼2 i gþ 1 g U þ U  132 U i 2 eg

i h

1

ð30Þ

Note that iγ ≈ 1 for all cases. The literature uses a further approximation to Eq. (30) to derive the mass-dependent approximations (i.e., Eq. (7)), namely, that the dimensionless number 32U itself is a small number (Bigeleisen and Mayer, 1947; Schauble, 2004). Because RT = 208.5 cm− 1 at 300 K, 32 U b 1 assumes ωe b~ 208.5 cm− 1. This approximation is not on the same level as the first approximation (ΔU « 1). In fact, there are two different limits to Eq. (30), the high 32U and low 32 U limit. It can be readily shown that:   1 1
b Y 2  i g þ 1  i g 32 U 2

i h

b Y 1þ

i h

  1
1  i g i g 32 U Þ2 12

U r1

U T1

ð31Þ

ð32Þ

Because 32U = hc32ωe / kT, Eqs. (31) and (32) predict that iβh will vary as 1 / T for low T and high ωe (high U), and as 1 / T2 for high T and low ωe (low U). Fig. 7 shows the variation of harmonic RPFRs with 32U for different masses, M, bonded to the sulfur atom: for example, M = 1 amu for H, 7 amu for Li, 12 amu for C, 16 amu for O, 19 amu for F, 32 amu for S, and 56 amu for Fe. The mass of atoms bonded to the sulfur atom, M, will turn out to be a rather important quantity as we shall see below. M affects the reduced mass via:
 i A ¼ i mM= i m þ M ð33Þ The lines and curves in Fig. 7 are given by Eqs. (31) and (32). It is clear that Eq. (31) predicts the behavior fairly well if 32 U ≥ ~ 3, while Eq. (32) is fairly accurate if 32 U ≤ ~ 3, regardless of mass, M. The fundamental frequencies of all except the heaviest molecules are in the 400–4000 cm− 1 range. The boundary between the 1 / T and 1 / T2 regions for iβh implies that (using U = hcω / kT = 3):  1
Tboundary ð K Þi xe cm1 2
 1=T region at T bTboundary ; 1=T 2 region at T NTboundary ð34Þ Thus, for 32ωe N 1000 cm− 1, the linear portion dominates the isotope fractionation at temperatures below 200 °C. However,

373

for 32ωe b 600 cm− 1, the quadratic portion dominates in the same temperature range (0–200 °C). For molecules with ωe both above and below 1000 cm− 1, then the iβ will behave as a + b / T + c / T2 (e.g., see Table 3). As discussed above, the magnitude of harmonic RPFRs (iβ) is a strong function of temperature (T), the mass (M), and the harmonic frequency (ωe). Fig. 8a and b show the dependency of 1000(34β − 1) values as a function of M and ωe at T = 50 °C, and Fig. 8c shows the 1000(34β − 1) values as a function of ωe and T for a fixed M = 16 amu. In general, they increase with increasing M and ωe and decreasing T. For example, a decrease in 1000 (34β − 1) value from 40 to 20 may be caused for molecules with M = 16 amu by a decrease in ωe from ~ 2200 to ~ 1400 cm− 1 at T = 50 °C (Fig. 8b), or by an increase in T from 0° to 150 °C at ωe = ~ 1800 cm− 1 (Fig. 8c); the same change in 34β value may be caused by a decrease in M from 16 to 7 amu at T = 50 °C and ωe =~ 2200 cm− 1 (Fig. 8a). Another application of Eqs. (26), (27), and Fig. 8 is to estimate the error in 34β value from errors in ωe and, therefore, errors in 32U and 34U. For a 10 cm− 1 error in ωe, the error in 34β (as 1000Δ34β in ‰) would be only 0.006–0.017‰ when M = 1 amu and ωe = 500 cm− 1 (for T = 0 to 200 °C); the error in 34 β increases to 0.013–0.024‰ for ωe = 2000 cm− 1. However, as the mass (M) increases, the error in 34β also increases. A 10 cm− 1 error for M = 16 amu would lead to a 0.07–0.186‰ error in 34β (for T = 0 to 200 °C) if ωe = 500 cm− 1 and to a 0.149–0.270‰ error in 34β if ωe = 2000 cm− 1. Of further interest is the role of these three variables (T, M, and ωe) on the mass-dependent relationship, (33β − 1) / (34β − 1) value (Fig. 8d–f). Fig. 8d–f show the (33β − 1) / (34β − 1) values under same conditions as Fig. 8a–c, respectively. The (33β − 1) / (34β − 1) values converge to 0.516 at higher T (N 600 °C), lower ωe (b500 cm− 1), and smaller M (b 7 amu); however, they decrease at lower T, higher ωe, and larger M. For example, the (33β − 1) / (34β − 1) value becomes as low as 0.503 at T = 50 °C, ωe = 2500 cm− 1, and M = 56 amu (Fig. 8d and e). Therefore, these three variables significantly affect both 34β values and (33β − 1) / (34β − 1) value. By decreasing T, increasing ωe and M, the 1000(34β − 1) values increase and the (33β − 1) / (34β − 1) values decrease (i.e., deviate from “expected” 0.515). We have also analyzed anharmonic effects on both 34β and 33 ( β − 1) / (34β − 1) values; the results are presented in Table S11 in the Supplementary material. As discussed earlier, the anharmonicity uniformly decreases the RPFRs. This decrease is larger for molecules with higher masses. The effect of ωexe on the RPFRs is very close to linear (Table S11): 34

b ¼ 34 b h  a  xe xe

ð35Þ

The value of the slope, a, varies with temperature and mass but is little affected by frequency. The decrease in 34β value for ωexe = 100 cm− 1 in the temperature range of 0 to 200 °C is 0.1 to 0.2‰ for M = 1 amu, 1.5 to 2.7‰ for M = 16 amu, and 2.2 to 4.2‰ for M = 32 amu. In contrast, the effect of anharmonicities is rather small for the (33β − 1) / (34β − 1) value. Anharmonicities tend to increase the (33β − 1) / (34β − 1) value, with the increase being more significant for the cases of large mass (Table S11).

374

T. Otake et al. / Chemical Geology 249 (2008) 357–376

For example, for M = 32 amu, ωe = 2500 cm− 1, and at T = 50, the (33β − 1) / (34β − 1) value is 0.5060 under harmonic conditions and increase only to 0.5067 at ωexe = 200 cm− 1. 7. Summary Important findings from this study are summarized below: 1. Among the various ab initio methods, both HF and hybrid HF-DFT methods yield harmonic vibrational frequencies of simple sulfur-bearing compounds, RPFRs and fractionation factors (34α) that are reasonably consistent with experimentally determined data (b 50 cm− 1 for vibrational frequencies and b 5‰ for fractionation factors). In particular, the HF method with a 6-31G(d) basis set and scaling is preferable for the calculations in sulfur isotope system. 2. Multiple sulfur isotope fractionation factors (33α, 34α, and 36 α) calculated by the first-principle methods were examined for verification of the mass-dependent approximations: ( 33 α − 1) / ( 34 α − 1) =0.515 and ( 36 α − 1) / ( 34 α − 1) = 1.89. Our results show that these approximations are accurate at temperatures higher than 500 °C, but become less accurate with decreasing temperatures. The (33α − 1) / (34α − 1) and (36α − 1) / (34α − 1) values range from 0.505 to 0.517 and from 1.88 to 1.96, respectively, at 0 °C. The larger variations occur at temperature below ~ 250 °C, especially when heavier species (e.g., SO42−(aq)) are involved. 3. Crossover of sulfur isotope fractionation occurs in several species couples (e.g., H2S–S2, H2S–S8, H2S–CS2, SO3–SO42 − (aq)), and possibly in many other systems where the fractionation factors between the pairs are small and one of the pairs contains hydrogen atoms. Differences in the crossover temperatures for different isotope molecules result in large deviations from the mass-dependent approximations within a 10 °C range about the crossover temperatures. However, Δ33S (and Δ36S) values expected from the crossover are typically small (b 0.01‰) because fractionation factors themselves are very small (α ≈ 1.0). 4. Fractionation factors between gaseous and their solvated species (e.g., SO2(g)–SO2(aq); H2S(g)–H2S(aq); H2S(g)– H2S(oil)) were estimated as a function of dielectric constant (ε) of solvent using the IEF-PCM model. The calculated fractionation factors (34α) are as large as 4‰ in the case of sulfate ion. The fractionation factors do not change considerably in the range of ε = 20 to 80; however, they rapidly decrease with decreasing ε for ε b 20. The IEF-PCM model, compared to other SCRF methods (e.g., the Onsager and SCI-PCM models), shows the maximum solvent effects on the RPFRs. When 4 water molecules are added to the Onsager model, no significant hydration effects on the fractionation factors and mass-dependent relationships were observed for the other sulfur-bearing compounds, except for H2S and SO3. 5. When anharmonicity is corrected without scaling for HSand H2S, the DFT methods yield frequency values closer to the experimental values for both molecules. Anharmonicity can significantly affect the RPFRs (0.5‰ for the 1000

(34βH2S − 1) value and 0.1‰ for the 1000(34βHS- − 1) value) and fractionation factors; however, it has little effect on the fractionation factor ratios. Interestingly, the scaled 34β values actually agree well with the full anharmonic 34β values. 6. Temperature (T), harmonic frequency (ωe), and mass of atoms bonded to sulfur (M) are the three important parameters affecting the RPFRs, fractionation factors, and mass-dependent relationships. By decreasing T and increasing ωe and M, the 1000( 34 β − 1) values increase and ( 33 β − 1) / ( 34 β − 1) values decrease (i.e., deviate from “expected” 0.515). The (33β − 1) / (34β − 1) values converge to 0.516 only at higher T (N600 °C), lower ωe (b 500 cm− 1), and smaller M (b7) and becomes as low as 0.503 at T = 50 °C, ωe = 2500 cm− 1, and M = 56. Acknowledgments We thank P. Deines, J. Kubicki, T. Fujii, Y. Watanabe, E. Altinok, and D. Bevacqua for the comments and suggestions on the early manuscript. We also thank J. Tossell and E. Schauble who provided constructive reviews of the early manuscript. The ab initio calculations were carried out using computers in HPC (High Performance Computing) group at Penn State. This research was supported by grants from NASA Astrobiology Institute (NCC2-1057l; NNA04CC06A) and NSF (EAR0229556) to HO. Appendix A. Supplementary data Supplementary data associated with this article can be found, in the online version, at doi:10.1016/j.chemgeo.2008.01.020. Appendix B. Anharmonic calculations using ab initio energies for a diatomic molecule For a diatomic molecule, the two new factors, G0, which corrects the ZPE, and ωexe, arise if the potential energy as a function of deviation from equilibrium for the one given normal mode, Δq, is not quadratic but actually follows: V ðrÞ ¼ V0 þ aDq2 þ bDq3 þ cDq4

ðA1Þ

where the deviation, Δq, is the shift in the particular normal coordinate from equilibrium (Wolfsberg, 1969). When the equilibrium geometry is indeed a minimum, there is essentially no Δq term (linear term). If the two atoms are along the z-axis, the normalized normal mode coordinate is given by (only including the z-direction): 1

Q¼

½m2 =ðm1 þ m2 Þ2 1 ½m1 =ðm1 þ m2 Þ2

! ðA2Þ

where m1 and m2 are the two masses of the atoms. Thus, Dq Dz1 ¼ pffiffiffiffiffiffi 1 m1 ½m2 =ðm1 þ m2 Þ2 ;

ðA3Þ

T. Otake et al. / Chemical Geology 249 (2008) 357–376 1 Dq Dz2 ¼ pffiffiffiffiffiffi½m1 =ðm1 þ m2 Þ2 ; m2

ðA4Þ

Dq and Dr ¼ Dz1  Dz2 ¼ pffiffiffi A

ðA5Þ

where the usual reduced mass, μ, is defined as μ = m1m2 / (m1 + m2) and Δq is the shift in Åamu1/2. By fitting a quartic polynomial to the data of energy versus Δq, the coefficients a, b, and c of Eq. (A1) can be determined. For each sulfur mass, the data are turned into a V versus Δq plot using Eq. (A5), which yields directly the a, b, and c coefficients. The coefficients yield the required parameters in Eq. (16) in the text. The a coefficient should yield the usual harmonic frequency. In the case that the energy is calculated in Hartrees, distance is in bohrs, and the mass is in amu, then pffiffiffiffiffi xe ¼ 5140:48707 2a

and

G0

¼

21 1 gþ d 2 2

xe xe ¼ 90g  2d

ðA6Þ

ðA7Þ

ðA8Þ

where g ¼ 3:5029187  1015 b2 =x4e ðcm1 Þ, d ¼ 2:386129575 109 c=x2e ðcm1 Þ (Wolfsberg, 1969). References Anbar, A.D., Jarzecki, A.A., Spiro, T.G., 2005. Theoretical investigation of iron 2+ isotope fractionation between Fe(H2O)3+ 6 and Fe(H2O)6 : implications for iron stable isotope geochemistry. Geochimica et Cosmochimica Acta 69, 825–837. Barone, V., 2004. Vibrational zero-point energies and thermodynamic functions beyond the harmonic approximation. Journal of Chemical Physics 120, 3059–3065. Barone, V., 2005. Anharmonic vibrational properties by a fully automated second-order perturbative approach. Journal of Chemical Physics 122, 014108. Becke, A.D., 1993. Density-functional thermochemistry. III. The role of exact exchange. Journal of Chemical Physics 98, 5648–5652. Bigeleisen, J., Mayer, M.G., 1947. Calculation of equilibrium constants for isotopic exchange reactions. Journal of Chemical Physics 15, 261–267. Cancès, E., Mennucci, B., Tomasi, J., 1997. A new integral equation formalism for the polarizable continuum model: theoretical background and applications to isotropic and anisotropic dielectrics. Journal of Chemical Physics 107, 3032–3041. Canfield, D.E., Thamdrup, B., 1994. The production of 34S-depleted sulfide during bacterial disproportionation of elemental sulfur. Science, 266, 1973–1975. Chacko, T., Cole, D.R., Horita, J., 2001. Equilibrium oxygen, hydrogen and carbon isotope fractionation factors applicable to geologic systems. In: Valley, J.W., Cole, D.R. (Eds.), Stable Isotope Geochemistry. . Reviews in Mineralogy and Geochemistry. The Mineralogical Society of America, Washington, DC, pp. 1–81. Cossi, M., Barone, V., Mennucci, B., Tomasi, J., 1998. Ab initio study of ionic solutions by a polarizable continuum dielectric model. Chemical Physics Letters 286, 253–260. Cotton, F.A., Wilkinson, G., 1972. Advanced Inorganic Chemistry. Interscience, New York.

375

Curtiss, L.A., Raghavachari, K., Redfern, P.C., Rassolov, V., Pople, J.A., 1998. Gaussian-3 (G3) theory for molecules containing first and second-row atoms. Journal of Chemical Physics 109, 7764–7776. Deines, P., 2003. A note on intra-elemental isotope effects and the interpretation of non-mass-dependent isotope variations. Chemical Geology 199, 179–182. Detmers, J., Brüchert, V., Habicht, K.S., Kuever, J., 2001. Diversity of sulfur isotope fractionations by sulfate-reducing prokaryotes. Applied and Environment Microbiology 67, 888–894. Ding, T., et al., 2001. Calibrated sulfur isotope abundance ratios of three IAEA sulfur isotope reference materials and V-CDT with a reassessment of the atomic weight of sulfur. Geochimica et Cosmochimica Acta 65, 2433–2437. Driesner, T., Ha, T.K., Seward, T.M., 2000. Oxygen and hydrogen isotope fractionation by hydration complexes of Li+, Na+, K+, Mg2+, F−, Cl−, and Br−: a theoretical study. Geochimica et Cosmochimica Acta 64, 3007–3033. East, A.L.L., Radom, L., 1997. Ab initio statistical thermodynamical models for the computation of third-law entropies. Journal of Chemical Physics 106, 6655–6674. Edraki, M., Golding, S.D., Baublys, K.A., Lawrence, M.G., 2005. Hydrochemistry, mineralogy and sulfur isotope geochemistry of acid mine drainage at the Mt. Morgan mine environment, Queensland, Australia. Applied Geochemistry 20, 789–805. Farquhar, J., Wing, B.A., 2003. Multiple sulfur isotopes and the evolution of the atmosphere. Earth and Planetary Science Letters 213, 1–13. Farquhar, J., Bao, H., Thiemens, M., 2000. Atmospheric influence of Earth's earliest sulfur cycle. Science, 289, 756–758. Farquhar, J., Savarino, J., Airieau, S., Thiemens, M.H., 2001. Observation of wavelength-sensitive mass-independent sulfur isotope effects during SO2 photolysis: implications for the early atmosphere. Journal of Geophysical Research 106 (E12), 32829–32839. Farquhar, J., et al., 2003. Multiple sulphur isotopic interpretations of biosynthetic pathways: implications for biological signatures in the sulphur isotope record. Geobiology, 1, 27–36. Foresman, J.B., Frisch, A.E., 1996. Exploring Chemistry with Electronic Structure Methods. Gaussian, Inc., Pittsburgh, PA. Frisch, M.J., et al., 2004. Gaussian 03, Revision C.01. Gaussian, Inc., Wallingford, CT. Hehre, W.J., Radom, L., Schleyer, P.V.R., Pople, J.A., 1986. An initio molecular orbital theory. 548 pp. Herzberg, G., 1966. Electronic spectra and electronic structure of polyatomic molecules. Van Nostrand, New York. 745 pp. Hester, R.E., Plane, R.A., 1963. Raman spectra of aqueous solutions of indium sulfate, nitrate, and perchlorate. Journal of Chemical Physics 38, 249–250. Hoering, T.C., 1989. The isotopic composition of bedded barites from the Archean of southern India. Journal of the Geological Society of India, 34, 461–466. Huber, K.P., Herzberg, G., 1979. Molecular spectra and molecular structure IV. Constants for Diatomic Molecules. Van Nostrand Reinhold, New York. 716 pp. Hulston, J.R., Thode, H.G., 1965. Variations in the S33, S34, and S36 contents of metorites and their relation to chemical and nuclear effects. Journal of Geophysical Research, 70, 3475–3484. Jamieson, J.W., Wing, B.A., Hannington, M.D., Farquhar, J., 2006. Evaluating isotopic equilibrium among sulfide mineral pairs in Archean ore deposits: case study from the Kidd Creek VMS deposit, Ontario, Canada. Economic Geology and the Bulletin of the Society of Economic Geologists 101, 1055–1061. Jensen, F., 1999. Introduction to Computational Chemistry. Wiley, New York. 429 pp. Johnston, D.T., et al., 2005. Active microbial sulfur disproportionation in the Mesoproterozoic. Science, 310, 1477–1479. Lasaga, A.C., 1990. Atomic treatment of mineral–water surface reactions. In: HochellaJr. Jr., M.F., White, A.F. (Eds.), Mineral–Water Interface Geochemistry. . Reviews in Mineralogy. The Mineralogical Society of America, Washington, DC, pp. 17–85. Lasaga, A.C., 1998. Kinetic Theory in the Earth Sciences. Princeton University Press, Princeton, NJ. 811 pp.

376

T. Otake et al. / Chemical Geology 249 (2008) 357–376

Lee, C., Yang, W., Parr, R.G., 1988. Development of the Colle-Salvetti correlation-energy formula into a functional of the electron density. Physical Review B, 37, 785–789. Levine, I.N., 2000. Quantum Chemistry. Prentice Hall, NJ. 739 pp. Liu, Y., Tossell, J.A., 2005. Ab initio molecular orbital calculations for boron isotope fractionations on boric acid and borates. Geochimica et Cosmochimica Acta 69, 3995–4006. McGinnety, J.A., 1972. Redetermination of the structure of potassium sulphate and potassium chromate: the effect of electrostatic crystal force upon observed bond lengths. Acta Crystallographica Section B: Structural Science 28, 2845–2852. Mennucci, B., Tomasi, J., 1997. Continuum solvation models: a new approach to the problem of solute's charge distribution and cavity boundaries. Journal of Chemical Physics 106, 5151–5158. Mennucci, B., Cancès, E., Tomasi, J., 1997. Evaluation of solvent effects in isotropic and anisotropic dielectrics and in ionic solutions with a unified integral equation method: theoretical bases, computational implementation, and numerical applications. Journal of Physical Chemistry B 101, 10506–10517. Meyer, V., Sutter, D., Dreizler, H., 1991. The centrifugally induced pure rotational spectrum and the structure of sulfur trioxide. A microwave Fourier transform study of a nonpolar molecule. Zeitschrift fur Naturforschung A: Journal of Physical Sciences 46a, 710–714. Ohmoto, H., Goldhaber, M.B., 1997. Isotopes of Sulfur and Carbon. In: Barnes, H.L. (Ed.), Geochemistry of Hydrothermal Ore Deposits. John Wiley & Sons, New York, pp. 517–611. Ohmoto, H., Lasaga, A.C., 1982. Kinetics of reactions between aqueous sulfates and sulfides in hydrothermal systems. Geochimica et Cosmochimica Acta 46, 1727–1745. Oi, T., 2000. Calculations of reduced partition function ratios of monomeric and dimeric boric acids and borates by the ab initio molecular orbital theory. Journal of Nuclear Science and Technology 37, 166–172. Oi, T., Yanase, S., 2001. Calculations of reduced partition function ratios of hydrated monoborate anion by the ab initio molecular orbital theory. Journal of Nuclear Science and Technology 38, 429–432. Oi, T., Taguma, N., Kakihana, H., 1985. Calculations of thermodynamic sulfur isotope effect. Journal of Nuclear Science and Technology 22, 818–832. Ono, S., Wing, B., Johnston, D., Farquhar, J., Rumble, D., 2006. Massdependent fractionation of quadruple stable sulfur isotope system as a new tracer of sulfur biogeochemical cycles. Geochimica et Cosmochimica Acta 70, 2238–2252. Ono, S., Shanks III, W.C., Rouxel, O.J., Rumble, D., 2007. S-33 constraints on the seawater sulfate contribution in modern seafloor hydrothermal vent sulfides. Geochimica et Cosmochimica Acta 71, 1170–1182. Ottonello, G., Zuccolini, M.V., 2005. Ab-initio structure, energy, and stable Cr isotopes equilibrium fractionation of some geochemically relevant H–O– Cr–Cl complexes. Geochimica et Cosmochimica Acta 69, 851–874. Riccardi, A.L., Arthur, M.A., Kump, L.R., 2006. Sulfur isotopic evidence for chemocline upward excursions during the end-Permian mass extinction. Geochimica et Cosmochimica Acta 70, 5740–5752. Richet, P., Bottinga, Y., Javoy, M., 1977. A review of hydrogen, carbon, nitrogen, oxygen, sulphur, and chlorine stable isotope fractionation among gaseous molecules. Annual Review of Earth and Planetary Sciences, 5, 65–110.

Robinson, R.A., Stokes, R.H., 1959. Electrolyte solutions. Academic Press, New York. 559 pp. Rye, R.O., 2005. A review of the stable-isotope geochemistry of sulfate minerals in selected igneous environments and related hydrothermal systems. Chemical Geology 215, 5–36. Sakai, H., 1968. Isotopic properties of sulfur compounds in hydrothermal processes. Geochemical Journal, 2, 29–49. Schauble, E., 2004. Applying stable isotope fractionation theory to new systems. In: Johnson, C.M., Beard, B.L., Albarède, F. (Eds.), Geochemistry of NonTraditional Stable Isotopes. . Reviews in Mineralogy and Geochemistry. The Mineralogical Society of America, Washington, DC, pp. 65–111. Schauble, E., Rossman, G.R., Taylor Jr., H.P., 2004. Theoretical estimates of equilibrium chromium-isotope fractionations. Chemical Geology 205, 99–114. Scott, A.P., Radom, L., 1996. Harmonic vibrational frequencies: an evaluation of Hartree-Fock, Møller-Plesset, Quadratic Configuration Interaction, Density Functional Theory, and semiemprical scale factors. Journal of Physical Chemistry 100, 16502–16513. Seo, J.-H., Lee, S.-K., Lee, I., 2007. Quantum chemical calculations of equilibrium copper (I) isotope fractionations in ore-forming fluid. Chemical Geology 243, 225–237. Shimanouchi, T., 1972. Tables of Molecular Vibrational Frequencies Consolidated, Vol. I, 1. National Bureau of Standards, Washington, DC. 160 pp. Tanger, J.C., Helgeson, H.C., 1988. Calculation of the thermodynamic and transport properties of aqueous species at high pressures and temperatures: revised equations of state for the standard partial molal properties of ions and electrolytes. American Journal of Science, 288, 19–98. Thode, H.G., Rees, C.E., 1971. Measurement of sulphur concentrations and the isotope ratios 33S/32S, 34S/32S and 36S/32S in Apollo 12 samples. Earth and Planetary Science Letters 12, 434–438. Thode, H.G., Cragg, C.B., Hulston, J.R., Rees, C.E., 1971. Sulphur isotope exchange between sulphur dioxide and hydrogen sulphide. Geochimica et Cosmochimica Acta 35, 35–45. Tomasi, J., Persico, M., 1994. Molecular interactions in solution: an overview of methods based on continuous distributions of the solvent. Chemical Reviews 94, 2027–2094. Tossell, J.A., 2005. Calculating the partitioning of the isotopes of Mo between oxidic and sulfidic species in aqueous solution. Geochimica et Cosmochimica Acta 69, 2981–2993. Tudge, A.P., Thode, H.G., 1950. Thermodynamic properties of isotopic compounds of sulphur. Canadian Journal of Research. Section B 28, 567–578. Urey, H.C., 1947. The thermodynamic properties of isotopic substances. Journal of Chemical Society 562–581. Wolfsberg, M., 1969. Correction to the effect of anharmonicity on isotopic exchange equilibria. In: Spindel, W. (Ed.), Isotope effects in chemical processes. . Advances in Chemical Series. American Chemical Society, Washington, D.C., pp. 185–191. Yamaji, K., et al., 2001. Theoretical estimation of lithium isotopic reduced partition function ratio for lithium ions in aqueous solution. Journal of Physical Chemistry A 105, 602–613. Zeebe, R.E., 2005. Stable boron isotope fractionation between B(OH)3 and B (OH)−4 . Geochimica et Cosmochimica Acta 69, 2753–2766.