Atomistic simulation of trace element incorporation into garnets—comparison with experimental garnet-melt partitioning data

Geochimica et Cosmochimica Acta, Vol. 64, No. 9, pp. 1629 –1639, 2000 Copyright © 2000 Elsevier Science Ltd Printed in the USA. All rights reserved 0016-7037/00 $20.00 ⫹ .00

Pergamon

PII S0016-7037(00)00336-7

Atomistic simulation of trace element incorporation into garnets— comparison with experimental garnet-melt partitioning data W.

VAN

WESTRENEN,1,* N. L. ALLAN,2 J. D. BLUNDY,1 J. A. PURTON,3 and B. J. WOOD1

1

CETSEI, Department of Earth Sciences, University of Bristol, Wills Memorial Building, Bristol BS8 1RJ, UK 2 School of Chemistry, University of Bristol, Cantock’s Close, Bristol BS8 1TS, UK 3 CLRC, Daresbury Laboratory, Warrington, Cheshire, WA4 4AD, UK (Received July 21, 1999; accepted in revised form January 3, 2000)

Abstract—We have studied the energetics of trace element incorporation into pure almandine (Alm), grossular (Gros), pyrope (Py) and spessartine (Spes) garnets (X3Al2Si3O12, with X ⫽ Fe, Ca, Mg, Mn respectively), by means of computer simulations of perfect and defective lattices in the static limit. The simulations use a consistent set of interatomic potentials to describe the non-Coulombic interactions between the ions, and take explicit account of lattice relaxation associated with trace element incorporation. The calculated relaxation (strain) energies Urel are compared to those obtained using the Brice (1975) model of lattice relaxation, and the results compared to experimental garnet-melt trace element partitioning data interpreted using the same model. Simulated Urel associated with a wide range of homovalent (Ni, Mg, Co, Fe, Mn, Ca, Eu, Sr, Ba) and charge-compensated heterovalent (Sc, Lu, Yb, Ho, Gd, Eu, Nd, La, Li, Na, K, Rb) substitutions onto the garnet X-sites show a near-parabolic dependence on trace element radius, in agreement with the Brice model. From application of the Brice model we derived apparent X-site Young’s moduli EX(1⫹, 2⫹, 3⫹) and the ‘ideal’ ionic radii r0(1⫹, 2⫹, 3⫹), corresponding to the minima in plots of Urel vs. radius. For both homovalent and heterovalent substitutions r0 increases in the order Py–Alm–Spes–Gros, consistent with crystallographic data on the size of garnet X-sites and with the results of garnet-melt partitioning studies. Each end-member also shows a marked increase in both the apparent EX and r0 with increasing trace element charge (Zc). The increase in EX is consistent with values obtained by fitting to the Brice model of experimental garnet-melt partitioning data. However, the increase in r0 with increasing Zc is contrary to experimental observation. To estimate the influence of melt on the energetics of trace element incorporation, solution energies (Usol) were calculated for appropriate exchange reactions between garnet and melt, using binary and other oxides to simulate cation co-ordination environment in the melt. Usol also shows a parabolic dependence on trace element radius, with inter-garnet trends in EX and r0 similar to those found for relaxation energies. However, r0(i⫹) obtained from minima in plots of Usol vs. radius are located at markedly different positions, especially for heterovalent substitutions (i ⫽ 1, 3). For each end-member garnet, r0 now decreases with increasing Zc, consistent with experiment. Furthermore, although different assumptions for trace element environment in the melt, e.g., REE3⫹ (VI) vs. REE3⫹ (VIII), lead to parabolae with differing curvatures and minima, relative differences between end-members are always preserved. We conclude that: 1. The simulated variation in r0 and EX between garnets is largely governed by the solid phase. This stresses the overriding influence of crystal local environment on trace element partitioning. 2. Simulations suggest r0 in garnets varies with trace element charge, as experimentally observed. 3. Absolute values of r0 and EX can be influenced by the presence and structure of a coexisting melt. Thus, quantitative relations between r0, E and crystal chemistry should be derived from well-constrained systematic mineral-melt partitioning studies, and cannot be predicted from crystal-structural data alone. Copyright © 2000 Elsevier Science Ltd recently, experimental evidence (e.g., Beattie, 1993; LaTourrette et al., 1993) suggested that only garnet meets this requirement. Similarly, Sm-Nd and Lu-Hf systematics (e.g., Salters, 1996) and MORB REE patterns (Shen and Forsyth, 1995) have been explained by the presence of garnet in the MORB source, and several percent melting in the garnet peridotite stability field has been incorporated in several MORB melting models (e.g., McKenzie and O’Nions, 1991; Shen and Forsyth, 1995). However, Robinson and Wood (1998) show that melting anhydrous garnet peridotite requires unrealistically high mantle potential temperatures, and Blundy et al. (1998) and Wood et al. (1999) conclude from partitioning experiments that clinopyroxene in the spinel stability field can impose a garnet-like

1. INTRODUCTION

How trace elements partition between minerals and silicate melts is of considerable current interest. For example, garnetmelt and clinopyroxene-melt partition coefficients (D’s) play a critical role in the ongoing debate about the postulated involvement of garnet in the formation of mid-ocean ridge basalt (MORB). Specifically, U-Th disequilibria in present-day MORB (Bourdon et al., 1996) require a mineral phase in the MORB source that retains U preferentially over Th, and, until

*Author to whom correspondence should be addressed ([email protected]). 1629

1630

W. van Westrenen et al.

signature on coexisting melts, possibly obviating the need for any garnet in the MORB source. Alternatively, Hirschmann and Stolper (1996), Blichert-Toft et al. (1999) and others address the possible role of garnet pyroxenite veins (comprising up to 5% of the upper mantle) in causing the garnet signature. A proper evaluation of the roles of garnet and pyroxene requires detailed knowledge of the physicochemical principles governing the distribution of trace elements in mineral-melt systems. Considerable progress has been made in recent years (e.g., Wood and Blundy, 1997). First, significant improvements in analytical techniques have facilitated accurate experimentation involving large suites of trace elements (e.g., Hart and Dunn, 1993). Secondly, the development of a predictive crystal-melt partitioning model by Blundy and Wood (1994) provides a convenient framework for the interpretation of partitioning data. Thirdly, increases in both computer speed and program methodology have made computer simulations of the incorporation of trace elements in crystal structures a feasible option, providing an independent insight into the energetics accompanying substitution in silicates (e.g., Purton et al., 1996; Purton et al., 1997a; Purton et al., 1997b). In this paper, we report computer simulations on trace element incorporation in garnet end-members almandine (Alm, Fe3Al2Si3O12) grossular (Gros, Ca3Al2Si3O12), pyrope (Py, Mg3Al2Si3O12) and spessartine (Spes, Mn3Al2Si3O12). We compare results from our simulations with recent experimental observations, in order to obtain a detailed understanding at the atomistic level of the controlling factors in garnet-melt partitioning. 2. BACKGROUND AND THEORETICAL METHODS The model of Blundy and Wood (1994) proposes that crystal chemistry dominates the partitioning behaviour in mineral-melt systems. Wood and Blundy (1997) successfully applied the model to predict clinopyroxene-melt D’s as a function of pressure (P), temperature (T) and bulk composition (X), and we are currently finalising an analogous parameterisation for garnet (Van Westrenen et al., 1999b). The model is based on the observation of Onuma et al. (1968) that mineral-melt D’s for series of isovalent trace elements show a near-parabolic dependence on trace element radius, and invokes the lattice strain model of Brice (1975) to quantify this near-parabolic dependence. Brice (1975) relates the lattice strain energy Ustrain, associated with the insertion of a trace element with radius ri into a spherical site with radius r0, to the elasticity of the site (Young’s modulus, E) and the size misfit (ri–r0) between substituent and host cations: U strain ⫽ 4 ␲ EN A

冋

1 1 r 共r ⫺ r 0兲 2 ⫹ 共r i ⫺ r 0兲 3 2 0 i 3

册

(1)

where NA is Avogadro’s constant. The larger the size misfit, or the higher E, the higher the strain energy associated with a substitution, and the lower the affinity of the host structure for a particular trace element. Assuming crystals are far more rigid than melts, the energy associated with the exchange reaction between melt and crystal is dominated by the strain energy Ustrain in the crystal. The Blundy and Wood (1994) model thus predicts that the trace element which produces the highest strain energy upon incorporation in a mineral will be the most incompatible, and therefore have the lowest mineral-melt partition coefficient. Van Westrenen et al. (1999a) have applied the Blundy and Wood (1994) model to a series of isothermal, isobaric garnet-melt partitioning experiments along the pyrope-grossular join, involving a wide range of homovalent and heterovalent substitutions into the garnet X-site. Values of r0 derived for trivalent trace elements show a systematic dependence on garnet composition, consistent with the observation that the X-site radius increases from pyrope (⬃0.89 Å) to grossular (⬃1.03 Å). Figure 1 shows data from one of the experiments, together with curves

Fig. 1. Experimental garnet-melt partition coefficients D (Grt / Melt) for 1⫹, 2⫹ and 3⫹ trace elements entering the garnet X-site (data from Van Westrenen et al., 1999a). 1␴ errors are smaller than symbol size. Curves are weighted least-squares fits to data of the Blundy and Wood (1994) model. Note parabolic dependence of D(Grt / Melt) on trace element radius.

fitted using the Blundy and Wood (1994) model. The near-parabolic dependence of Ustrain on ri predicted by Eqn. 1 is clearly visible for divalent and trivalent trace elements. The apparent site modulus EX increases with increasing trace element charge, as observed for other minerals (e.g., Wood and Blundy, 1997). Interestingly, fitted values of r0 show a significant variation with trace element charge, with r0(3⫹) ⫽ 0.935 ⫾ 0.001 Å, r0(2⫹) ⫽ 0.99 ⫾ 0.03 Å, and r0(1⫹) ⬃ 1.2 Å (based on differences between DLi and DK), a feature not observed in clinopyroxene-silicate melt partitioning (Wood and Blundy, 1997), but recently documented in wollastonite-carbonate melt partitioning (Law et al., 2000). It is crucial that any predictive garnet-melt partitioning model can rationalise these observations (Van Westrenen et al., 1999b), and the primary aim of the simulations reported here is to suggest explanations for the trends observed in garnet-melt experiments. For the simulations a Born (ionic) model was used, assigning integral ionic charges, based on the accepted chemical valence rules, to all species. To describe non-Coulombic interactions between adjacent ions we used the consistent set of potentials listed in Table 1 for all the compounds studied. Cation-oxygen interactions were described by 2-body Buckingham potentials. Account of the oxide ion polarisability was taken using the shell model of Dick and Overhauser (1958), and a three-body O-Si-O bond bending term incorporated. Purton et al. (1996) and Purton et al. (1997a) have shown the applicability of these potentials for simulations of oxides (CaO and MgO) and silicates (olivine, orthopyroxene, clinopyroxene). We used the General Utility Lattice Program (GULP; Gale, 1997) for all simulations. Static simulations of perfect lattices give the lattice energy and crystal structure of the low-temperature garnets. In the static limit, the lattice structure is determined by the condition ⭸U/⭸Xi ⫽ 0, where U is the internal energy, and the variables {Xi} define the structure (i.e., the lattice vectors, the atomic positions in the garnet unit cell, and the oxygen shell displacements). Simulated structures for Py, Alm, Spes and Gros were used as a basis for defect energy calculations. In every computational run, one or more defects are introduced into the crystal, e.g., for homovalent substitution, one divalent cation at the X-site of a perfect garnet lattice was replaced by one trace element divalent cation. Initial, unrelaxed defect energies Udef(i) were calculated without allowing any atoms to move. The total energy of the defective system was then minimised by allowing the surrounding ions to relax to accommodate the misfit cation(s). Positions of cores and shells of atoms around the introduced defect(s) were optimised using the two-region approach (Catlow and Mackrodt, 1982). In the inner region, containing the defect(s) and typically total-

Atomistic simulation of partitioning into end-member garnets Table 1. Interatomic potential parameters and calculated lattice energies. Interaction Li⫹/O2⫺ O2⫺/O2⫺ Na⫹/O2⫺ Mg2⫹/O2⫺ Al3⫹/O2⫺ Si4⫹/O2⫺ K⫹/O2⫺ Ca2⫹/O2⫺ Sc3⫹/O2⫺ Mn2⫹/O2⫺ Fe2⫹/O2⫺ Co2⫹/O2⫺ Ni2⫹/O2⫺ Rb⫹/O2⫺ Sr2⫹/O2⫺ Ba2⫹/O2⫺ La3⫹/O2⫺ Nd3⫹/O2⫺ Eu2⫹/O2⫺ Eu3⫹/O2⫺ Gd3⫹/O2⫺ Ho3⫹/O2⫺ Yb3⫹/O2⫺ Lu3⫹/O2⫺

A (kJ mol⫺1)

␳ (Å)

25331.2 2196384 122231.1 137829 107571.5 123878 65652.1 105207 125372.6 97199 116515 143927 152688 88706.4 132667 89895 138909.5 133139.7 120462 131026.6 128981.1 130274 126356.8 129974.9

0.3476 0.1490 0.3065 0.2945 0.3118 0.3205 0.3798 0.3437 0.3312 0.3262 0.3084 0.2951 0.2882 0.3772 0.3500 0.3949 0.3651 0.3601 0.3556 0.3556 0.3551 0.3487 0.3462 0.3430

⫺2865 ⫺2533 ⫺3948 ⫺15534 ⫺13125 ⫺2171 ⫺3445 ⫺14017 ⫺3722 ⫺3851 ⫺3909 ⫺3979 ⫺1958 ⫺3248 ⫺3023 ⫺12228 ⫺12527 ⫺3258 ⫺12754 ⫺12813 ⫺13073 ⫺13262 ⫺13341

KB(kJ mol⫺1 rad⫺2)

␪0 (°)

202

109.47

Shell charge (ⱍeⱍ)

k (kJ mol⫺1 Å⫺2)

⫺2.86902

7229

O2⫺/Si4⫹/O2⫺

O2⫺

For heterovalent substituents, to ensure overall charge balance, a second trace element was inserted simultaneously into the crystal, as near as possible to the first defect. Previous work (Purton et al. 1997a) has demonstrated the importance of defect association in forsterite and diopside, and so we present results here only for associated defects. For trivalent trace cations, a Na or Li cation was placed on an adjacent X-site, or an Al cation on an adjacent Z-site (YAG-type substitution). For univalent trace cations, compensation was achieved by a trivalent cation on the X-site or an Si replacing Al on the Y-site (majorite-type substitution). This choice of compensating defect is consistent with the results of Purton et al. (1997a), who evaluated the energetics of a wide range of possible substitutions in forsterite and diopside. For trivalent trace cations, the favoured mechanism is charge-balancing with a Li on the X-site in all end-member garnets. For univalent cations, the most energy-effective charge-balancing mechanism involves a trivalent cation on the X-site, but the trivalent cation differs from end-member to end-member (Lu for Py, Yb for Alm and La for Spes and Gros). In pyrope, but not the other garnets, majorite-type substitution is comparable in energy. Although all the simulations reported in this paper are in the static limit, defect energies in this limit have been shown to be a good approximation to defect enthalpies at elevated temperatures (Taylor et al., 1997). For more details on defect energy calculations in silicates the reader is referred to Patel et al. (1991); Purton et al (1996); Purton et al. (1997a).

Ulat C (kJ mol⫺1 Å6) (kJ mol⫺1)a 0 2690 0 0 0 1028 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

1631

3. RESULTS AND DISCUSSION

3.1. Defect and Relaxation Energies Table 2 shows calculated lattice parameters for all endmember garnets, which are reproduced to within 1.6% of the experimental values. Simulated bulk moduli K are 8 –20% larger than those measured at room pressures and temperatures, but the differences are largely due to: 1. The difference in temperature between the simulations, which relate to the static limit (i.e., T ⫽ 0 in the absence of lattice vibrations), and experiment. Garnet bulk moduli increase by up to 7 GPa per 300 K decrease in temperature (e.g., Anderson and Isaak, 1995). 2. Small differences between measured and simulated lattice parameters.

Interatomic potentials of the form ␸(r) ⫽ Aexp(-r/␳) ⫺ C/r6, with r the inter-ionic distance (2-body) and ␸(␪) ⫽ KB(␪ ⫺ ␪0) with ␪ the O-Si-O bond angle (3-body). Core-shell interactions given by 0.5kx2 where x is the core-shell separation. Parameters from Purton et al., (1996; 1997a). a Lattice energies are for corresponding simple oxides, e.g., Li2O, MgO, Al2O3. REE oxides were all assigned to the same space group (bixbyite structure).

Table 3 lists initial (unrelaxed) and final (relaxed) defect energies (Udef(i) and Udef(f)) for all substitutions in the four endmember garnets considered in this study. First, for a given valence, there is a linear increase in Udef(f) with increasing ionic radius (Purton et al., 1996; Purton et al., 1997a). Secondly, for heterovalent substitutions the magnitude of the defect energies depends on the charge-balancing mechanism. Also tabulated in Table 3 is the relaxation energy Urel, which is the energy released by the garnet lattice as the ions surrounding

ing 320 –380 ions, the lattice relaxation is largest and the elastic force equations are solved explicitly (i.e., the force on all ions in this region is required to be zero). In the outer region, comprising 1200 –1700 ions, lattice relaxation was assumed to be much smaller, and estimated using the Mott-Littleton approximation (Mott and Littleton 1938). The number of ions is sufficient to ensure convergence of defect energies with region size. After convergence, final, relaxed defect energies Udef(f) were obtained.

Table 2. Comparison between calculated and experimental lattice parameters and bulk moduli. a (Å)

Bulk modulus (GPa)

Mineral

Formula

Calculated

Experimental

Calculated

Experimentalb

Pyrope Almandine Spessartine Grossular

Mg3Al2Si3O12 Fe3Al2Si3O12 Mn3Al2Si3O12 Ca3Al2Si3O12

11.281 11.386 11.529 11.874

11.459 11.531 11.612 11.845

214 208 209 194

171–179 ⬃175 172–179 166–174

a b

a

Compilation by Smyth and Bish (1988). At room P and T, range from compilation of Conrad et al. (1999).

1632

W. van Westrenen et al. Table 3. Defect, relaxation and solution energies in end-member garnets. Elements listed in order of increasing ionic radius. Pyrope

Substitution Ni Mg Co Fe Mn Ca Eu Sr Ba Sc ⫹ A1Si⫺1 Lu ⫹ A1Si⫺1 Yb ⫹ A1Si⫺1 Ho ⫹ A1Si⫺1 Gd ⫹ A1Si⫺1 Eu ⫹ A1Si⫺1 Nd ⫹ A1Si⫺1 La ⫹ A1Si⫺1 Sc ⫹ LiX⫺1 Lu ⫹ LiX⫺1 Yb ⫹ LiX⫺1 Ho ⫹ LiX⫺1 Gd ⫹ LiX⫺1 Eu ⫹ LiX⫺1 Nd ⫹ LiX⫺1 La ⫹ LiX⫺1 Sc ⫹ NaMn⫺1 Lu ⫹ NaMn⫺1 Yb ⫹ NaMn⫺1 Ho ⫹ NaMn⫺1 Gd ⫹ NaMn⫺1 Eu ⫹ NaMn⫺1 Nd ⫹ NaMn⫺1 La ⫹ NaMn⫺1

U def (i) U def (f)

Almandine Urel

Usol

U def (i) U def (f)

Urel

Spessartine Usol

Divalent defects 0.22 ⫺2.36 ⫺125 ⫺129 3.86 ⫺1.09 0.00 0.00 ⫺94.5 ⫺96.8 2.26 0.23 31.8 31.4 0.46 ⫺7.64 ⫺65.8 ⫺66.6 0.74 ⫺8.56 102 100 2.39 2.97 0.00 0.00 257 243 13.6 17.3 143 139 4.36 9.82 683 587 95.6 84.0 534 471 63.5 65.0 1202 940 262 250 1013 813 200 220 1191 923 268 223 1001 796 204 193 1931 1433 497 508 1698 1292 406 464 ⫺33.1 ⫺33.4

2224 2575 2625 2748 2936 2982 3173 3459 ⫺103 248 298 421 609 655 846 1132

Li ⫹ SiA1⫺1 ⫺1562 Na ⫹ SiA1⫺1 ⫺1247 K ⫹ SiA1⫺1 ⫺461 Rb ⫹ SiA1⫺1 ⫺41 Li ⫹ REEX⫺1 248 Na ⫹ REEX⫺1 564 K ⫹ REEX⫺1 1348 Rb ⫹ REEX⫺1 1768 Li ⫹ LaMn⫺1 Na ⫹ LaMn⫺1 K ⫹ LaMn⫺1 Rb ⫹ LaMn⫺1

1803 2174 2222 2338 2507 2547 2706 2928 ⫺382 ⫺23 25 138 303 341 496 713

421 401 403 410 429 436 467 531 279 271 273 284 307 315 350 420

211 244 253 275 314 324 369 442 163 185 193 210 246 254 296 363

⫺1971 409 ⫺1777 530 ⫺1269 807 ⫺1072 1030 ⫺23 271 166 398 665 683 859 909

165 193 521 611 184 207 525 612

2129 2453 2499 2613 2787 2829 3005 3270 ⫺319 4.10 50.7 164 338 380 557 821

⫺1673 ⫺1392 ⫺659 ⫺274 50.6 332 1071 1456

Trivalent defectsa 1693 437 198 2054 399 222 2102 398 230 2214 398 248 2379 407 283 2418 411 292 2572 433 333 2789 480 400 ⫺586 266 153 ⫺237 241 164 ⫺190 241 171 ⫺81.2 245 186 79.5 259 216 116 264 224 267 289 261 478 343 323

U def (i) U def (f) ⫺241 ⫺214 ⫺189 ⫺129

⫺257 ⫺226 ⫺197 ⫺133

349 776 763 1406

317 644 628 1105

2021 2310 2352 2454 2611 2648 2806 3044 ⫺589 ⫺301 ⫺258 ⫺157 0 37.6 196 434 ⫺349 ⫺60 ⫺18 83.4 240 278 436 674

1549 1898 1944 2053 2212 2249 2399 2608 ⫺857 ⫺520 ⫺475 ⫺370 ⫺214 ⫺179 ⫺33.0 171 ⫺689 ⫺352 ⫺307 ⫺202 ⫺47.0 ⫺11.4 134 338

Univalent defectsa ⫺2064 391 169 ⫺1808 ⫺2187 ⫺1879 487 188 ⫺1567 ⫺2014 ⫺1290 631 596 ⫺903 ⫺1543 ⫺1004 730 776 ⫺560 ⫺1361 ⫺190 241 171 ⫺157 ⫺370 ⫺10.4 342 185 83.4 ⫺203 475 596 490 747 258 662 794 570 1090 436 433 169 674 335 1338 792 1680 969

Urel

Grossular Usol U def (i) U def (f)

15.45 0.24 12.13 0.24 8.23 ⫺9.93 4.10 ⫺4.03 0.00 0.00 31.6 40.3 133 180 135 154 301 406

⫺447 ⫺426 ⫺408 ⫺362 ⫺262 323 308 833

Urel

Usol

⫺523 ⫺495 ⫺468 ⫺410 ⫺288

76.1 10.8 68.9 8.25 60.0 ⫺4.27 48.2 ⫺3.76 25.7 ⫺11.1 0.00 0.00 295 27.7 108 280 28.3 82.9 718 115 296

472 412 408 401 398 399 408 436 268 220 217 213 214 216 229 262 340 292 290 286 287 289 302 336

183 195 201 215 245 252 288 348 140 139 144 155 181 187 219 273 142 141 146 157 182 188 220 274

1851 2069 2102 2178 2299 2328 2449 2631 ⫺1073 ⫺855 ⫺821 ⫺745 ⫺624 ⫺596 ⫺474 ⫺292

1252 1577 1619 1720 1867 1901 2039 2231 ⫺1421 ⫺1111 ⫺1069 ⫺973 ⫺830 ⫺797 ⫺663 ⫺475

598 492 483 459 433 427 411 400 348 256 248 228 206 201 189 183

164 150 153 159 177 181 205 248 130 102 104 106 119 123 143 181

379 446 641 801 213 287 489 654 265 339 546 712

176 183 472 547 155 156 436 507 271 271 548 617

⫺2061 ⫺1899 ⫺1377 ⫺1121 ⫺292 ⫺130 390 646

⫺2436 ⫺2285 ⫺1863 ⫺1702 ⫺475 ⫺330 82 240

375 386 486 580 182 200 308 406

203 188 429 484 182 160 391 443

Ionic radii in Å (from Shannon, 1976); energies in kJ/mol. First element always substituting for X-site cation. X ⫽ Mg for Py, Fe for Alm, Mn for Spes, Ca for Gros; REE ⫽ Lu for Py, Yb for Alm, La for Spes and Gros. a

the defect(s) move to accommodate the misfit. Urel, equivalent to Ustrain in Eqn. 1, is defined as Udef(i)–Udef(f), and is always positive. As expected, Urel is also dependent on the chargebalancing mechanism; in what follows we shall only discuss data pertaining to the charge-balancing mechanisms which give the lowest defect and relaxation energies. Figure 2 shows the variation of Urel for univalent, divalent and trivalent trace elements with trace cation (eight-fold coordinated) radius ri (Shannon, 1976) for the four garnets. The general characteristics of these plots are similar for each end-member: for each

series of trace elements with the same charge, Urel shows a near-parabolic dependence on ri. It is important to stress that the simulations do not involve the use of any concept of ionic radius. As such, Figure 2 provides independent evidence for the applicability of the form of the Brice (1975) model (Eqn. 1) to relaxation energies in defective garnet lattices. We fitted the Brice model (Eqn. 1) to each series of relaxation energies. Fitted values for apparent X-site Young’s modulus EX and ‘ideal’ ionic radius r0 are given in Table 4, and shown by the curves in Figure 2.

Atomistic simulation of partitioning into end-member garnets

1633

Fig. 2. Simulated relaxation energies (symbols) for 1⫹, 2⫹ and 3⫹ trace elements in (a) Pyrope, (b) Almandine, (c) Spessartine, (d) Grossular. Data are from Table 3. Curves are least-squares fits (Table 4) to data of the Brice (1975) model, Eqn. 1.

four compounds (Table 2) is such that further comparison with our calculations is not possible. Simulated pyrope garnet is the least compressible, which leads to a more marked discrimination between misfit cations at the pyrope X-site, and thus the highest EX in Table 4. For the heterovalent substitutions r0 increases in the order Py ⬍ Alm ⬍ Spes ⬍ Gros, as also observed for isovalent substitution, while the apparent EX decreases. As mentioned earlier, all heterovalent results are for the associated pair of defects involving the heterovalent substituent and the charge-compensating defect (Li⫹ for 3⫹ substituents and REE3⫹ for 1⫹

Differences in crystal-chemistry between the four garnets are reflected in trends of the fitted r0 and EX. For isovalent substituents (circles in Fig. 2) the position of the minimum increases from pyrope (r0 ⬇ rMg) to grossular (r0 ⬇ rCa). This trend in r0 is compatible with crystallographic data on garnet X-site radii (e.g., Smyth and Bish, 1988) and with the variation in cell parameters for simulated garnets (Table 2). The variation in apparent EX mirrors that of the calculated bulk moduli of the four garnets, which decrease from pyrope to grossular. The experimental bulk modulus of grossular is evidently lower that that of pyrope, however the range in values reported for the

Table 4. Fits of relaxation energy data (Table 3) to the Brice equation (Eqn. 1). 1 ⫹ defects

2 ⫹ defects

3 ⫹ defects

Garnet

r0

EX

r0

EX

r0

EX

Pyrope Almandine Spessartine Grossular

0.81(30) 0.85(20) 0.90(20) 1.03(10)

190(115) 179 (94) 159 (78) 120 (43)

0.89(1) 0.92(1) 0.96(1) 1.12(1)

401(14) 361(12) 309 (9) 289(14)

0.938(3) 0.974(1) 1.027(1) 1.160(5)

748(22) 705(16) 651(12) 541(19)

r0 in Å; EX in GPa. Values in parentheses are 1 standard deviation.

1634

W. van Westrenen et al.

Table 5. Influence of trace element charge on relaxation energies in spessartine. Relaxation energy Urel a

Substituent Ni Mg Co Fe Mn Ca Eu Sr Ba

Charge 1⫹

Charge 2⫹

Charge 3⫹

302 308 316 331 368 500 694 696 956

15.5 12.1 8.2 4.1 0 31.6 133 135 301

569 552 522 493 421 279 230 231 267

Energies in kJ/mol. All elements substituting for Mn.

a

substituents). For a given garnet, the apparent EX increases with increasing charge in each garnet (e.g., for Spes, EX ranges from 159 ⫾ 78 GPa for univalent defects to 651 ⫾ 12 GPa for trivalent defects), entirely consistent with experimental observations: The apparent EX(3⫹) derived from 33 published garnet-melt partitioning experiments ranges from 419 ⫾ 65 to 870 ⫾ 174 GPa, while values for EX(2⫹) are in the range 125 ⫾ 14 to 334 ⫾ 12 GPa. r0 increases systematically with defect ionic charge, with r0(2⫹) always equal to the radius of the host cation, as found for other minerals and oxides (Purton et al., 1996; Purton et al., 1997a). The increase in r0 with increasing trace element charge is opposite to the trend observed in garnet-melt partitioning data (Van Westrenen et al., 1998 and Fig. 1). To investigate this discrepancy further, we performed two additional series of defect calculations on spessartine, assigning charges of 1⫹ and 3⫹ respectively to the divalent trace elements, without varying the short-range interatomic potentials listed in Table 1. No charge-balancing substitutions were made. Resulting defect and relaxation energies are given in Table 5. Figure 3 shows a plot of Urel versus trace element radius r. For clarity, we assigned the same radii to 1⫹ and 3⫹ cations as to the corresponding 2⫹ cation. Again, r0(1⫹) ⬍ r0(2⫹) ⬍ r0(3⫹). Since the non-Coulombic interactions are the same for each series of defects, the variations in r0 with charge in Figure 3 must be caused solely by the variation in trace element charge. These charged defects impose an extra strain on the structure, since after relaxation the positions of the oxygens nearest the defect have changed from those adopted when divalent cations occupy the X-site. The result of increasing the charge of the substituent cation is to create a positively charged defect, which exerts an additional Coulombic attraction on these neighbouring oxygen anions, in this case on the 8 oxygens around the dodecahedral spessartine X-site. Decreasing the charge of the substituent ion creates a negatively charged defect, which has the opposite effect on the neighbouring oxygens. Consequently, the structure can more readily accommodate larger 3⫹ cations (positively charged defects) and smaller 1⫹ cations (negatively charged defects), leading to the observed increase in r0 with trace element charge seen in Figures 2 and 3. It is thus possible to explain the variation in r0 with trace element charge seen in the calculated relaxation energies plotted in Figures 2a– d. However, the discrepancy between simulation and

Fig. 3. Influence of charge on relaxation energies for trace elements in spessartine (symbols). Data are from Table 5. Curves are leastsquares fits to data using Eqn. 1. (Note increase in r0 (and EX) with charge.)

experiment, which indicates r0(1⫹) ⬎ r0(2⫹) ⬎ r0(3⫹), remains. A possible reason for this is that, so far, we have looked at energies calculated purely from solid garnet properties, ignoring the presence of an additional phase (i.e., melt) in the experiments. In the next section we address the influence of the melt, via the calculation of solution energies, and provide a possible explanation for the experimental trend in r0 seen in Figure 1. 3.2. Solution Energies Relaxation energies are a function only of the properties of the solid phase. In garnet-melt systems, however, partitioning of trace elements involves: (i) Removal of the element(s) from the melt; (ii) Incorporation of the element(s) into the garnet structure by a substitution mechanism (as modeled in our simulations); and (iii) Insertion of the host cation(s) into the melt. To assess the possible influence of melt presence on our results, we have calculated trace element solution energies (Usol), making a series of assumptions as to the environment of the trace-elements in the melt. These calculations are a closer approximation to the complete cation exchange process (Purton et al., 1996). For homovalent substitutions, processes (i) to (iii) can be described by the substitution reaction (in this case for defect J2⫹ replacing Mg2⫹ in pyrope) JO(m) ⫹ Mg3Al2Si3O12(s) ⫽ JMg2Al2Si3O12(s) ⫹ MgO(m), (2) We approximate Usol (i.e., the energy associated with reaction 2) by U sol ⫽ U def(J) ⫹ U lat(MgO) ⫺ U lat(JO)

(3)

As a first approximation, then, we have assumed that the local

Atomistic simulation of partitioning into end-member garnets

1635

Fig. 4. Simulated solution energies (symbols) for 1⫹, 2⫹ and 3⫹ trace elements in (a) Pyrope, (b) Almandine, (c) Spessartine, (d) Grossular. Data are from Table 3. Curves are least-squares fits (Table 6) to data of the Brice (1975) model, Eqn. 1.

environments of Mg and J in a melt are equivalent to their environments in the corresponding (solid) simple oxides MgO and JO. Thus, in addition to the defect energies tabulated in Table 3, differences between the lattice energies Ulat of the host and trace element oxides contribute to Usol. Explicit inclusion of the heats of fusion of these oxides (Aylward and Findlay, 1994) makes little difference. For example, considering molten rather than solid JO and MgO in Eqn. 2, the solution energy for Sr at garnet X-sites changes by only about 10 kJ mol⫺1. The corresponding reactions for trivalent (J3⫹) substitutions are 1/ 2 J2O3(m) ⫹ 1/ 2 Li2O(m) ⫹ Mg3Al2Si3O12(s) ⫽ 2 MgO(m) ⫹ (JLiMg)Al2Si3O12(s)

(4)

1/ 2 J2O3(m) ⫹ 1/2 Al2O3(m) ⫹ Mg3Al2Si3O12(s) ⫽ MgO(m) ⫹ SiO2(m) ⫹ (JMg2)Al2(Si2Al)O12(s)

(5)

for charge-compensation by Li on the X-site, and Al on the Z-site respectively. Corresponding solution energies in this case are approximated by U sol ⫽ U def共J ⫹ Li) ⫹ 2 U lat(MgO) ⫺ 1/ 2 U lat(J2O3) ⫺ 1/ 2 U lat(Li2O)

(6)

U sol ⫽ U def共J ⫹ Al兲 ⫹ U lat(MgO) ⫹ U lat(SiO2) ⫺ 1/ 2 U lat(J2O3) ⫺ 1/ 2 U lat(Al2O3)

(7)

A similar set of equations can be derived for partitioning of univalent trace elements (analogous to Eqn. 4 with J ⫽ REE3⫹). Lattice energies for J2O3, J2O and JO, tabulated in Table 1, were determined using the same potentials as used for defect calculations. Calculated values of Usol, which determine the favoured solution and charge-compensation mechanisms, are given in Table 3. For example, for trivalent trace elements exchange reaction (4) involving coupled J3⫹/Li⫹ substitution is lowest in energy, as mentioned in an earlier section. Solution energies associated with the lowest-energy substitutions are plotted against trace element radius in Figure 4. Clearly, solution energies, like relaxation energies, show a near-parabolic dependence on radius. This result stresses the overriding importance of crystal local environment for trace element partitioning: even when the influence of an (admittedly oversimplified) melt is included, the parabolic dependence on radius, predicted from crystal properties only (Eqn. 1), persists. Fits of simulated Usol data to Eqn. 1 (curves in Fig. 4) are given in Table 6. Differences in r0(1⫹, 2⫹, 3⫹) between the garnets follow trends identical to those discussed above for Urel, with r0 smallest for Py and largest for Gros. Apparent EX,

1636

W. van Westrenen et al. Table 6. Fits of solution energy data (Table 3) to the Brice equation (Eqn. 1). 1 ⫹ defects Garnet

r0

2 ⫹ defects EX

r0

3 ⫹ defects EX

r0

EX

Assuming simple oxides in melt Pyrope Almandine Spessartine Grossular

0.92(21) 0.94(21) 0.98(19) 1.04(17)

185 102) 0.85(3) 337(37) 0.857(13) 184(105) 0.88(3) 345(36) 0.884(10) 180(104) 0.92(2) 353(34) 1.914 (7) 168(104) 1.00(1) 364(32) 0.981 (3) Assuming YAG-type garnet environment in melt (3 ⫹ data only)

Pyrope Almandine Spessartine Grossular

0.733(31) 0.807(16) 0.887 (7) 1.047 (1)

551(39) 559(36) 567(32) 580(23) 211(19) 219(16) 228(11) 240 (2)

r0 in Å; EX in GPa. Values in parentheses are 1 standard deviation. All values are for lowest energy configuration only.

obtained from the curvature of the solution energy vs. ri plots, are now the same within error for all garnets for each series of trace elements with the same charge. For each garnet, the apparent EX still follow the order EX(3⫹) ⬎ EX(2⫹) ⬎ EX(1⫹). In contrast, absolute values of r0(1⫹), r0(2⫹) and r0(3⫹) are markedly different from those derived from relaxation energy data (Table 4 and Fig. 5). For example, in Gros, Na is the 1⫹ trace element with the lowest Usol (Fig. 4d), while Li has the lowest Urel(1⫹) (Fig. 2d). Similarly, the 3⫹ trace element with the lowest Usol in Gros is Yb (Fig. 4d), whereas La has the lowest Urel (Fig. 2d). The latter observation is in good agreement with experimental garnet-melt partitioning data, which always show a maximum close to Yb, and not La. The size of these shifts is not the same for each mineral group; for example, both Urel and Usol for trivalent cations in the Ca sites in diopside (Purton et al., 1997a) and wollastonite (K. Law, pers. comm.) have minima close to La. Interestingly, in a given end-member garnet variations in fitted r0 with trace element charge, obtained from fitting Usol, show a different trend from that seen in Urel (Fig. 3), which we have already discussed. This is illustrated in Figure 5. r0(3⫹) is now within error of or smaller than r0(2⫹), which in turn appears smaller than r0(1⫹). This is much closer to what is observed in fits to experimental garnet-melt partitioning data (Van Westrenen et al., 1998; Fig. 1), as well as in fits to data in other mineral-melt systems (e.g., calcite-melt; Law et al., 1999). Evidently the nature of the ‘other phase’ can influence which cation is most easily accommodated into a crystal lattice. Relaxation (strain) energy alone does not take this effect into account. The success of the Brice equation (Eqn. 1), which deals with crystal strain energy only, in rationalising trends in traceelement incorporation appears to be caused by a fortunate, but only partial cancellation of energy terms. For example, if we consider again the exchange reaction (2) used to approximate isovalent garnet-melt partitioning:

while simultaneously replacing J-O short-range interactions with Mg-O short-range interactions in the ‘melt’. Very roughly, these replacement energies may cancel, so that the net (solution) energy associated with reaction (2), is close to the energy released by the garnet lattice as it accommodates misfit cation J, i.e., the relaxation energy Urel. In this fortunate, albeit partial,

JO(m) ⫹ Mg3Al2Si3O12(s) ⫽ JMg2Al2Si3O12(s) ⫹ MgO(m),

Fig. 5. Fitted values of r0 as a function of trace element charge. Urel and Usol data (circles and triangles) taken from Tables 4 and 6, experimental garnet-melt partitioning data (open squares) from Van Westrenen et al. (1999a) (Gros: experiment 13, producing Py9Gros91; Py: experiment 11, producing Py84Gros16). Dashed and solid lines show trends in simulation and experiment, respectively.

(7) the exchange process involves replacing Mg-O short-range interactions with J-O short-range interactions in the garnet,

Atomistic simulation of partitioning into end-member garnets

1637

Table 7. Solution energies assuming J3Al5O12 in melt. Solution energy Usol Trace element J* Sc Lu Yb Ho Gd Eu Nd La

Ulat(J3Al5O12)

Py

Alm

Spes

Gros

⫺56975 ⫺57511 ⫺57904 ⫺58001 ⫺58432 ⫺58731 ⫺58859 ⫺59855

161 188 193 206 228 233 257 296

151 168 172 181 198 203 223 255

138 143 145 151 163 166 180 206

128 106 105 102 101 101 105 114

Energies in kJ/mol. Charge-balanced by Li on X-site (Eqn. 8).

cancellation of terms lies the success of empirical mineral-melt partitioning models, which consider only Ustrain ⫽ Urel. This discussion, however, shows that this cancellation is only partial, as shown by the differences between the relaxation and solution energy curves, and will depend on the similarity of melt and crystal environments of the ions involved. The exchange reactions studied so far make a particularly simple assumption about the nature of the local melt environment of trace elements. We decided to study the effect different melt environments could have on Usol(3⫹) by assuming a YAG-type environment for trivalent trace elements in the melt. Exchange reactions and corresponding energies then become 1/2 Li2O(m) ⫹ 1/3 J3Al5O12(m) ⫹ Mg3Al2Si3O12(s) ⫽ 2 MgO(m) ⫹ (JLiMg)Al2Si3O12(s) ⫹ 5/6 Al2O3(m), U sol ⫽ U def(J ⫹ Li) ⫹ 2 U lat(MgO) ⫹ 5/6 U lat(Al2O3) ⫺ 1/3 U lat(J3Al5O12) ⫺ 1/2 U lat(Li2O)

(8)

1/3 J3Al5O12(m) ⫹ 1/2 Al2O3(m) ⫹ Mg3Al2Si3O12(s) ⫽ MgO(m) ⫹ SiO2(m) ⫹ (JMg2)Al2(Si2Al)O12(s) ⫹ 5/6 Al2O3(m), U sol ⫽ U def (J ⫹ Al) ⫹ U lat (MgO) ⫹ U lat (SiO2) ⫹ 1/3 U lat (Al2O3) ⫺ 1/3 U lat (J3Al5O12).

(9)

The coordination number of J3⫹ in our ‘melt component’ (J3Al5O12) is now 8, as opposed to 6 in the case of J2O3 melt components. This is obviously an extreme case, and we do not propose that trivalent trace elements have a YAG-type local environment in real silicate melts. Nevertheless, it is illustrative to assess the effect on the calculated solution energies. Lattice energies for YAG-type components, again taken from GULP simulations, and calculated solution energies are given in Table 7, and shown in Figure 6 for Py and Gros. Fits of Eqn. 1 to the data are given in Table 6. The parabolic dependence of Usol on radius is again preserved, and r0(3⫹) increases systematically from Py to Gros. However, Figure 6 also shows that r0(3⫹) and apparent EX(3⫹) are different from those assuming simple oxide melt environments. The most compatible trivalent cation in Gros changes from Yb (Fig. 4) to Gd (Fig. 6), and values of apparent EX(3⫹) (211–240 GPa) are unrealistically low compared with experimental data. Although trends in r0 and apparent EX are directly related to trends in crystal-chemistry, abso-

Fig. 6. Calculated solution energies for trivalent trace elements (J3⫹) in Pyrope and Grossular using two different assumptions about trace element melt environment as discussed in the text: J2O3 (open symbols, data from Table 3) and J3Al5O12 (filled symbols, data from Table 7). Curves are least-squares fits (Table 6) to data of the Brice (1975) model, Eqn. 1.

lute values of r0 and apparent EX therefore depend on the nature and coordination environment of the phase with which the trace element is being exchanged. We have not so far explicitly discussed the relative magnitudes of the partition coefficients for differently charged ions, since we have concentrated on the variation of D with ion size for a given charge. The experimental variation in D with charge for given size (Fig. 1) is not always reproduced in the simulations. For example, all calculated solution energies for 3⫹ ions indicate that 2⫹ ions are more soluble in garnet than 3⫹ cations, but experimentally some small 3⫹ cations have higher values of D. Experimentally, the values of the maxima in the parabolae in Onuma diagrams are very sensitive to temperature and pressure (Wood and Blundy, 1997), and all the calculations reported here are in the static limit and zero pressure. In addition, calculated solution energies reflect small differences between a number of large quantities, which differ with trace element charge since the solution mechanism changes. Calculated variations in solubility with size for a given charge would therefore be expected to be more reliable than those with charge for a given size. Both experimental data (Fig. 1) and simulated Usol data (Figs. 4, 5, 6) strongly suggest that r0 in mineral-melt systems varies considerably with trace element charge. The nature of this variation depends both on crystal structure and on the nature of the melt. The successful application of the Blundy and Wood (1994) model to trace element partitioning between clinopyroxene (Wood and Blundy, 1997), garnet (Van Westrenen et al., 1999b) and a wide range of anhydrous silicate melt compositions, suggest that the local environment of trace elements in these melts does not vary significantly. However, our simulations suggest that it is possible that relationships between mineral composition and r0 and E are not transferable to, for example, systems with hydrous or carbonatitic melts. This implies that quantitative relations between r0, E and crystal chemistry should be derived from well-constrained systematic mineral-melt partitioning studies, and cannot be predicted

1638

W. van Westrenen et al.

from crystal-structural data alone (e.g., Bottazzi et al., 1999; Tiepolo et al., 1998). 4. CONCLUSIONS

In this paper, we have demonstrated the ability of static lattice energy simulations on end-member garnets to reproduce trends in high-temperature, high-pressure garnet-melt partitioning experiments. The following conclusions can be drawn from the integration of simulations and experiments: 1. Variations in r0 and EX between different garnet end-members result from variations in crystal-chemistry. This confirms the importance of crystal local environment on trace element partitioning, and allows rationalisation of garnetmelt partitioning data using the Brice (1975) model (Eqn. 1). 2. Simulations strongly suggest r0 in garnets varies with trace element charge, as observed experimentally. 3. Absolute values of r0 and EX can be influenced by the presence and structure of a coexisting melt. This precludes the use of crystal-structural data alone to predict r0. 4. Trace element incorporation in mineral-melt systems is an exchange process and therefore solution energies (which take account of both crystal lattice relaxation and this exchange with the melt) are more relevant to partitioning studies than relaxation energies (which are based on lattice relaxation only). Further work should be directed towards high-temperature simulations that explicitly take account of lattice vibrations and cation ordering in both solid and molten phases (e.g., Purton et al., 2000). Acknowledgments—This work was supported by the EU (contract ERBFMBICT 971991), NERC grants GR3/09772 and GR9/02621, EPSRC grant GR/L31340 and Dutch VSB Foundation. JDB is grateful to the Royal Society for a University Research Fellowship. A thoughtful review by Bill Minarik, and comments by an anonymous reviewer, Thomas Acda, Paul de Munnik and Fraukje Brouwer, greatly improved the manuscript. REFERENCES Anderson O. L. and Isaak D. G. (1995) Elastic constants of mantle minerals at high temperature. In Mineral physics and crystallography—A handbook of physical constants (ed. T. J. Ahrens), AGU Reference Shelf Vol. 2, pp. 64 –97. American Geophysical Union. Aylward G. and Findlay T. (1994) SI Chemical data, 3rd ed. Wiley. Beattie P. (1993) Uranium-thorium disequilibria and partitioning on melting of garnet peridotite. Nature 363, 63– 65. Blichert-Toft J., Albare`de F., and Kornprobst J. (1999) Lu-Hf isotope systematics of garnet pyroxenites from Beni Bousera, Morocco: Implications for basalt origin. Science 283, 1303–1306. Blundy J. D. and Wood B. J. (1994) Prediction of crystal-melt partition coefficients from elastic moduli. Nature 372, 452– 454. Blundy J. D., Robinson J. A. C., and Wood B. J. (1998) Heavy REE are compatible in clinopyroxene on the spinel lherzolite solidus. Earth Plan. Sci. Lett. 160, 493–504. Bottazzi P., Tiepolo M., Vannucci R., Zanetti A., Brumm R., Foley S. F., and Oberti R. (1999) Distinct site preferences for heavy and light REE in amphibole and the prediction of Amph/LDREE Contrib. Mineral. Petrol. 137, 36 – 45. Bourdon B., Zindler A., Elliott T., and Langmuir C. H. (1996) Constraints on mantle melting at mid-ocean ridges from global 238U– 230Th disequilibrium data. Nature 384, 231–235. Brice J. C. (1975) Some thermodynamic aspects of the growth of strained crystals. J. Crystal Growth 28, 249 –253. Catlow C. R. A. and Mackrodt W. C. (1982) Theory of simulation methods for lattice and defect energy calculations in crystals. In

Computer simulation of solids (eds. C. R. A. Catlow and W. C. Mackrodt), pp. 3–20. Springer-Verlag. Conrad, P. G., Zha, C-S., Mao, H., and Hemley, R. J. (1999) The high-pressure, single-crystal elasticity of pyrope, grossular, and andradite. Am. Min. 84, 374 –383. Dick B. G. and Overhauser A. W. (1958) Theory of the dielectric constants of alkali halide crystals. Phys. Rev. 112, 90 –103. Gale J. D. (1997) GULP—a computer program for the symmetry adapted simulation of solids. Faraday Trans. 93, 629. Hart S. R. and Dunn T. (1993) Experimental cpx/melt partitioning of 24 trace elements. Contrib. Mineral. Petrol. 113, 1– 8. Hirschmann M. M. and Stolper E. M. (1996) A possible role for garnet pyroxenite in the origin of the “garnet signature” in MORB. Contrib. Mineral. Petrol. 124, 185–208. LaTourrette T. Z., Kennedy A. K., and Wasserburg G. J. (1993) Thorium-Uranium fractionation by garnet: Evidence for a deep source and rapid rise of oceanic basalts. Science 261, 739 –742. Law K. M., Blundy J. D., Wood B. J., and Ragnarsdottir V. (2000) Trace element partitioning between wollastonite and silicate-carbonate melt. Min. Mag., in press. Law K. M., Blundy J. D., Wood B. J., and Ragnarsdottir V. (1999) An investigation of the crystal chemical controls on trace element partitioning between carbonated melts and wollastonite, witherite and calcite. Eos 80, S357 (abstr.). McKenzie D. and O’Nions R. K. (1991) Partial melt distributions from inversion of rare earth element concentrations. J. Petrol. 32, 1021– 1091. Mott N. F. and Littleton M. T. (1938) Conduction in polar crystals. I. Electrolytic conduction in solid salts. Trans. Faraday Soc. 34, 485– 499. Onuma N., Higuchi H., Wakita H., and Nagasawa H. (1968) Trace element partition between two pyroxenes and the host lava. Earth Plan. Sci. Lett. 5, 47–51. Patel A., Price G. D., and Mendelssohn M. J. (1991) A computer simulation approach to modelling the structure, thermodynamics and oxygen isotope equilibria of silicates. Phys. Chem. Min. 17, 690 – 699. Purton J. A., Allan N. L., Blundy J. D., and Wasserman E. A. (1996) Isovalent trace element partitioning between minerals and melts—a computer simulation model. Geochim. Cosmochim. Acta 60, 4977– 4987. Purton J. A., Allan N. L., and Blundy J. D. (1997a) Calculated solution energies of heterovalent cations in forsterite and diopside: Implications for trace element partitioning. Geochim. Cosmochim. Acta 61, 3927–3936. Purton J. A., Allan N. L., and Blundy J. D. (1997b) Impurity cations in the bulk and the {001} surface of muscovite. J. Mat. Chem. 7, 1947–1951. Purton J. A., Blundy J. D., and Allan N. L. (2000) Computer simulation of high temperature forsterite-melt partitioning. Am. Min. In press. Robinson J. A. C. and Wood B. J. (1998) The depth of the garnet/spinel transition in fractionally melting peridotite. Earth Plan. Sci. Lett. 164, 277–284. Salters V. J. M. (1996) The generation of mid-ocean ridge basalts from the Hf and Nd isotope perspective. Earth Plan. Sci. Lett. 141, 109 –123. Shannon R. D. (1976) Revised effective ionic radii and systematic studies of interatomic distances in halides and chalcogenides. Acta Cryst. A32, 751–767. Shen Y. and Forsyth D. W. (1995) Geochemical constraints on the initial and final depth of melting beneath mid-ocean ridges. J. Geophys. Res. 100, 2211–2237. Smyth J. R. and Bish D. L. (1988) Crystal structures and cation sites of the rock-forming minerals, Allen & Unwin. Taylor M. B., Barrera G. D., Allan N. L., Barron T. H. K., and Mackrodt W. C. (1997) Free energy of formation of defects in polar solids. Faraday Discuss. 106, 377–387. Tiepolo M., Vannucci R., Zanetti A., Brumm R., Foley S. F., Bottazzi P., and Oberti R. (1998) Fine-scale structural control of REE site-preference: The case of amphibole. Min. Mag. 62A, 1517–1518.

Atomistic simulation of partitioning into end-member garnets Van Westrenen W., Blundy J. D., Purton J. A., and Wood B. J. (1998) Towards a predictive model for garnet-melt trace element partitioning: Experimental and computational results. Min. Mag. 62A, 1580 – 1581. Van Westrenen W., Blundy J. D., and Wood B. J. (1999a) Crystalchemical controls on trace element partitioning between garnet and anhydrous silicate melt. Am. Min. 84, 838 – 847. Van Westrenen W., Blundy J. D., and Wood B. J. (1999b) Prediction

1639

of garnet-silicate melt trace element partition coefficients. Eos 80, S357 (abstr.). Wood B. J. and Blundy J. D. (1997) A predictive model for rare earth element partitioning between clinopyroxene and anhydrous silicate melt. Contrib. Mineral. Petrol. 129, 166 –181. Wood B. J., Blundy J. D., and Robinson J. A. C. (1999) The role of clinopyroxene in generating U-series disequilibrium during mantle melting. Geochim. Cosmochim. Acta 63, 1612–1620.