Loading next page...
/lp/ou_press/controls-on-trace-element-distribution-in-oxides-and-silicates-WV7OfwSBgz
- Publisher
- Oxford University Press
- Copyright
- © The Author(s) 2018. Published by Oxford University Press. All rights reserved. For permissions, please e-mail: journals.permissions@oup.com
- ISSN
- 0022-3530
- eISSN
- 1460-2415
- DOI
- 10.1093/petrology/egy027
- Publisher site
- See Article on Publisher Site
Abstract
ABSTRACT Understanding and quantifying the partitioning of elements at low concentrations is important for petrology, as minor and trace elements are used as tracers of many geological processes. The lattice strain model has consequently been found to be very useful, as it links mineral/melt partition coefficients to the elastic properties of the mineral and to the difference in ionic radius between the trace element and the cation it replaces. However, this model has limitations, particularly in terms of describing crystal strain, due to the form of the equation and from the choice of a hard-sphere model. After a review of the thermodynamics of trace element incorporation into crystal sites and equilibrium between mineral phases, we present classical atomistic modelling using transferable empirical potentials. Following incorporation of one (or more) mismatching element(s), crystal strain appears strongly related to the environment of the site of exchange, with anisotropic, vacancy-rich minerals deforming more than densely-packed minerals, due to anisotropic deformation, rotation and tilting of the surrounding polyhedra. As shown by the computed cation–oxygen length in strained crystal sites, bond strain is smaller than the difference in cation radii, and varies between structures, with densely-packed minerals being less strained. Consequently, computed relaxation energies are smaller in less compressible minerals. This has implications for modelling the partitioning of trace elements, and highlights the limits of using continuum mechanics below crystal cell scale. Neither oxides nor silicates have isotropic elastic properties, and at nanoscale they do not deform like continuous material. Predictions of partitioning between mineral and melt (or fluid) remain hampered by several factors, amongst which are (i) knowledge of the thermodynamics properties of the dissolved species (in melts and fluids); (ii) the precision of estimated defect and strain energies; and (iii) the presence of solid solution in natural systems, even when limited to a few percent. The latter may have orders-of-magnitude effects on the calculated partition coefficients due to interactions between strain fields around minor cations, leading to a chemical mixing regime at low concentration different from solid solutions where phase components are found in high proportions. Partitioning between mineral phases in an assemblage is described with similar equations as for mineral/liquid equilibria. In cases where terms linked to fluid species disappear, such as incorporation with similar substitutions in two minerals, the assumption that crystal strain energy is the preponderant term during cation exchange becomes unnecessary and it is preferable to use defect energies rather than strain energies. Application and discussion of these concepts are presented for mineral/melt partitioning and for partitioning between minerals, using garnet/clinopyroxene equilibria. The advantage of atomistic modelling is that it does not rely on fitting. Agreement with experimental data shows predictive accuracy within an order of magnitude, which is poorer than what may be achieved by fitting the lattice strain model to experimental data, but validating the theoretical approach and sufficient for application to some petrology problems. Improvements of the modelling for better application to minerals like augite involve systematic work on the effect of solid solutions on the strain field around defects and imply solving ordering problems. INTRODUCTION Minor and trace elements have been used for decades as tracers of geologic processes in a broad sense, with considerable success for extremely varied materials and case studies (see for example Allègre et al., 1977; Putirka, 2008). Understanding and quantifying the trace element distribution between fluids and co-existing minerals is consequently of huge interest for petrology, not only in magmatic settings, but also for metamorphic petrology where fluids have an H2O or CO2-rich solvent, as opposed to silicate melts. There is interest in using elements in low concentration for studying metamorphic processes, and particularly in quantifying partitioning for complex minerals such as amphibole and phyllosilicates (chlorite, micas, serpentine). These have been shown to be key for the cycling of light, fluid-mobile elements as well as large-ion lithophile elements (Bebout et al., 1999, 2007, 2013; El Korh et al., 2009; Pabst et al., 2012; Smye et al., 2013; Scambelluri et al., 2015) in addition to being used as proxies for crystallization environment (Deschamps et al., 2012) and for thermobarometric purposes relying on thermodynamic modelling (e.g. Plunder et al., 2012; Dragovic et al., 2012; Baxter & Caddick, 2013). However thermodynamic models for amphibole and phyllosilicates suffer from complexities arising from numerous solid solutions, polytypes, iron speciation, and metastability, all of which apply to various extents in chlorite (Holland et al., 1998; Vidal et al., 2006; De Andrade et al., 2011), phengite (Massonne & Szpurka, 1997; Parra et al., 2002), illite (Dubacq et al., 2010; Bourdelle et al., 2013), smectite (Vidal & Dubacq, 2009; Andrieux & Petit, 2010) serpentine (Evans, 2004) and amphibole (Diener & Powell, 2012). Additional difficulties such as identifying mineralogical assemblages at equilibrium in deformed, brecciated or partially molten rocks and estimating the composition of the corresponding fluids hamper thermobarometry and reactive transport modelling. An important difference between trace and major elements is their sensitivity to reservoir effects, where small fractions of minerals in the rock may control the chemical budget of trace elements, as for Nd in monazite and Hf in zircon (Garçon et al., 2011). The reservoir effect may also induce very strong zoning in minor and trace elements in minerals showing otherwise homogeneous contents of major elements (e.g. Ginibre et al., 2002) and fractionate element ratios in the manner of Rayleigh distillation (e.g. Kohn, 2003). This sensitivity may be used to reconstruct the petrologic history of magmatic and metamorphic assemblages, with the advantage of being independent from thermodynamic modelling based on pseudo-sections and multi-equilibrium thermobarometry: the trace element signature of parageneses provides an independent record of (dis)equilibrium (e.g. Zack et al., 2002; Rubatto & Hermann, 2003). This requires quantifying the extent of trace element incorporation in complex and ubiquitous crystals such as amphibole and phyllosilicates, and their partitioning with other phases. Estimating partition coefficients between crystals and melts in the absence of measurements has been shown to be possible, at least within an order of magnitude, first through the use of Onuma diagrams (Onuma et al., 1968) and then with a simple equation by Blundy & Wood (1994), which links crystal/melt partitioning to the elastic properties of the crystal and to the size misfit between the trace element and host crystal lattice during cation exchange. This approach known as the ‘lattice strain model’ has achieved considerable success and has been refined and inverted in various ways, such as recalculating the elastic properties of crystal sites by fitting the lattice strain model to experimentally-derived trace element partition coefficients (e.g. van Westrenen & Draper, 2007; Frei et al., 2008; Pickles et al., 2016). Since the thermodynamic expression of partitioning between a mineral and a melt or a metamorphic fluid is essentially identical, the lattice strain model should provide useful constrains on partitioning in metamorphic rocks and mineral assemblages (see e.g. Pickles et al., 2016), and was tested with this aim in the present study. Here we discuss the limits of lattice strain theory after a brief review of the mathematical formulation of trace element partitioning in rocks, with and without fluids, for magmatic and metamorphic systems. From this theoretical starting point, we advance the idea that using continuum mechanics approaches is inappropriate to study trace element incorporation in minerals, due to elastic anisotropy and to the site dependency of crystal strain, consistent with previous analysis of strain fields around point defects (e.g. Eshelby, 1951; Carpenter et al., 2009). Atomistic simulations are used to study strain around trace elements incorporated in various minerals. These results illustrate why using differences in cation radii as a measure of crystal strain is inaccurate. It is further shown that the local environment of the defect-bearing sites, not only their elastic properties, plays a major role in crystal strain: ‘flexible’ environments, such as interlayer positions in phyllosilicates and in the vicinity of vacancies, are more deformed by point defects (and, therefore, have higher strain energy) than stiff densely-packed crystals. Applications to the prediction of partition coefficients for mineral/melt and in a classical garnet/clinopyroxene assemblage are shown. While agreement with experimental data is imperfect and only within an order of magnitude, the theoretical approach is validated. Further improvements are possible with more detailed studies of ordering in minerals and of the thermodynamics of the fluid species. APPROACH FOR MODELLING AND THERMODYNAMIC BASIS The approach for modelling trace element incorporation within a mineral structure is briefly reviewed below, applying general equations to examples from the literature. The aim is to show the links between minerals and fluids, the complexities arising from the numerous possible substitutions, as well as to provide rationalization and simplification methods. Throughout, trace element incorporation in mineral structures is modelled via substitution within crystallographic sites, and expressed for cations only. The presented equations and conclusions also broadly apply to anions such as halogens, which follow similar incorporation mechanisms. Noble gases and neutral species are not considered here, although they may follow some of these mechanisms e.g. the filling of large pre-existing vacancies such as the ring sites of amphibole is thought to be of major importance for noble gases (Jackson et al., 2016). Other mechanisms for trace element incorporation (e.g. adsorption on surfaces, trapping in crystals within fluid inclusions) are not investigated here, although they may also be of importance in some settings for budgeting trace element contents (see e.g. Tipper et al., 2012; Kumagai et al., 2014, respectively). Incorporation of atoms within H2O-bearing crystallographic sites, such as the interlayer sheet of clay minerals or micro-porous structures as zeolites, is also not considered here, because the properties of intra-crystalline H2O molecules are difficult to model within the framework of rigid structures and the atomistic modelling carried out below cannot reproduce them. However, there is clear evidence that the properties of water clusters in the interlayer sheet of clay minerals depend on the ionic potential of cations (e.g. Vidal & Dubacq, 2009), suggesting broadly similar trace element dependencies over charge and ionic radius, although likely different in magnitude to anhydrous crystal sites. Mineral/liquid partitioning The equation for trace element incorporation as defects within a host mineral structure, from a melt or from a fluid, is written as an exchange between the host crystal and the trace element initially in solution (McIntire, 1963; Banno & Matsui, 1973; Blundy & Wood, 1994, 2003, and references therein). For instance, Ni incorporation into forsterite (Fo) is given by: Mg2SiO4+2 Ni2+(L)⇔2 Mg2+(L)+Ni2SiO4 (1) Forsterite+2 Ni2+(L)⇔2 Mg2+(L)+Ni-doped Forsterite where (L) indicates species within a liquid phase (when used in this study, the term ‘liquid’ refers to melts and metamorphic fluids, even though the latter may be in supercritical form). The corresponding partition coefficient D is simply DNiFo/L=xNiFo/xNiL (following Morse, 2015, for the definition of the partition coefficient D and of the distribution – or exchange – coefficient KDbelow), with x the molar concentration of the element in the component ( xNiL is the Ni content of the liquid and xNiFo that of the Ni-doped forsterite). The equilibrium constant K1 of reaction 1 is: K1=(aMg2+L)2.aNi2SiO4aMg2SiO4.(aNi2+L)2 (2) where a denotes the activity of the thermodynamic component, aNi2SiO4 is the activity of the Ni-olivine end-member which increases with Ni content, aMg2SiO4 in Ni-forsterite approaching unity at low Ni contents. The use of molecular end-members for the solid phase here is for simplification of notation, it is not necessarily assumed that discrete units of each end-member are found in the mineral structure. Molecular end-members have often been preferred for silicate melts to 1M ideal solute components in publications, including thermodynamic datasets (e.g. Hirschmann & Ghiorso, 1994; Holland & Powell, 2011; Green et al., 2012). However the activity of a species may be equivalently expressed with both 1M components or molecular end-members provided that activity-concentrations relations are accounted for (Evans & Powell, 2006; Dubacq et al., 2013). The choice of end-members impacts the relationship between the equilibrium constant K1 and the corresponding distribution coefficient KD, defined as: KD=(nNi/nMg)S(nNi/nMg)L=xNi2+S.xMg2+LxMg2+S.xNi2+L (3) Within an ionic solution model where Ni2+ is incorporated on one site of multiplicity m = 2, aNi2SiO4 = (xNi2+SγNi2+)m= (xNi2+SγNi2+)2 and aMg2SiO4 = (xMg2+SγMg2+)2, where γ is an activity coefficient (Wood & Nicholls, 1978; Ganguly & Saxena, 1987; Spear, 1993; Mäder et al., 1994, amongst others). Assuming that the components behave ideally both in the liquid and in the solid (γ = 1) leads, therefore, to K1 = (KD)2. Choosing MgSi0·5O2 and NiSi0·5O2 as end-members yields: MgSi0·5O2+Ni2+(L)⇔ Mg2+(L)+NiSi0·5O2 (4) with equilibrium constant K4=aMg2+L.aNiSi0·5O2aMgSi0·5O2.aNi2+L (5) where K4 = KD in an ideal mixing scheme. In natural melts or metamorphic fluids, deviation from ideal behaviour may be intuited to be larger for the host cation Mg2+(L) than for the trace element Ni2+(L) if the latter is very dilute (Henry’s law). Dependency of partitioning on composition and solid solutions There is a wealth of discussion of the impact of composition on trace element partitioning, already detailed in the early review of McIntire (1963). Solid solutions lead to the exchange of cations over large proportions (≥ 1%) in the host crystal, such as incorporation of a fayalite component (Fa: Fe2SiO4) into olivine [similar to eq. (1)]. These solid solutions are linked to a change in the mineral/liquid partition coefficient, which may range over orders of magnitude even in the case of ideal mixing (e.g. Lagache & Dujon, 1987). For the olivine example, this translates as: DNiolivine/L varies because DNiFo/L≠DNiFa/L. An exchange reaction in the form of equation (1) may be written for the Fa component with equilibrium constant K1(Fa)=(aFe2+L)2.aNi2SiO4aFe2SiO4.(aNi2+L)2. The variation of the equilibrium constant with composition is linear in an ideal model; this is not so in the case of Ni incorporation in olivine as DNiolivine/L varies non linearly with the Fo content (Schreiber, 1979; Matzen et al., 2013): even at low Ni contents where it may be argued that the activity of trace species nears ideality (e.g. Häkli & Wright, 1967), aMg2SiO4 and aMg2+L may not be equal to xMg2SiO4 and xMg2+L, so that K1 is affected. For minor Ni contents (0·5 mol. %), Hirschmann & Ghiorso (1994) showed that aNiSi0·5O2 is not equal to xNiSi0·5O2 and varies along the forsterite-fayalite solid solution. Substitutions over several sites For minerals with several crystallographic sites prone to exchange, such as forsterite which has two octahedral sites (M1 and M2) where divalent cations may be hosted, trace element partitioning may differ between sites, in the manner of Fe–Mg ordering reactions (e.g. Rinaldi et al., 2000; Redfern et al., 2000). In forsterite the corresponding exchange reaction is: (Ni2+)M1+(Mg2+)M2⇔(Mg2+)M1+(Ni2+)M2 (6) with equilibrium constant K6=aMg2+M1.aNi2+M2aNi2+M1.aMg2+M2, which may be similarly expressed in terms of the molecular end-members NiM1MgM2SiO4 and MgM1NiM2SiO4. Two reactions in the form of equation (1) appear: Mg2SiO4+Ni2+(L)⇔Mg2+(L)+NiM1MgM2SiO4 (7) Mg2SiO4+Ni2+(L)⇔Mg2+(L)+MgM1NiM2SiO4 (8) where only two out of the three reactions [eq. (6)–(8)] are independent. For forsterite, small divalent cations (such as Co2+ and Ni2+) have been shown to favour the M1 site, and larger cations (such as Mn2+) the M2 site, with the consequence that K6 ≈ 10 (Rajamani et al., 1975; Annersten et al., 1984; Miyake et al., 1987). This approach is useful for rationalizing exchange reactions in complex minerals like phyllosilicates: the Ni2+-Mg2+ exchange in clinochlore (Mg5Al2Si3O10(OH)8) may take place in up to four octahedral sites occupied by Mg (Joswig et al., 1980; Zanazzi et al., 2006), with varying bond types (4 Mg-O and 2 Mg-OH bonds in M1 and M2, 6 Mg-OH bonds in M3 and M4) and bond lengths (from ∼1·98 Å to ∼2·10 Å) and, therefore, varying bond enthalpy (e.g. Vieillard, 1994a). Coupled substitutions and vacancies Coupled substitutions (e.g. Tschermak) as well as creation of vacancies during exchange allow maintaining charge balance, with contrasted efficiency as discussed below. Coupled trace element exchange: Li+ + La3+ in forsterite Heterovalent substitutions may be achieved by simultaneous incorporation of two defects, with identical or different charges. For example trivalent REE may be exchanged with Mg2+ in forsterite concomitantly to substitution of another Mg2+ for Li+: Mg2SiO4+Li+(L)+La3+(L)⇔2 Mg2+(L)+LiLaSiO4 (9) with equilibrium constant K9=(aMg2+L)2.aLiLaSiO4aMg2SiO4.aLi+L.aLa3+L. This reaction is energetically favourable to straight La3+ incorporation into the solid (e.g. Allan et al., 2003), however its efficiency (as measured via the activity of the trace element mineral end-member LiLaSiO4) may be expected to be small unless large proportions of Li are available in the liquid phase (aLiLaSiO4 is proportional to the product of trace element activities: aLi+L.aLa3+L). This leads McIntire (1963) to state ‘it is reasonable to conclude that two or more trace elements will partition between a liquid and a solid solution as if the other trace elements were not present’ (p.1220). Defect and substitution with a major element: Ba2+ in muscovite The incorporation of heterovalent trace elements is eased by the presence of a major element in the substitution mechanism. For example Ba2+ incorporation in muscovite (KAl3Si3O10(OH)2) can take place through the reaction: Ba2+L+Al3+L+KXIIAl2,□VI(Al,Si3)IVO10(OH)2⇔K+L+Si4+L+BaXIIAl2,□VI(Al2,Si2)IVO10(OH)2 (10) where □ is a vacancy and IV, VI and XII are coordination numbers (Grapes, 1993; Hetherington et al., 2003; Ibhi et al., 2005). Ba2+ is incorporated in an interlayer position replacing K+, and Al3+ is exchanged with Si4+ in an adjacent tetrahedral site. In metamorphic fluids at equilibrium with muscovite, Al3+ is readily available [e.g. Verlaguet et al., 2006), therefore, not limiting unlike Li in eq. (9)]. Vacancy creation: Ba2+ in muscovite and La3+ in forsterite Another possibility of Ba2+ incorporation into muscovite is the creation of interlayer vacancies □XII, which by definition are not limiting either: Ba2+L+2 KXIIAl2,□VI(Al,Si3)IVO10(OH)2⇔2 K+L+2 (Ba0·5,□0·5)XIIAl2,□VI(Al,Si3)IVO10(OH)2 (11) The creation of vacancies may account for incorporation of many cations through heterovalent substitutions, such as trivalent cations like La3+ into forsterite (here expressed as Mg4Si2O8 for convenience): Mg4Si2O8+2 La3+(L)⇔3 Mg2+(L)+La2□MgSi2O8 (12) The prevalence of one incorporation mechanism over the other [eq. (10) and (11)] may be concentration-dependent, with the possibility that one of the two substitutions is energetically favoured in the dilute limit but not at higher concentrations, involving short- and long-range ordering and favouring the other substitution. Mineral/mineral partitioning Even though fluids are major actors in redistributing species within metamorphic rocks of all grades, remnants of fluid phases are often elusive and harder to characterize than the minerals with which they equilibrated. Expressing trace element partitioning between minerals is, therefore, important to reconstruct the trace element budget of metamorphic rocks and restite, and to verify whether equilibrium was achieved in the considered assemblage. Taking the example of Ni and Mn in forsterite MgSi0·5O2 and enstatite (En, MgSiO3 orthopyroxene), Ni incorporation in enstatite from a liquid L is in the form of equation (4): MgSiO3+Ni2+(L)⇔ Mg2+(L)+NiSiO3 (13)Equations (4) and (13) combine to yield: MgSi0·5O2+NiSiO3⇔ NiSi0·5O2+MgSiO3 (14) which describes equilibrium between minerals, via a solvent (melt or metamorphic fluid) now implicit. Note that the terms for liquid species disappear only if the incorporation mechanism in the mineral phases involves the same host cation. Mn incorporation leads similarly to: MgSi0·5O2+MnSiO3⇔ MnSi0·5O2+MgSiO3 (15)Equations (14) and (15) yield: MnSi0·5O2+NiSiO3⇔ NiSi0·5O2+MnSiO3 (16) which is the final partition equation for both trace elements Mn and Ni in enstatite and forsterite. The corresponding distribution coefficient KD(MnNi)Fo/En is: KD(MnNi)Fo/En=xMnFo.xNiEnxNiFo.xMnEn=DMnFo/EnDNiFo/En (17) where xMnFo is the concentration of Mn in forsterite and so on. DMnFo/En and DNiFo/En are individual partition coefficients for Mn and Ni, which if known can be used to calculate KD(MnNi)Fo/En. Assuming that the considered trace elements behave ideally (i.e. non-ideal interactions are sufficiently small that activity equals concentration), KD(MnNi)Fo/En is the equilibrium constant of reaction 16 and can be directly measured on the mineral assemblage (an assumption already described by Berthelot, 1872: ‘Deux corps étant mis en présence simultanément de deux dissolvants se partagent entre eux comme si chacun de ces corps agissait isolément’, p. 417). The interest of this approach is obvious for petrology (as noted by Häkli & Wright, 1967), from melting models to lower-temperature processes such as using tracers of equilibrium in metamorphic rocks and quantifying the partitioning of pollutants. However, accurate estimation of partition coefficients between minerals is not straightforward for reasons of precision of measurements and assumptions over equilibrium (as discussed by Zack et al., 2002, for eclogitic conditions). A DISCUSSION OF THE LATTICE STRAIN MODEL Role of strain energy Distinction between the free enthalpy of fusion of the host mineral [Mg2SiO4 in eq. (1)] and the free energy of exchange (e.g. ΔGr of reaction 4) is fundamental to the description of partitioning (Blundy & Wood, 1994), through the equilibrium relation ΔGr04=-R.T.ln(K4) where T is temperature and R the ideal gas constant. The free energy of exchange ΔGr04 is assumed to be controlled by the energy required to deform the forsterite lattice when incorporating Ni as a defect into one of the Mg sites of forsterite. This translates as assuming that the difference in standard state chemical potentials of the liquid species is equal to the difference in standard state chemical potentials of their solid counterpart plus the strain energy of the host crystal UεNi→Fo, or in terms of equation (4): μMg2+(L)0-μNi2+(L)0=μMgSi0·5O20-μNiSi0·5O20+UεNi→Fo (18) Assuming first that the activity of components equals their concentration (e.g. aMgFo=xMgFo), then that UεNi→Fo is the main control on partitioning, combines with equation (5) to yield: DNiFo/L=xMgFoxMgLexp(-UεNi→FoR.T) (19) with xMgFo the fraction of Mg in forsterite. Concentrations of host and trace cations may be expressed in weight or mole fractions in the bulk mineral or liquid, but not as fractions of crystallographic site occupancies, to account for the multiplicity of sites in minerals. In this expression, the composition and properties of the liquid appear only through its cation concentration (here xMgL). This implies direct limitations, for example on modelling calcium partitioning in olivine (more sensitive to the amount of alumina, alkali and ferrous iron in the melt than to temperature, e.g. Libourel, 1999) and on modelling the effect of pressure on the melt structure (which impacts speciation and partitioning: Mysen & Kushiro, 1979; Cochain et al., 2015, amongst others). It is also noteworthy that strain energy ( UεNi→Fo) is by definition a positive quantity, whether the defect is smaller than the host ion, larger, or a vacancy. Therefore, in equation (19), DNiFo/L increases with temperature (as the term exp(-UεNi→FoR.T) approaches unity) up to xMgFoxMgL. As values larger than xMgFoxMgL have been reported for DNiFo/L above 1000 °C (e.g. Häkli & Wright, 1967; Duke, 1976; Mysen, 1976), it is arguable that the first assumption of ideality of the host cation may be erroneous in some systems. For trace element partitioning between two minerals (taking eq. (14) as an example), similar reasoning leads to the disappearance of the terms involving liquid species. For example the Gibbs energy of reaction for Ni incorporation into enstatite [eq. (13)] is assumed to follow ΔGr013=UεNi→En, and Ni partitioning between forsterite and enstatite [eq. (14)] simplifies as ΔGr014=UεNi→Fo-UεNi→En. In this view: DNiFo/En=xMgFoxMgEnexp( -UεNi→Fo+UεNi→EnR.T) (20) which is independent of the liquid species terms but still relies on assuming that the difference in properties of melts and mineral species is subordinate to the energy of crystal strain [an assumption which yielded eq. (19) from eq. (18)]. The lattice strain model Following the work of Brice (1975), the model of Blundy & Wood (1994) assumes that the size difference between cations leads to a local strain in the crystal, with a strain energy depending on the Young’s modulus of the crystal, equivalent to a ‘local elastic pressure’ (Brice, 1975). The following equation has been proposed for strain energies Uε (Brice, 1975; Blundy & Wood, 1994, 2003; Wood & Blundy, 1997, amongst others): Uε=4.π.E.NAr02(rM-r0)2+13(rM-r0)3 (21) where E is the Young’s modulus of the crystal, rM the ionic radius of the defect element M, r0 the ‘strain-free’ radius of the host cation in the crystal lattice, NA Avogadro’s constant. Brice (1975) originally calculated the strain energy [eq. (21)] resulting from replacing the host ion by an ion with similar charge, but different radius, as the work required to expand a sphere of radius r0 to rM. Injecting this expression into a partitioning equation such as equation (19) allows estimating partition coefficients using only the charge and radius of the trace cations. This expression of strain energy has been found to partly match relaxation energies. Relaxation energies are defined as the difference between the defect energy before and after relaxation, e.g. for Ni in forsterite the difference between the energetic level of the forsterite structure with one Ni atom replacing one Mg atom, and that of the altered forsterite structure with atomic positions optimized around the Ni defect to minimize the overall energy. Such relaxation energies have been calculated with semi-empirical atomistic models (Allan et al., 2003; van Westrenen et al., 2003): the shapes of the dependence of the strain energy [eq. (21)] and relaxation energies to the cation radius appear globally identical, but estimated Young’s moduli are higher by an order of magnitude than their measured equivalents for forsterite and diopside for example, and must be altered as well as the ionic radius of the host cation to improve the fit. Differences between partition coefficients measured and estimated with equation (21) have been explained by variations in optimum site radius and site valence (e.g. Wood & Blundy, 2001; Adam & Green, 2006), to account for solid solutions and sites with mixed occupancies. Hammouda et al. (2009) have shown systematic differences between calculated and experimentally-derived clinopyroxene/melt and garnet/melt partition coefficients, which they attributed to a combination of vacancies in clinopyroxene and of the peculiarities of the structure of carbonated melts. For minerals with several crystal sites, the assumption that the elastic properties of the bulk crystal are relevant to describe interactions on varied sites has been mitigated with the concept of ‘site compressibility’ where using an elastic property per site has been preferred (e.g. Allan et al., 2003; Frei et al., 2008; Cartier et al., 2014), with Young’s modulus or bulk modulus depending on publications. The incorporation of ions in sites where the host cation has a different charge, involving a coupled substitution mechanism as presented above to maintain charge balance, has also been studied in the light of the lattice strain model (see Wood & Blundy, 2001), where the apparent Young’s modulus has been identified as a fundamental control on the curvature of the partition equation when plotted in Onuma diagrams. The misfit between observables and equation (21) is, therefore, often explained by the properties of the crystal, of the exchanged cations or of the melts, however these discrepancies likely originate from the form of equation (21), which is questionable. Brice (1975) assumed that the strain integrates over the sphere radius as ∫rM∞ε dr=rM-r0. This is equivalent to the sum of the displacements of all particles within the strained material but inconsistent with the definition of strain, which is expressed as a dimensionless quantity, scaled to the original length (or volume) of the material: ε=rM-r0r0 or ε=VM-V0V0, with the consequence that strain such as that due to displacement of atoms around the exchanged cation tends to zero with increasing size of the material, unlike in the definition of Brice (1975). The use of Young’s modulus as the appropriate elastic property to solve a radial strain problem is also questionable. Assuming that cation exchange takes place at the origin of an isotropic body of infinite size, and forcing a larger cation defined as an elastic sphere of radius (1+α)r0 into a spherical hole of radius r0, the elastic displacement ui along i is defined as ui=δxi/r3, with δ=αr03(1+ν)/3(1-ν) (see equation (20) of Eshelby, 1951) where ν is Poisson’s ratio. ν depends on the shear modulus µS and on the bulk modulus BM of the material such as ν=(3BM-2 μS)/(2(3BM+μS)). This is illustrated in Figure 1 using a square grid deforming isotropically by insertion of a larger sphere within the material: lines are deformed into curves, showing that shearing takes place even if the displacement is radial and Young’s modulus is the same in all directions. Two out of the four elastic parameters (E, µS, BM and ν) are, therefore, necessary for complete description of strain in an isotropic mineral. Corresponding equations for traction forces within inclusions and heterogeneities are given in Eshelby (1951, 1956), with the implication that forces are not identical whether α is positive (the included sphere is larger than r0) or negative (the included sphere is smaller). It follows that the strain energy is a strongly asymmetrical function of the size difference between host cation and defect. This is important when considering substitution of a series of defects larger or smaller than the host cation, such as substitution of divalent defects with Ca2+, where, for example, Ni2+ and Mg2+ are smaller, and Ba2+ larger. In other words, strain energy into isotropic material is expected to vary between substitution of Mg2+ into Ca-rich material, and substitution of Ca2+ into Mg-rich material, as quantified below for oxides. Fig. 1. View largeDownload slide Schematic effect of homovalent cation substitution (An+ ⇔ Bn+) within an isotropic material symbolised by a square grid where atoms are located on the nodes. The initial grid around A (solid lines) is homogeneously radially deformed to the dotted coloured curves when the bigger B cation is substituted for A (white arrow). Black arrows (1–4) show the displacement of the nodes of the grid and decrease in length away from the centre (1 > 2 > 3 > 4). Shearing is highlighted by the curving of the originally straight grid, a consequence of radial strain. Fig. 1. View largeDownload slide Schematic effect of homovalent cation substitution (An+ ⇔ Bn+) within an isotropic material symbolised by a square grid where atoms are located on the nodes. The initial grid around A (solid lines) is homogeneously radially deformed to the dotted coloured curves when the bigger B cation is substituted for A (white arrow). Black arrows (1–4) show the displacement of the nodes of the grid and decrease in length away from the centre (1 > 2 > 3 > 4). Shearing is highlighted by the curving of the originally straight grid, a consequence of radial strain. Whether defects in minerals (and particularly oxides and silicates) may be accurately accounted for with an isotropic model is discussed below, as well as less obvious assumptions in the use of the ionic radii derived by Shannon (1976) in equation (21) when estimating strain due to cation exchange. These assumptions are linked to the use of a hard-sphere model with fixed size for cations and oxygen. First, atomic radii (and, therefore, bond lengths) of defects are assumed to be the same in the strained crystal and their oxide, second the values used for cation radii must be accurate both in absolute value and relative to each other, as equation (21) depends exponentially on them. A simple consequence is that predicted defect incorporation will be energetically favourable to ‘soft’ crystals or sites (low Young’s modulus) at the expense of less compressible crystals or sites (high Young’s modulus), regardless of the effect of strain, when it is arguable that ‘soft’ crystals are more strained than less compressible ones. Any of these assumptions limits the use of the lattice strain model. They are tested below using atomistic modelling and in regard of continuum mechanics, with the conclusion that any equation based on a single elastic parameter for estimating partitioning is simplistic at best, even for minerals as simple as periclase (cubic MgO). ATOMISTIC MODELLING: LINKING STRAIN, DEFECT ENERGIES AND RELAXATION ENERGIES In this section, strain energies are evaluated using atomistic modelling based on force field methods with the GULP code (Gale, 1997; Gale & Rohl, 2003) and empirically-derived potentials. GULP is used to find a local minimum to the lattice energy of the mineral by varying atomic positions within the crystal configuration, a method well adapted to the study of structural defects (see in particular Allan et al., 2003; van Westrenen et al., 2003), ordering, and solid solution in silicates (Bosenick et al., 2001; Vinograd & Sluiter, 2006, amongst many others), including phyllosilicates (Sainz-Diaz et al., 2001a, 2001b; Palin et al., 2001; Palin & Dove, 2004; Dubacq et al., 2011). Defect energies and relaxation energies The strain energy Uε is estimated from the relaxation energy (the difference in defect energy before and after relaxation) in previously-optimized crystal structures. The defect energy is first calculated as the difference between the optimized defect-free structure and the same structure with the defect. Relaxation consists in re-optimizing the defect-bearing structure to obtain the relaxed defect energy. Therefore, in exchange reactions, the relaxed defect energy, when expressed per mole of defect, may be related to the difference between the standard state chemical potential of the host end-member and that of the defect-bearing end-member. For Ni incorporation into olivine [eq. (4)], this yields μNiSi0·5O20=μMgSi0·5O20+UrelaxedNi→Fo; with UrelaxedNi→Fo the relaxed defect energy. It is noteworthy that this definition of the standard state chemical potential of the NiSi0·5O2 species is that of a Ni-doped forsterite end-member, not of liebenbergite (Ni2SiO4) which has a different structure, different bond lengths and a different standard state chemical potential. The present definition of μNiSi0·5O20 is only valid when Ni is a trace element in forsterite. Simple access to defect energies is an advantage of atomistic methods as the use of defect energy is superior to that of relaxation energies in some cases, and particularly for mineral partitioning where the assumption of strain energy being the main control becomes unnecessary. For Ni incorporation into olivine: ΔGr04=μMg2+(L)0-μNi2+(L)0+UrelaxedNi→Fo (22) which is of limited use for mineral/melt partitioning as the μ(L)0of the melt species are only poorly known, however yields an expression in the form of equation (20) for Ni partitioning between forsterite and enstatite: DNiFo/En=xMgFoxMgEnexp(-UrelaxedNi→Fo+UrelaxedNi→EnR.T) (23) which this time is independent of the properties of the melt species and does not rely on the assumption that strain energy is preponderant. Computational method Similarly to Allan et al. (2003), a self-consistent collection of Buckingham potentials derived from previous studies (Lewis & Catlow, 1985; Purton et al., 1996; van Westrenen et al., 2000; Sainz-Diaz et al., 2001a) was used assigning integral positive charges Q to cations (others than H+, such as QFe2+ = 2, QFe3+ = 3 and so on), a core/shell model for oxygen (QOcore = 0.86902, QOshell = -2.86902 with a spring constant) and partial charges for hydroxyl groups (QOcore = -1.426, QH+ = 0.426). A Morse potential is used for hydroxyl groups (see Palin et al., 2001). Values for potentials are given in Supplementary Data Electronic Appendix 1; supplementary data are available for downloading at http://www.petrology.oxfordjournals.org, with values for three-body interactions used to ensure correct cation coordination. Defect energies are calculated with the Mott-Littleton method (Mott & Littleton, 1938) with inner regions between 12 and 18 Å and outer regions between 20 and 35 Å depending on structure complexity. When assigning defects to sites with partial occupancies (such as Aliv and Siiv in the tetrahedral sheet of muscovite), cations were individually positioned within super-cells of appropriate size (for muscovite, respecting Aliv avoidance). The mean field approximation has been used for these sites when defects were assigned to other sites. This approach is efficient for a local (in terms of energy) optimization of a given crystal structure and well adapted here where the initial configuration is known. Computed optimized structures compare very favourably with more sophisticated modelling: within a few percent of the calculated value at worst for bond lengths in both unstrained and strained crystals (see for example ab initio studies of Freyria-Fava et al., 1993; Orlando et al., 1994; Israel et al., 2003, on MgO, ScO, BaO, CaO), which is sufficient for the discussion in the present study. In terms of evaluating defect energies, comparison between for example the ab initio work of Orlando et al. (1994) and the present approach gives a discrepancy of ∼0·15 eV for a Ca2+ defect in periclase MgO, which is about 2% of the calculated value. Comparison with more complicated structures and chemical systems is difficult because of the sparsity of ab initio data. A 15% uncertainty on defect energies would leave the discussion in this section unaffected, but estimation of partition coefficients asks for relative precision better than 5%, as defect energies vary little within similar structures (see garnet-clinopyroxene example below). A number of published room-temperature structure refinements have been selected to compare relaxation energies, starting with simple oxides with rocksalt structure (SrO, BaO, CoO, NiO, FeO, EuO, MgO, CaO, MnO: Wyckoff, 1963; Hazen, 1976; Pacalo & Graham, 1991; Fiquet et al., 1999); α-corundum (Lewis et al., 1982); eskolaite (Finger & Hazen, 1980) and Sc2O3 (Schleid & Meyer, 1989). The following room-temperature structure refinements for silicates have been used, simplified when necessary: actinolite (Evans & Yang, 1998, using their refinement with the highest Mg proportion and neglecting minor Fe2+, Mn2+, Ti4+ and Al3+); almandine and pyrope (Armbruster et al., 1992); spessartine (Rodehorst et al., 2002); magnesiochloritoid (Ivaldi et al., 1988, replacing Fe2+ with Mg2+ and assigning H+ bonded to O1A symmetrically to H+ bonded to O1B); clinochlore (Welch & Marshall, 2001); diopside (Levien & Prewitt, 1981); enstatite (Hugh-Jones & Angel, 1994); fayalite (Kudoh & Takeda, 1986); forsterite (Birle et al., 1968); grossular (Meagher, 1975); jadeite (Prewitt & Burnham, 1966); kushiroite (Okamura et al., 1974); lizardite (Mellini & Zanazzi, 1987); muscovite (Richardson & Richardson, 1982); phlogopite (Redhammer & Roth, 2002); talc (Diego et al., 2012, measured at 223 K, setting α and γ angles at 90°). Mg-celadonite (the phengite end-member as defined by Vidal & Parra, 2000) has been generated from the phlogopite structure, replacing ivAl with ivSi, M1Mg with a vacancy M1□, and filling the M2Mg2+ site with 0·5 M2Mg2+ and 0·5 M2Al3+ for stoichiometry and following Al-avoidance. Results: Homovalent substitutions Results are first presented in terms of structural aspects of lattice strain after incorporation of a Ca2+ defect within the Mg site of periclase MgO (Fig. 2) and forsterite (Fig. 3). In these figures the final lattice configuration is shown, indicating both the distance and direction over which each cation is displaced with reference to the initial lattice (before defect incorporation), by comparison between initial and final coordinates. Strain is considered as terminated and the path followed by each atom between initial and final lattice configuration is not studied here. The displacement of atoms has been linearly interpolated between atoms after each simulation to help figure readability, without altering the results on the nodes (where atoms are located). Fig. 2. View largeDownload slide Displacement in periclase MgO around a Ca2+ defect (grey sphere) along (a) the <001> plane and (b) the <-110 > plane. Planes are shown in pink on the grey cube by the side of each drawing. Arrows show the direction of displacement, with fixed size for readability. Colours show the length of displacement on a log scale from yellow (largest displacement) to blue (negligible and within rounding on coordinates) as on the colour bar [same for (a) and (b)]. Displacement lengths have been linearly interpolated between atoms. Final atom positions are shown as yellow and black dots for Mg and O. Note that neither the direction nor the length of displacement are radial. (c) Displacement on a log scale versus initial distance of all atoms from the defect. Points overlap due to symmetry (e.g. the six closest oxygen atoms around 2·1 Å from the defect all have a calculated displacement of about 0·1 Å). Fig. 2. View largeDownload slide Displacement in periclase MgO around a Ca2+ defect (grey sphere) along (a) the <001> plane and (b) the <-110 > plane. Planes are shown in pink on the grey cube by the side of each drawing. Arrows show the direction of displacement, with fixed size for readability. Colours show the length of displacement on a log scale from yellow (largest displacement) to blue (negligible and within rounding on coordinates) as on the colour bar [same for (a) and (b)]. Displacement lengths have been linearly interpolated between atoms. Final atom positions are shown as yellow and black dots for Mg and O. Note that neither the direction nor the length of displacement are radial. (c) Displacement on a log scale versus initial distance of all atoms from the defect. Points overlap due to symmetry (e.g. the six closest oxygen atoms around 2·1 Å from the defect all have a calculated displacement of about 0·1 Å). Fig. 3. View largeDownload slide (a) Displacement in forsterite around a Ca2+ defect (grey sphere) along the <001> plane. As in Figure 2, arrows indicate the direction of displacement with fixed size; colours show the length of displacement on a log scale as on the colour bar. Displacements lengths have been linearly interpolated between atoms. (b) Bond strain (as defined in the text) in planar view around the same defect for atoms around the <001> plane. Shades of red indicate stretching and blue compression. The first octahedron (around the defect) is shown. Atomic positions are shown as dots in (a) and (b). No O and Si atoms are exactly on the <001> plane passing through the defect and, therefore, O and Si appear on (b), not (a). Fig. 3. View largeDownload slide (a) Displacement in forsterite around a Ca2+ defect (grey sphere) along the <001> plane. As in Figure 2, arrows indicate the direction of displacement with fixed size; colours show the length of displacement on a log scale as on the colour bar. Displacements lengths have been linearly interpolated between atoms. (b) Bond strain (as defined in the text) in planar view around the same defect for atoms around the <001> plane. Shades of red indicate stretching and blue compression. The first octahedron (around the defect) is shown. Atomic positions are shown as dots in (a) and (b). No O and Si atoms are exactly on the <001> plane passing through the defect and, therefore, O and Si appear on (b), not (a). Periclase MgO has a cubic NaCl rocksalt structure (space group Fm-3 m). In Figure 2, results of Ca2+ incorporation are presented with two cross-sections through the Ca2+ defect along different planes, and the length of displacement of each atom is plotted with respect to the initial distance of the atom to the defect. The calculated displacement is maximal for the oxygen atoms closest to the defect, then displacement decreases very rapidly with distance to the defect (Fig. 2c). The calculated displacement falls below 0.01 Å within 6 Å from the defect. Interestingly, the displacement is anisotropic as well as the directions of displacement of atoms, which are not all pointing away from the defect, even considering only atoms within 6 Å from the defect. Subsequently there are differences between Figure 2a and b, with less displacement along the axes. The same patterns are observed with similar implications for the other cations and in the other simple oxides investigated here, with differences in the magnitude and length of displacement as a function of defect and mineral nature. Similar calculations have been carried out for forsterite, which has a considerably more complicated structure (orthorombic–dipyramidal, space group Pbnm). A Ca2+ ion has been exchanged with Mg2+ in the M1 site. Figure 3a shows that the calculated displacement is strongly anisotropic, and that the Mg atoms closest to the defect have been shifted not away but towards the defect, even though Ca2+ is larger than Mg2+. Note that this displacement is small (on the order of 10–2 Å) compared to the displacement of the closest oxygen atoms (on the order of 10–1 Å). Figure 3b shows the strain of individual bonds linking cations and oxygen around the plane shown in Figure 3a. The Ca2+-oxygen bonds in the central octahedron show the largest amount of stretching (∼9% of the original bond length). The Ca2+-oxygen bond lengths do not reach the values observed in monticellite MgCaSiO4: the smallest and longest Ca2+-oxygen bond lengths are about 2.29 Å and 2·5 Å in monticellite (after Sharp et al., 1987), they are computed at 2·22 Å and 2·35 Å in the Ca-bearing forsterite. Oxygen-cation bonds in the surrounding Si4+ tetrahedra and Mg2+ octahedra are variably deformed, with individual polyhedra being variably strained. Stretching some bonds and shortening others induces tilting of the polyhedra close to the defect-bearing site, with the consequence that Mg2+ atoms move towards the defect as seen on Figure 3a. The incorporation of other divalent defects in the M1 site of the forsterite lattice has been computed to investigate the evolution of strain intensity with defect size, as shown in Figure 4. Polyhedra sharing a corner with the defect-bearing octahedron are generally more tilted (Fig. 4a), but less strained (Fig. 4b), than polyhedra sharing an edge (two oxygen atoms) with the defect site. Expansion (positive strain) of one of the edge-sharing octahedra is explained by a combination of expansion of the shared edge (following expansion of the defect site) and tilt of the surrounding polyhedra. Tilting and strain increase with the size mismatch between the exchanged cations. The defect–O–cation angle appears to vary linearly with the difference in cation radii (Fig. 4a). Variations in strain of the polyhedra in contact with the defect site (Fig. 4b) are more complicated and more site-dependent than variations in defect–O–cation angle. Fig. 4. View largeDownload slide Strain of the forsterite lattice around point defects versus ionic radius (after Shannon, 1976), for tetrahedra (triangles) and octahedra (circles) adjacent to the defect-bearing site. The original forsterite lattice is plotted with the radius of Mg2+. Grey and open symbols are used for polyhedra sharing one bridging oxygen atom (a corner) and two bridging oxygen atoms (an edge) with the defect site, respectively. (a) Tilt of the polyhedra around the defect (see displacement of atoms in Fig. 3). (b) Volumetric strain of the same polyhedra as (a), expressed in comparison with the forsterite lattice. Fig. 4. View largeDownload slide Strain of the forsterite lattice around point defects versus ionic radius (after Shannon, 1976), for tetrahedra (triangles) and octahedra (circles) adjacent to the defect-bearing site. The original forsterite lattice is plotted with the radius of Mg2+. Grey and open symbols are used for polyhedra sharing one bridging oxygen atom (a corner) and two bridging oxygen atoms (an edge) with the defect site, respectively. (a) Tilt of the polyhedra around the defect (see displacement of atoms in Fig. 3). (b) Volumetric strain of the same polyhedra as (a), expressed in comparison with the forsterite lattice. To illustrate the effect of varying mineral structure, elastic properties, cation size and charge, Figure 5 illustrates the results of calculations carried out for eight simple oxides (SrO, BaO, CoO, NiO, FeO, EuO, MgO, CaO) with rocksalt structure where the incorporation of each divalent cation has been simulated as in Figure 2. Results of similar calculations are presented for three isostructural sesquioxides (α-corundum Al2O3, eskolaite Cr2O3 and Sc2O3) where the cation has been exchanged with Al, Fe, La, Sc, Y, Nd, Eu, Cr, Gd, Ho, Yb, Lu, Pu and Mn (all in trivalent form). Computations for selected silicates are also shown in Figure 5a and b: divalent cations have been exchanged with Mg2+ into the M site of enstatite, the M1 sites of forsterite, Mg-celadonite and phlogopite, and the M1B site of chloritoid; trivalent cations have been exchanged into the three M sites of chloritoid hosting Al3+. Results are also shown for divalent cation exchange in the X site of almandine (hosting Fe2+). Fig. 5. View largeDownload slide Links between strain, cation size mismatch and relaxation energy. (a) Variation of cation–oxide bond length within the defect-hosting polyhedron, versus the difference in cation size between host and defect from Shannon (1976), for selected oxides and silicates. Mean bond lengths are used when defect-bearing sites are deformed. Divergence from the 1:1 line shows that lattice strain is not well represented with a hard-sphere model. Correlations appear for oxides with the same structure. (b) Corresponding relaxation energy versus the variation of cation–oxide bond length within the defect-hosting polyhedron, highlighting asymmetric behaviour and differences between oxides, silicates and host site within silicates. (c) Relaxation energy versus bond length variation for Ca2+ incorporated into Mg sites only, for selected minerals, showing a positive correlation and that stiff minerals are less strained than more deformable ones. (d) Relaxed defect energy ( Urelaxedcation→host) as a function of mean bond length for divalent cations in oxides and selected silicates. See text for discussion. Fig. 5. View largeDownload slide Links between strain, cation size mismatch and relaxation energy. (a) Variation of cation–oxide bond length within the defect-hosting polyhedron, versus the difference in cation size between host and defect from Shannon (1976), for selected oxides and silicates. Mean bond lengths are used when defect-bearing sites are deformed. Divergence from the 1:1 line shows that lattice strain is not well represented with a hard-sphere model. Correlations appear for oxides with the same structure. (b) Corresponding relaxation energy versus the variation of cation–oxide bond length within the defect-hosting polyhedron, highlighting asymmetric behaviour and differences between oxides, silicates and host site within silicates. (c) Relaxation energy versus bond length variation for Ca2+ incorporated into Mg sites only, for selected minerals, showing a positive correlation and that stiff minerals are less strained than more deformable ones. (d) Relaxed defect energy ( Urelaxedcation→host) as a function of mean bond length for divalent cations in oxides and selected silicates. See text for discussion. The change of the cation-oxygen bond length (Δbond length) in the site of defect incorporation is plotted in Figure 5aversus the difference in cation size between host and defect, using the values given by Shannon (1976). When sites are deformed (different cation-oxygen bond lengths in the same site), average bond lengths are reported. Positive values for Δbond length mean that the bond length has increased due to defect incorporation, and conversely for negative values. The generally positive correlation shows that bond lengths vary consistently with the size difference of cations. Several trends appear, with simple oxides, sesquioxides, silicates and sites within chloritoid forming individual trends approaching linearity. The divergence from the 1:1 line is largest for simple oxides. This difference implies that bond length in strained crystals is not the sum of fixed atomic radii for defects and oxygen, in other words that atomic radii vary with the crystal structure. This is consistent with ab initio studies on silicates (see Gibbs et al., 2013, for a review) where the electronic structure of oxygen is shown to adapt to the bond length and type. It follows that estimating crystal strain from the difference in atomic radius is inaccurate and should be taken only as a rule of thumb. Figure 5b presents the relaxation energy estimated for each substitution shown in Figure 5a. Strikingly, simple oxides and sesquioxides show a third-degree polynomial dependence of relaxation energy on the change in bond length, as well as silicates and sites within silicates. Relaxation energy is asymmetric with respect to the change in bond length, where the incorporation of small cations implies less strain energy than larger cations (consistent with the analysis of Eshelby, 1951). Relaxation energies within similar structures give at first order comparable dependencies on the change in bond length, regardless of the elastic properties of the host mineral. The asymmetry of Figure 5b is explained by comparing a Ca2+ defect in MgO and a Mg2+ defect in CaO. In MgO, the initial bond length is 2·10 Å and it stretches to 2·20 Å (Δ = 0.10 Å) for a relaxation energy around 90 kJ/mol after incorporating Ca2+. In CaO, the initial bond length is 2·40 Å and it is reduced to 2·27 Å (Δ = -0·13 Å) for a relaxation energy around 60 kJ/mol after incorporating Mg2+. The difference between incorporating larger and smaller cations is noteworthy in terms of energy. For crystal strain this difference explains the (second-order) divergence from a straight line of the trend observed for simple oxides in Figure 5a. Calculations for octahedral sites of silicates show that their response varies between silicates and between sites in the same silicate, as exemplified by chloritoid where a fourfold increase in relaxation energy is observed between divalent cation incorporation in the Mg2+ M1B site and trivalent cation incorporation in the Al3+ M2A or M2B site. The Al3+ M1A site appears intermediate. These relaxation energies are uncorrelated with the elastic properties of individual sites: bulk moduli for octahedral sites of chloritoid have been estimated by Comodi et al. (1992) around 53 GPa for M1A, 63 GPa for M1B, 144 GPa for M2A and 116 GPa for M2B. However, some general trends are observed: first, OH-bearing sites are more strained (greater Δbond length on Fig. 5b) and show higher relaxation energies (see M1B of chloritoid, celadonite, phlogopite). Second, defects incorporated in sites in the vicinity of vacancies are also more strained, with correspondingly higher relaxation energies (compare for example dioctahedral Mg-celadonite to trioctahedral phlogopite). Figure 5c compares relaxation energies to strain after incorporation of a Ca2+ defect into Mg-sites, for the Mg-bearing minerals of Figure 5a and b to which pyrope, clinochlore (M1 and M3), lizardite, actinolite and talc have been added. Periclase MgO and pyrope are less strained and show lower relaxation energies than anhydrous silicates like enstatite and forsterite; phyllosilicates and chloritoid are more deformed and consequently have high relaxation energies. As for Mg-celadonite, defect-bearing sites are very stretched in lizardite and in the M3 site of clinochlore. In lizardite, the Mg-bearing octahedral sites are located in contact with a large interlayer space, in clinochlore the M3 site is part of the interlayer ‘brucite’ sheet. In both these sites associated with interlayer positions, strain is greater than in their mica-like counterpart where octahedral layers are surrounded by stiff tetrahedral layers with which each octahedron shares 4 oxygen atoms (phlogopite, talc and the M1 site of clinochlore in Fig. 4c). The response to strain of the surroundings of the defect-bearing site is an important control on relaxation, which does not solely depend on the elastic properties of a given site but also on the ability of the surrounding polyhedra to deform via tilting. Corresponding defect energies, after relaxation, are presented for simple oxides and selected silicates in Figure 5d. Again, defect energies correlate positively with strain (as measured by the change in bond length or by the difference in cation radius, as these two quantities scale almost linearly in oxides as shown in Fig. 5a). Comparison of defect energies for MgO and NiO, which have cations of similar radii, shows little difference in their response although they have very different elastic properties (for Young’s modulus: MgO = ∼310 GPa, NiO = ∼130 GPa; Soga & Anderson, 1966; Bass, 1995; de Jong et al., 2015). These observations are summarized below: the use of the cation radii derived by Shannon (1976) with fixed-size oxygen is questionable in itself and particularly as a measure of crystal strain for point defects. The excellent agreement of experimentally-derived radii for Mg2+, Co2+, Ni2+, Mn2+ in their constituent oxide of Sasaki et al. (1979) with ab initio modelling (e.g. Freyria-Fava et al., 1993; Orlando et al., 1994; Israel et al., 2003) highlights the interest of accounting for variations in the size of the oxygen atom with bond length; relaxation energies are poorly correlated with the bulk modulus of the defect site itself, mostly because the immediate surroundings of the defect sites are important for strain propagation; relaxation energies are higher in ‘flexible’ environments (where strain is larger) than in densely-packed minerals deforming like stiff elastic wire mesh. This has the counter-intuitive consequence that relaxation energies are higher in hydrated minerals with vacancies (such as the octahedral sheet of phengite) or interlayer spaces (serpentine, chlorite) than in stiff anhydrous minerals like garnet; although linking continuum mechanics and atomistic approaches on point defects has been proposed earlier (e.g. Sahoo, 1984; Garikipati et al., 2006), continuum mechanics is inappropriate to apply to study trace element incorporation in oxides and silicates because of the site dependency of the relaxation energy. Using elastic properties derived for the whole cell, even if elastic anisotropy is considered, is incorrect because the incorporation mechanism matters as well as the nature of the immediate surroundings of the defect site. Results: Heterovalent substitutions A difficulty in modelling and quantifying partitioning with heterovalent substitutions is that the number of possible substitutions may be very high. Heterovalent substitutions necessitate swapping of several ions or creating vacancies to maintain charge balance, which may imply considering distinct crystallographic sites and a large range of cations. Some constraints may be identified to reduce the number of substitutions to investigate. First, defect-bearing polyhedra are thought to be distributed close to one another in order to minimize the charge imbalance distribution (see e.g. Allan et al., 2003). An example of joined Li+ and La3+ substitution in periclase MgO is given below, to compare with Ca2+. Second, as noted previously for equation (9), the efficiency of substitutions involving two trace elements may be expected to be small for statistical reasons. Third, the smallest number of cations involved in the substitution mechanism is likely to minimize crystal strain and, therefore, be energetically favourable. Atomistic simulations for Ba2+ in muscovite are used below to compare two likely substitution mechanisms. Li++La3+=2Mg2+in periclase MgO Figure 6 shows the results of exchanging two Mg2+ ions with two defect cations with different charges and contrasted cation radii (Li+: rM = 0·76 Å; La3+: rM = 1·03 Å; for Mg2+: r0 = 0·72 Å; Shannon, 1976). The defects have been positioned in adjacent Mg octahedra along the <001> plane. The computed displacement (Fig. 6a) is strongly anisotropic, both in terms of direction and amplitude of displacement. It is greatest in close vicinity of the defects. Bond lengths in the Li+-bearing site are compressed by a few percent (Fig. 6b). Bonds in the La3+-bearing site are stretched up to 17% to accommodate the cation 40% larger than Mg2+. Away from the defects, displacement is greatest for atoms aligned with the defects, and decreases quickly (Fig. 6c) but less so than for Ca2+ (compare with Fig. 2). Consequently the displacement is not radial and there is a complex pattern of bond strain with many polyhedra being twisted (Fig. 6b), consistently with observations on homovalent substitutions. Fig. 6. View largeDownload slide Joined incorporation of Li+ (small blue sphere) and La3+ (larger purple sphere) within periclase MgO. (a) Displacement of atoms along the <001> plane as in Figure 2: fixed-size arrows indicate the direction of displacement, interpolated colours show the length of displacement from yellow (large) to blue (negligible). (b) Bond strain for bonds along the <001> plane. Stretching (red tones) is evident around La3+, compression (shades of blue) is clearest around Li+. Note that strain intensity is coloured on a log scale. (c) Displacement of atoms on a log scale versus distance to the Li+ defect (centred). Some points overlap due to symmetry. Fig. 6. View largeDownload slide Joined incorporation of Li+ (small blue sphere) and La3+ (larger purple sphere) within periclase MgO. (a) Displacement of atoms along the <001> plane as in Figure 2: fixed-size arrows indicate the direction of displacement, interpolated colours show the length of displacement from yellow (large) to blue (negligible). (b) Bond strain for bonds along the <001> plane. Stretching (red tones) is evident around La3+, compression (shades of blue) is clearest around Li+. Note that strain intensity is coloured on a log scale. (c) Displacement of atoms on a log scale versus distance to the Li+ defect (centred). Some points overlap due to symmetry. The apparent complexity of displacement and strain highlights the fact that elastic deformation in minerals around atom-scale point defects differs from the radial point of view of continuum mechanics due to the structuring nature of cation-oxygen bonds: periclase is seen as deforming like a grid, not a continuous material, and far from isotropically. Ba2+ in muscovite: vacancy creation vs coupled substitution Muscovite has no divalent cation in its end-member form. It is composed of two tetrahedral layers enclosing a dioctahedral layer, forming sheets separated by interlayer potassium. From the previous considerations, Ba2+ incorporation likely involves the interlayer site, where Ba2+ may substitute either: (1) with K+, and charge balance is ensured through substitution of a neighbouring Si4+ with Al3+ [eq. (10)], or (2) with two neighbouring K+, keeping charge balance via a vacancy [eq. (11)]. Using a previously optimised muscovite structure respecting Al avoidance in tetrahedral layers (Palin & Dove, 2004), Ba2+ and Al3+ have been exchanged with neighbouring K+ and Si4+, with a resulting relaxation energy for equation (10) of 579 kJ/mol. This operation has been repeated exchanging two neighbouring K+, filling one position with Ba2+ and leaving the other empty [eq. (11)]. The corresponding relaxation energy of 196 kJ/mol is almost three times lower. Therefore, equation (11) is predicted to be very significantly favourable to partitioning over equation (10), as partitioning depends exponentially on strain energy [eq. (19)–(23)]. Although equation 11 is predicted to be energetically favourable in the dilute limit, Ba concentrations as high as 20 wt% BaO have been observed and unambiguously attributed to the substitution described in equation (10) (Grapes, 1993; Hetherington et al., 2003). Concentrations in BaO of 18 wt% may theoretically be reached via equation (11), at which point all of the interlayer space will be filled by Ba2+ and vacancies. The observed 20 wt% BaO via equation (10) implies an increase in tetrahedral Al content from one AlIV atom per formula unit in muscovite to 1.5 AlIV atom. Tetrahedral ordering being commonplace in micas where it leads to Al-avoidance below at least 1500 K (e.g. Palin et al., 2001), its ordering pattern must be modified in Ba-rich muscovite with the result of greater stability compared to a vacancy-rich interlayer space. This is consistent with previous studies showing that interlayer vacancy-rich micas are unstable due to long-range ordering (Dubacq et al., 2011, and references therein). From this it is concluded that in muscovite, the efficiency of interlayer vacancy creation is greater than that of the mechanism involving exchange in tetrahedral positions in the dilute limit. At higher Ba concentration, Ba-rich micas are explained by equation (10), and not equation (11) as long-range ordering becomes more potent. COMPARING ATOMISTIC MODELLING TO EXPERIMENTS AND PREDICTING PARTITION COEFFICIENTS Mineral/liquid partitioning Garnet/melt partitioning, and the influence of mineral composition The effect of mineral composition on mineral/liquid partitioning is illustrated here with Sr exchange between garnet and anhydrous melt, using the experimental data of van Westrenen et al. (1999). Figure 7 shows the partition coefficient for Sr calculated from their ∼1550 °C experiments as a function of the pyrope content of garnet along the pyrope–grossular solid solution ((Ca-Mg)3Al2Si3O12). A striking observation is that the partition coefficient is maximal around 65% pyrope, and decreases very steeply with increasing pyrope content. Assuming the measurements are representative of equilibrium, these steep changes must indicate a change in the thermodynamics of the melt species: the pyrope–grossular solid solution cannot explain these abrupt changes as it has been shown to vary smoothly and with near ideality at these temperatures (see for example activity-compositions relations in Vinograd & Sluiter, 2006). Fig. 7. View largeDownload slide Influence of composition on garnet/melt partitioning for Sr along the pyrope–grossular solid solution at 1550 °C. (a) Values of exp(-UεSr→Grt/R.T) [in the form of eq. (19)] calculated from the experimental data of van Westrenen et al. (1999) (open circles), and with relaxation energies from atomistic modelling (filled circles) for pyrope and grossular. Grey circles are calculated from the end-members with an ideal mixing model (dotted line). Note the y-axis is in a log scale. (b) Corresponding garnet/melt partition coefficient for Sr. Experimental uncertainties are smaller than symbols when not shown. Fig. 7. View largeDownload slide Influence of composition on garnet/melt partitioning for Sr along the pyrope–grossular solid solution at 1550 °C. (a) Values of exp(-UεSr→Grt/R.T) [in the form of eq. (19)] calculated from the experimental data of van Westrenen et al. (1999) (open circles), and with relaxation energies from atomistic modelling (filled circles) for pyrope and grossular. Grey circles are calculated from the end-members with an ideal mixing model (dotted line). Note the y-axis is in a log scale. (b) Corresponding garnet/melt partition coefficient for Sr. Experimental uncertainties are smaller than symbols when not shown. The melt composition is needed in the partitioning equation [eq. (19)] to calculate the mineral/melt ratio for the host cation, therefore, a continuous partitioning equation from pyrope to grossular cannot be calculated without knowledge of the associated melt composition. Comparison of experiments with atomistic modelling of relaxation energies is given in Figure 7a through the calculation of exp(-UεSr→grt/R.T) [see eq. (19)], which is possible both from experimental results and modelling. Figure 7a shows the value of exp(-UεSr→grt/R.T) calculated using the data of van Westrenen et al. (1999) and with atomistic modelling for the pyrope and grossular end-members, linked by an ideal mixing model along the pyrope–grossular solid solution. The computed value for grossular exp-UεSr→GroR.T=0·17 is much higher than that for pyrope exp-UεSr→PyR.T=8·9*10-8. This results in an extremely steep decrease of the exponential term—and, therefore, of the expected partition coefficient—over five orders of magnitude between about 95 and 100% pyrope in garnet. At lower pyrope content, the calculated partitioning is controlled by the value for grossular, which is in the range of the low-pyrope experiment. Calculated values fall within 1.3 order of magnitude of the measurements (well above estimated analytical uncertainties, with a reduced χ2 ≈ 2000 calculated following Bevington & Robinson, 2002). It has been suggested through review that uncertainties in the study of van Westrenen et al. (1999) were possibly underestimated. The calculated value for χ2 would be lower (better fit) using larger uncertainties, and we estimate that the current discussion remains unaffected unless uncertainties are underestimated by more than two orders of magnitude, which is unlikely.Corresponding partition coefficients are shown on Figure 7b as grey circles and exhibit the same misfit. As expected, ideal mixing in garnet does not reproduce the abrupt changes in exp(-UεSr→grt/R.T) for experiments with intermediate compositions. In this view, the increase of the partition coefficient to a maximum of around 70% pyrope (Fig. 7b) is best explained by variations of the chemical potential of melt species. This example highlights the difficulty of comparing work on synthetic systems, natural systems and atomistic calculations for end-members: not only the mineral/liquid ratio of the host cation must be known, but also a small fraction of a cation with relatively large ionic radius in the exchange site of the host crystal predominantly with small ionic radius (such as a few percent Ca in a Mg-site) may affect the partition coefficient for large cations over orders of magnitude. This is a reminder of the ‘plateau effect’, where the effect of chemical mixing at small concentrations has been shown to be non-linear and linked to the size of strain fields around defects (e.g. Salje, 1995). The plateau effect is interpreted as a saturation effect (Carpenter et al., 2009) and the non-linearity as the transition from a classical chemical mixing regime to a regime dominated by the defect content and the overlap between strain clouds around defects. In Figure 7b, pyrope fractions below 95–98 % correspond to the first regime when at high pyrope content, the calculated partition coefficients are extremely sensitive to the grossular content. Away from the plateau effect, the presence in the host crystal of a cation with which exchange is more difficult alters the partition coefficient less dramatically (such as several percent Mg in the Ca-site of grossular, Fig. 7a). This is consistent with the asymmetric shape of strain energy with respect to defect size (shown in Fig. 5b and in the analysis of Eshelby, 1951, 1956). Olivine/melt partitioning, and deviation from ideality Figure 8 shows the application of equation (19) to the partitioning of divalent cations between olivine and melt. Strain energy has been calculated for defects in forsterite and fayalite, as shown in Figure 8a at 1150 °C versus cation radius. Defects were incorporated in the M1 site for Mn2+ and smaller ions, in the M2 site for Ca2+ and larger ions. In both cases, the exponential term controlling partitioning in equation (19) shows a maximum at unity for the host cation ( UεMg→Fo=0) and again with strongly asymmetric shape. For cations larger than Ca2+, high relaxation energies will result in extremely small olivine / melt partition coefficients, regardless of the melt composition (on the order of 10–10 and smaller). For this reason, Figure 8b focuses on olivine/melt partition coefficients DR2+ol/melt estimated for cations with radii between that of Ni and Ca, calculated using experimental melt composition and an ideal activity model for olivine along the forsterite–fayalite join. Experimental results have been selected as representative from Mysen (1976), Duke (1976) and Häkli & Wright (1967) for Ni, and from the compilation of Libourel (1999) for Ca. Experimental values are reasonably well reproduced for Ni, but are under-estimated by a factor of 30 for Ca in forsterite (highlighted on Fig. 8b). This discrepancy is too large to originate solely from variations in the activity of the host cation ( aMgFo/ aMgL). The mismatch between the Gibbs energy of reaction for equation (19) for Ca in forsterite and strain energy may rather come from: Fig. 8. View largeDownload slide Strain energy and mineral/melt partition coefficients for homovalent substitution in olivine. (a) Calculated values of exp(-UεR2+→Ol/RT) in forsterite and fayalite, at T = 1150 °C. Note y-axis is in a log scale. (b) Olivine/melt partition coefficient ln(DR2+Ol/melt) for Ca and smaller elements [red box in (a)], calculated with equation (19) for experimental mineral/melt compositions. Squares show the corresponding measured partition coefficients. Note the mismatch for Ca in forsterite, see text for discussion. Fig. 8. View largeDownload slide Strain energy and mineral/melt partition coefficients for homovalent substitution in olivine. (a) Calculated values of exp(-UεR2+→Ol/RT) in forsterite and fayalite, at T = 1150 °C. Note y-axis is in a log scale. (b) Olivine/melt partition coefficient ln(DR2+Ol/melt) for Ca and smaller elements [red box in (a)], calculated with equation (19) for experimental mineral/melt compositions. Squares show the corresponding measured partition coefficients. Note the mismatch for Ca in forsterite, see text for discussion. incorporation of Ca via a different substitution mechanism. The Mg–Ca exchange equilibrium has been thoroughly studied by many authors and this is discarded as unlikely; poor estimation of the strain energy. This is discarded as the sole explanation because strain energy would have to be overestimated by a factor of 1.8 to explain the mismatch, probably much larger than modelling uncertainty; strong non-ideality of Ca in forsterite. The forsterite – monticellite solid solution is notoriously non-ideal, with Margules parameters (controlling non-ideality in solid solutions) on the order of 50 kJ/mol (e.g. Adams & Bishop, 1985) but not large enough to explain the difference in partitioning for low Ca concentrations. Typical Ca content of olivine in experiments are in the range 0·5–1·5 wt% CaO, corresponding to an average distance between Ca defects in the range ∼15 to ∼10 Å if distributed homogenously into the olivine lattice. These distances are in the high range for strain fields around Ca in forsterite (see Fig. 3), but it is arguable that the discrepancy is partly due to the plateau effect. This is consistent with lower values for DCa2+ol/melt in Ca-poor melts, however Libourel (1999) has shown that this decrease originates mostly from the melt structure, at least for melts in the range 2–30 wt% CaO; strong non-ideality of Ca in the melt. As Ca is a major melt component, aCamelt directly affects and is affected by melt structure (e.g. Watson, 1979). Libourel (1999) unambiguously show that DCa2+ol/melt is sensitive to melt speciation, especially to variations in Si and Na acting respectively as network-forming and network-modifying cations. A combination of the last three sources of uncertainties is retained to explain the mismatch, and further investigations would be of interest to quantify deviation from the assumption that strain energy is predominant on partitioning as a function of concentration. In the compilation of Libourel (1999), it is shown that DCa2+ol/melt increases by a factor of ∼10 from forsterite to fayalite. This is consistent with the difference in the term exp(-UεOl/R.T) in forsterite and fayalite (Fig. 8a), which is a factor of 10 at 1500 °C and 20 at 1100 °C. It remains clear that strain energy and its variations are part of the controlling factors in mineral / melt partitioning. Mineral/mineral partitioning For partitioning between minerals, it was shown above that partition and distribution coefficients between solid phases may be calculated independently from the composition of the associated melt (or metamorphic fluid), with the help of differences in the incorporation energy of trace elements into mineral lattices [see eq. (13)–(17)], or Pickles et al., 2016). It is also shown that in some cases, using defect energies rather than relaxation energies to estimate the incorporation energy into the exponential term [see for example eq. (23)] releases the assumption of strain energy being the preponderant term in the exchange reaction. Using relaxation energy or defect energy may alter the calculated partition coefficient over orders of magnitude. This is investigated here, studying the partitioning of some divalent cations and (trivalent) rare earth between garnet and clinopyroxene. Many experimental results are available for REE, generally showing preferential incorporation of the largest REE into clinopyroxene at the expense of garnet (e.g. Shimizu, 1975). We selected the experimental results of Tuff & Gibson (2007) on picritic melts. Their work presents partition coefficients broadly consistent with studies on basaltic compositions and has been selected over more recent studies such as Pickles et al. (2016) due to the simpler composition of the reported mineral assemblages, with minor Cr and negligible Fe3+ in garnet and clinopyroxene. Apart from Sc, REE are incorporated into the X site of garnet and into the M2 site of clinopyroxene (e.g. Blundy & Wood 1994, 2003; van Westrenen et al., 1999; Gaspar et al., 2008). These sites are occupied by eight-coordinated, chiefly divalent cations (Ca2+ in grossular and diopside): trivalent cation incorporation implies a Tschermak-like substitution for charge balance. In the modelling, Sc3+ replaces octahedral Al3+; divalent cations exchange with homovalent X-site cations of garnet; in clinopyroxene the M2 site hosts Sr2+ and Ba2+, and Co2+ is incorporated into the M1 site. The resulting exchange equation between garnet and clinopyroxene for eight-coordinated trivalent defects, here written for the Ca end-members, is: Ca3Al2Si3O12+R3+Mg(Al, Si)O6 ⇔ Ca2R3+Al2(Al,Si2)O12+CaMgSi2O6 (24) Grossular+defect-rich Cpx ⇔ defect-rich garnet+diopside Equation (24) is independent of the liquid composition and a partition coefficient may be written as for the olivine/orthopyroxene exchange [eq. (20)]. The present atomistic approach requires knowledge of the cation distribution into sites and must use discrete positions for the exchanged cations to ensure correct charge distribution. This is a major problem for augite where the M1 site is filled with Mg2+, Fe2+ and Al3+ and the M2 site is distorted following incorporation of cations with dissimilar charge, radius and coordination number (Na+, Ca2+, Mg2+). Ordering in augite consists of: (1) short-range ordering minimizing charge imbalance (for example due to a jadeite component, Na+ in M2 and Al3+ in M1 are adjacent); and (2) long-range ordering (e.g. Al avoidance) which is more temperature-sensitive (see for example Wood, 1974; Bosenick et al., 2001; Flemming & Luth, 2002). Intermediate compositions may, therefore, be modelled either: (1) through the use of end-members and activity models; (2) via large super-cells where ordering is known with the aid of Monte Carlo simulations; or (3) using the mean field approximation where site properties are an average of individual cation properties. Ordering in clinopyroxene being far beyond the scope of this study, the mean field approximation has been selected for simplicity: additional charge brought by trivalent defects has been compensated by equivalent exchange of Al3+ with Si4+ in adjacent tetrahedral sites, in minerals with the experimental compositions reported by Tuff & Gibson (2007). In their P-511 experiment, garnet has the composition Na.01Ca.36Mg1.94 Fe·692+Mn·01Cr.02Al1·92Ti·04Si3O12 and has been modelled as Py65Alm23Grs12 (abbreviations after Whitney & Evans, 2010). Clinopyroxene has an augitic composition Na.09Ca·35Mg1·06 Fe·282+ Al·18VITi·01 Al·09IVSi1·91O6; this has been simplified to Na·10Ca·35Mg1·10 Fe·252+ Al·20VI Al·10IVSi1·90O6. For comparison, calculations are also reported for end-members. Figure 9 presents calculated partition coefficients for the experimental compositions and for the Ca and Mg end-members. Experimental data from Pickles et al. (2016) have been added for comparison and information on divalent cations. As expected from garnet/melt partitioning (Fig. 7) and discussed above, changes in mineral composition dramatically affect partitioning: 4 to 6 orders of magnitude of difference between Ca and Mg end-members are shown for most elements in Figure 9a. For minerals with experimental compositions, results are intermediate (mixing is not necessarily linear, as in Fig. 7a). The partition coefficients estimated using defect energy (Fig. 9a) also differ from the ones obtained using strain energy (Fig. 9b) by orders of magnitude. The overall shape defined by the experimental results is well reproduced in Figure 9a (apart for Ba). The strain energy approach (Fig. 9b) shows much poorer agreement. Fig. 9. View largeDownload slide Comparison of garnet/pyroxene partition coefficients from experimental work and atomistic calculations. Values in (a) are calculated using defect energy, in (b) with strain energy. Co, Sr and Ba are calculated for homovalent substitution in divalent state, all other calculations are for trivalent cations. Sc3+ incorporation is modelled as homovalent in the Al-site (octahedral). Larger REE are modelled as incorporating eight-coordinated sites via coupled substitutions, as explained in the text. Green circles labelled Grt/Cpx show partitioning calculated for the experimental composition of Tuff & Gibson (2007) given in legend. Selected data from Pickles et al. (2016) allow comparison for divalent elements. Results are also shown for the same substitutions into the Mg and Ca end-members of garnet and pyroxene (Fe not shown for readability). Ca end-members: grossular and diopside (for cations with CN = 8) plus Ca–Tschermack pyroxene (kushiroite, for cations with CN = 6). Mg end-members are pyrope and enstatite. CN, Coordination number. Fig. 9. View largeDownload slide Comparison of garnet/pyroxene partition coefficients from experimental work and atomistic calculations. Values in (a) are calculated using defect energy, in (b) with strain energy. Co, Sr and Ba are calculated for homovalent substitution in divalent state, all other calculations are for trivalent cations. Sc3+ incorporation is modelled as homovalent in the Al-site (octahedral). Larger REE are modelled as incorporating eight-coordinated sites via coupled substitutions, as explained in the text. Green circles labelled Grt/Cpx show partitioning calculated for the experimental composition of Tuff & Gibson (2007) given in legend. Selected data from Pickles et al. (2016) allow comparison for divalent elements. Results are also shown for the same substitutions into the Mg and Ca end-members of garnet and pyroxene (Fe not shown for readability). Ca end-members: grossular and diopside (for cations with CN = 8) plus Ca–Tschermack pyroxene (kushiroite, for cations with CN = 6). Mg end-members are pyrope and enstatite. CN, Coordination number. For REE, Co and Sr, results are in agreement with experimental values within an order of magnitude when using defect energy. However, discrepancies are clear both in terms of absolute values and variations ( DR3+grt/cpx too low for La3+, La3+/Yb3+ ratio too high). These differences likely arise from the simplified augite model used here, and from the precision of the present atomistic modelling. Partitioning values estimated for Ba2+ are five orders of magnitude lower than values reported by Tuff & Gibson (2007). Assuming that the Ba content was well measured in garnet, this very large difference is not straightforward to explain. It is noteworthy that Ba2+ is much larger than other defects investigated here, including lanthanides (La3+ ∼ 1·16 Å; Ba2+ ∼1·42 Å when eight-coordinated according to Shannon, 1976). The discrepancy may first indicate that the incorporation mechanism for Ba is not homovalent substitution in the Ca site [as in eq. (24)], at least in one of the two minerals. The Ca site being the largest site in garnet and clinopyroxene, it is a likely host for Ba2+ incorporation, but under-estimated values for DBa2+grt/cpx suggest that the defect energy should be lower for Ba in garnet. This may be the case by coupled substitution of several defects, in the manner of incorporating a Li6BaLa2Ta2O12 end-member as synthetized for lithium-ion batteries (e.g. Thangadurai & Weppner, 2005). Second, it is possible that Ba incorporation causes migration of atoms around the defect site that would be not reproduced here. Further discussion and perspectives Using defect energy rather than strain energy when calculating partitioning between minerals better reproduces experimental results, both in terms of absolute values and variations. This is explained by the greater adequacy of defect energy to represent bond enthalpy over strain energy. In other words, cation exchange is not only controlled by cation charge and radius: for example octahedral Al3+ and Fe3+ have comparable octahedral size (0·535 Å and 0·55 Å, respectively, i.e. within 3%) and Vieillard (1994a, b) estimated that their contribution to the formation enthalpy of silicates generally varies within 10% (using an empirically-derived parameter Δ HOcation2-(comp)). More generally, the agreement of calculations with experimental data shown here is much less good than what is achieved by fitting the lattice strain model to experimental values, for example in the study of Pickles et al. (2016). It must be emphasized that the good agreement of the lattice strain model to garnet/clinopyroxene partitioning experiments is only obtained after fitting, by altering the values of the site’s apparent Young modulus and host cation radius. Agreement between the present calculations and experiments may be significantly improved by the use of Monte Carlo simulations for better description of ordering. Moreover, it is noteworthy that defect energies computed for garnet and clinopyroxene are close in absolute values (within 10% for REE). As the partition coefficient depends exponentially on the difference between these values, accuracy within 1% would be welcome. Unfortunately such accuracy is not guaranteed with atomistic simulations based on force fields methods as carried out here, and partitioning simulations would benefit from more sophisticated modelling. Specific points of interest for more rigorous thermodynamic analysis include systematic experimental and modelling studies of strain fields around defects and their variations within solid solutions, to evaluate concentrations at which defects truly behave like trace elements, without interacting with each other, and to estimate effects due to the addition of defects for comparison with natural systems. It remains that the present approach has the advantage of being relatively simple, computationally fast and not based on fitting. SUMMARY AND CONCLUSIONS In an attempt to explain and quantify trace element partitioning in metamorphic assemblages, the controlling processes of cation incorporation into oxides and silicates have been reviewed. Cation exchange in crystal sites successfully explains crystal/melt partitioning, and may be modelled with exchange reactions involving end-members and thermodynamic calculations. This approach has met considerable success through the lattice strain model of Blundy & Wood (1994). However the approach of the lattice strain model is limited. The use of the Brice (1975) equation is inadequate to solve a radial strain problem because of an inconsistency in the definition of strain, with the result that both the shear and bulk modulus should appear in the equation for correct modelling of strain, even in isotropic material. Applying continuum mechanics to study deformation at atomic scale in silicates is also fundamentally an oversimplification, due to the non-continuous nature of crystal strain around point defects. With the help of classical atomistic modelling, the limitations of hard-sphere models (with fixed radii for cations and oxygen) have been highlighted to show that on average the cation–oxygen bond length in crystals is less altered by cation exchange in trace amounts than the difference in cation radii as calculated with the values of Shannon (1976). Crystal strain depends not only on the difference in cation radii, but also on the environment of the exchange site, where more anisotropic and, or, less densely-packed structures allow more deformation, by anisotropic strain and tilting of the surrounding polyhedra. Phyllosilicates are particularly prone to large strain values as they frequently contain vacancies and ‘flexible’ interlayer sites. These limits explain part of the differences between measured partition coefficients and values estimated from the lattice strain model. They add up to the assumption that crystal strain is the preponderant control on mineral/fluid partitioning: neglecting the chemical potential of fluid species prevents modelling of several mechanisms altering partition coefficients, such as the effect of pressure and composition on the structure and polymerization of melts and speciation for metamorphic fluids, effects that are explained by reaction, mixing and possible non-ideality between fluid phase components. Atomistic modelling gives fast access to relaxation energies (equivalent to strain energies) and defect energies. Using relaxation energies solely to estimate partition coefficients between mineral and melt again relies on the assumption that crystal strain energy is preponderant in the mineral–melt exchange reaction. This assumption may be lifted to estimate partition coefficients between two minerals by using defect energies. However, defect energy may only be used to estimate mineral/fluid partition coefficients if the chemical potentials of the fluid species are known. The defect and relaxation energies presented here give useful constraints on mineral/melt and mineral/mineral partitioning, but more calculations and more sophisticated computation techniques (such as ab initio calculations) are required in conjunction with better knowledge of the structure of the liquid. First, the precision of calculations is important as defect energies between similar sites in different crystals may be within a few percent, yet partitioning equations depend exponentially on their difference. Second, it is shown that when a small fraction of a cation with which exchange is energetically favourable are present in a crystal site less apt to deform, such as a few percent grossular into pyrope, the calculated partition coefficient is altered over orders of magnitude. In this example, the energetics of exchange between eight-coordinated cations and a divalent defect such as Sr2+ are very sensitive to Ca content close to the pyrope end-member. The importance of considering the length scale of strain fields around defects is emphasized, in addition to evaluating short and long-range ordering in minerals for precise estimation of variations of partition coefficients along solid solutions. This would benefit petrological studies of fluid–rock interactions at crustal metamorphic temperatures, where phyllosilicates are ubiquitous and ordering in crystals is more important than at magmatic conditions. ACKNOWLEDGEMENTS B.D. is indebted to B. Villemant, A. Verlaguet, O. Shorttle and C. Chopin for helpful, lengthy discussions and for commenting on the paper. Detailed and constructive reviews by T. Hammouda, V. van Hinsberg, an anonymous reviewer and editor J. Hermann benefited the manuscript and are gratefully acknowledged. FUNDING This work was supported by INSU (CNRS). A.P. acknowledges funding through the ERC project SINK (306810). SUPPLEMENTARY DATA Supplementary data for this paper are available at Journal of Petrology online. REFERENCES Adam J. , Green T. ( 2006 ). Trace element partitioning between mica- and amphibole-bearing garnet lherzolite and hydrous basanitic melt: 1. Experimental results and the investigation of controls on partitioning behaviour . Contributions to Mineralogy and Petrology 152 , 1 – 17 . Google Scholar CrossRef Search ADS Adams G. , Bishop F. ( 1985 ). An experimental investigation of thermodynamic mixing properties and unit-cell parameters of forsterite-monticellite solid solution . American Mineralogist 70 , 714 – 722 . Allan N. L. , Du Z. , Lavrentiev M. Y. , Blundy J. D. , Purton J. A. , van Westrenen W. ( 2003 ). Atomistic simulation of mineral-melt trace-element partitioning . Physics of the Earth and Planetary Interiors 139 , 93 – 111 . Google Scholar CrossRef Search ADS Allègre C. J. , Treuil M. , Minster J. F. , Minster B. , Albarède F. ( 1977 ). Systematic use of trace element in igneous process . Contributions to Mineralogy and Petrology 60 , 57 – 75 . Google Scholar CrossRef Search ADS Andrieux P. , Petit S. ( 2010 ). Hydrothermal synthesis of dioctahedral smectites: the Al-Fe3+ chemical series: Part I: Influence of experimental conditions . Applied Clay Science 48 , 5 – 17 . Google Scholar CrossRef Search ADS Annersten H. , Adetunji J. , Filippidis A. ( 1984 ). Cation ordering in Fe-Mn silicate olivines . American Mineralogist 69 , 1110 – 1115 . Armbruster T. , Geiger C. A. , Lager G. ( 1992 ). Single-crystal X-ray structure study of synthetic pyrope almandine garnets at 100 and 293 K . American Mineralogist 77 , 512 – 521 . Banno S. , Matsui Y. ( 1973 ). On the formulation of partition coefficients for trace elements distribution between minerals and magma . Chemical Geology 11 , 1 – 15 . Google Scholar CrossRef Search ADS Bass J. D. ( 1995 ). Elasticity of Minerals, Glasses, and Melts . American Geophysical Union 45 – 63 . Baxter E. F. , Caddick M. J. ( 2013 ). Garnet growth as a proxy for progressive subduction zone dehydration . Geology 41 , 643 – 646 . Google Scholar CrossRef Search ADS Bebout G. E. , Bebout A. E. , Graham C. M. ( 2007 ). Cycling of B, Li, and LILE (K, Cs, Rb, Ba, Sr) into subduction zones: SIMS evidence from micas in high-P/T metasedimentary rocks . Chemical Geology 239 , 284 – 304 . Google Scholar CrossRef Search ADS Bebout G. E. , Ryan J. G. , Leeman W. P. , Bebout A. E. ( 1999 ). Fractionation of trace elements by subduction-zone metamorphism—effect of convergent margin thermal evolution . Earth and Planetary Science Letters 171 , 63 – 81 . Google Scholar CrossRef Search ADS Bebout G. E. , Agard P. , Kobayashi K. , Moriguti T. , Nakamura E. ( 2013 ). Devolatilization history and trace element mobility in deeply subducted sedimentary rocks: Evidence from Western Alps HP/UHP suites . Chemical Geology 342 , 1 – 20 . Google Scholar CrossRef Search ADS Bevington P. , Robinson D. ( 2002 ). Data Reduction and Error Analysis for the Physical Sciences , 3rd edn. New York : McGraw-Hill Higher Education , 352 p. Berthelot M. ( 1872 ). Sur les lois qui président au partage d’un corps entre deux dissolvants (théorie) . Annales de Chimie et de Physique 26 , 408 – 417 . Birle J. D. , Gibbs G. V. , Moore P. B. , Smith J. V. ( 1968 ). Crystal structures of natural olivines . American Mineralogist 53 , 807 – 824 . Blundy J. , Wood B. ( 1994 ). Prediction of crystal–melt partition coefficients from elastic moduli . Nature 372 , 452 – 454 . Google Scholar CrossRef Search ADS Blundy J. , Wood B. ( 2003 ). Partitioning of trace elements between crystals and melts . Earth and Planetary Science Letters 210 , 383 – 397 . Google Scholar CrossRef Search ADS Bosenick A. , Dove M. T. , Myers E. R. , Palin E. J. , Sainz-Diaz C. I. , Guiton B. S. , Warren M. C. , Craig M. S. , Redfern S. A. T. ( 2001 ). Computational methods for the study of energies of cation distributions: applications to cation-ordering phase transitions and solid solutions . Mineralogical Magazine 65 , 193 – 219 . Google Scholar CrossRef Search ADS Bourdelle F. , Parra T. , Beyssac O. , Chopin C. , Vidal O. ( 2013 ). Clay minerals as geo-thermometer: a comparative study based on high spatial resolution analyses of illite and chlorite in Gulf Coast sandstones (Texas, U.S.A.) . American Mineralogist 98 , 914 – 926 . Google Scholar CrossRef Search ADS Brice J. ( 1975 ). Some thermodynamic aspects of the growth of strained crystals . Journal of Crystal Growth 28 , 249 – 253 . Google Scholar CrossRef Search ADS Carpenter M. A. , McKnight R. E. A. , Howard C. J. , Zhou Q. , Kennedy B. J. , Knight K. S. ( 2009 ). Characteristic length scale for strain fields around impurity cations in perovskites . Physical Review B 80 , 214101 . Google Scholar CrossRef Search ADS Cartier C. , Hammouda T. , Doucelance R. , Boyet M. , Devidal J. L. , Moine B. ( 2014 ). Experimental study of trace element partitioning between enstatite and melt in enstatite chondrites at low oxygen fugacities and 5 GPa . Geochimica et Cosmochimica Acta 130 , 167 – 187 . Google Scholar CrossRef Search ADS Cochain B. , Sanloup C. , de Grouchy C. , Crépisson C. , Bureau H. , Leroy C. , Kantor I. , Irifune T. ( 2015 ). Bromine speciation in hydrous silicate melts at high pressure . Chemical Geology 404 , 18 – 26 . Google Scholar CrossRef Search ADS Comodi P. , Mellini M. , Zanazzi P. F. ( 1992 ). Magnesiochloritoid: compressibility and high pressure structure refinement . Physics and Chemistry of Minerals 18 , 483 – 490 . Google Scholar CrossRef Search ADS De Andrade V. , Susini J. , Salomé M. , Beraldin O. , Rigault C. , Heymes T. , Lewin E. , Vidal O. ( 2011 ). Submicrometer Hyperspectral X-ray imaging of heterogeneous rocks and geomaterials: applications at the Fe K-edge . Analytical Chemistry 83 , 4220 – 4227 . Google Scholar CrossRef Search ADS de Jong M. , Chen W. , Angsten T. , Jain A. , Notestine R. , Gamst A. , Sluiter M. , Ande C. K. , van der Zwaag S. , Plata J. J. , Toher C. , Curtarolo S. , Ceder G. , Persson K. , Asta M. ( 2015 ). Charting the complete elastic properties of inorganic crystalline compounds . Scientific Data 2 , 150009 . Google Scholar CrossRef Search ADS Deschamps F. , Godard M. , Guillot S. , Chauvel C. , Andreani M. , Hattori K. , Wunder B. , France L. ( 2012 ). Behavior of fluid-mobile elements in serpentines from abyssal to subduction environments: Examples from Cuba and Dominican Republic . Chemical Geology 312-313 , 93 – 117 . Google Scholar CrossRef Search ADS Diego G. G. , Merlini M. , Valdrè G. , Liermann H. P. , Nénert G. , Rothkirch A. , Kahlenberg V. , Pavese A. ( 2012 ). On the crystal structure and compressional behavior of talc: a mineral of interest in petrology and material science . Physics and Chemistry of Minerals 40 , 145 – 156 . Diener J. F. A. , Powell R. ( 2012 ). Revised activity-composition models for clinopyroxene and amphibole . Journal of Metamorphic Geology 30 , 131 – 142 . Google Scholar CrossRef Search ADS Dragovic B. , Samanta L. M. , Baxter E. F. , Selverstone J. ( 2012 ). Using garnet to constrain the duration and rate of water-releasing metamorphic reactions during subduction: an example from Sifnos, Greece . Chemical Geology 314–317 , 9 – 22 . Google Scholar CrossRef Search ADS Dubacq B. , Bickle M. J. , Evans K. A. ( 2013 ). An activity model for phase equilibria in the H2O-CO2-NaCl system . Geochimica et Cosmochimica Acta 110 , 229 – 252 . Google Scholar CrossRef Search ADS Dubacq B. , Vidal O. , De Andrade V. ( 2010 ). Dehydration of dioctahedral aluminous phyllosilicates: thermodynamic modelling and implications for thermobarometric estimates . Contributions to Mineralogy and Petrology 159 , 159 – 174 . Google Scholar CrossRef Search ADS Dubacq B. , Vidal O. , Lewin E. ( 2011 ). Atomistic investigation of the pyrophyllitic substitution and implications on clay stability . American Mineralogist 96 , 241 – 249 . Google Scholar CrossRef Search ADS Duke J. M. ( 1976 ). Distribution of the period four transition elements among olivine, calcic clinopyroxene and mafic silicate liquid: experimental results . Journal of Petrology 17 , 499 – 521 . Google Scholar CrossRef Search ADS El Korh A. , Schmidt S. T. , Ulianov A. , Potel S. ( 2009 ). Trace element partitioning in HP-LT metamorphic assemblages during Subduction-related metamorphism, Ile de Groix, France: a detailed LA-ICPMS study . Journal of Petrology 50 , 1107 – 1148 . Google Scholar CrossRef Search ADS Eshelby J. ( 1951 ). The force on an elastic singularity . Philosophical Transactions of the Royal Society of London A: Mathematical, Physical and Engineering Sciences 244 , 87 – 112 . Google Scholar CrossRef Search ADS Eshelby J. ( 1956 ). The continuum theory of lattice defects. In: Seitz F. , Turnbull D. (eds.) Solid State Physics, Volume 3 , New York : Academic Press , pp. 79 – 144 . Evans B. ( 2004 ). The serpentinite multisystem revisited: chrysotile is metastable . International Geology Review 46 , 479 – 506 . Google Scholar CrossRef Search ADS Evans B. W. , Yang H. ( 1998 ). Fe-Mg order-disorder in tremolite-actinolite ferro-actinolite at ambient and high temperature . American Mineralogist 83 , 458 – 475 . Google Scholar CrossRef Search ADS Evans K. , Powell R. ( 2006 ). A method for activity calculations in saline and mixed solvent solutions at elevated temperature and pressure: A framework for geological phase equilibria calculations . Geochimica et Cosmochimica Acta 70 , 5488 – 5506 . Google Scholar CrossRef Search ADS Finger L. W. , Hazen R. M. ( 1980 ). Crystal structure and isothermal compression of Fe2O3, Cr2O3, and V2O3 to 50 kbars . Journal of Applied Physics 51 , 5362 – 5367 . Google Scholar CrossRef Search ADS Fiquet G. , Richet P. , Montagnac G. ( 1999 ). High-temperature thermal expansion of lime, periclase, corundum and spinel . Physics and Chemistry of Minerals 27 , 103 – 111 . Google Scholar CrossRef Search ADS Flemming R. L. , Luth R. W. ( 2002 ). 29Si MAS NMR study of diopside—Ca-Tschermak clinopyroxenes: detecting both tetrahedral and octahedral Al substitution . American Mineralogist 87 , 25 – 36 . Google Scholar CrossRef Search ADS Frei D. , Liebscher A. , Franz G. , Wunder B. , Klemme S. , Blundy J. ( 2008 ). Trace element partitioning between orthopyroxene and anhydrous silicate melt on the lherzolite solidus from 1.1 to 3.2 GPa and 1, 230 to 1, 535°C in the model system Na2O-CaO-MgO-Al2O3-SiO2 . Contributions to Mineralogy and Petrology 157 , 473 – 490 . Google Scholar CrossRef Search ADS Freyria-Fava C. , Dovesi F. , Saunders V. R. , Leslie M. , Roetti C. ( 1993 ). Ca and Be substitution in bulk MgO: ’ab initio’ Hartree-Fock and ionic model supercell calculations . Journal of Physics: Condensed Matter 5 , 4793 – 4804 . Google Scholar CrossRef Search ADS Gale J. D. ( 1997 ). GULP: A computer program for the symmetry-adapted simulation of solids . Journal of the Chemical Society, Faraday Transactions 93 , 629 – 637 . Google Scholar CrossRef Search ADS Gale J. D. , Rohl A. L. ( 2003 ). The General Utility Lattice Program (GULP) . Molecular Simulation 29 , 291 – 341 . Google Scholar CrossRef Search ADS Ganguly J. , Saxena S. K. ( 1987 ). Mixtures and Mineral Reactions . In: Minerals, Rocks and Mountains, 19 . Heidelberg : Springer Science & Business Media . 291 p. Google Scholar CrossRef Search ADS Garçon M. , Chauvel C. , Bureau S. ( 2011 ). Beach placer, a proxy for the average Nd and Hf isotopic composition of a continental area . Chemical Geology 287 , 182 – 192 . Google Scholar CrossRef Search ADS Garikipati K. , Falk M. , Bouville M. , Puchala B. , Narayanan H. ( 2006 ). The continuum elastic and atomistic viewpoints on the formation volume and strain energy of a point defect . Journal of the Mechanics and Physics of Solids 54 , 1929 – 1951 . Google Scholar CrossRef Search ADS Gaspar M. , Knaack C. , Meinert L. D. , Moretti R. ( 2008 ). REE in skarn systems: a LA-ICP-MS study of garnets from the Crown Jewel gold deposit . Geochimica et Cosmochimica Acta 72 , 185 – 205 . Google Scholar CrossRef Search ADS Gibbs G. V. , Ross N. L. , Cox D. F. , Rosso K. M. , Iversen B. B. , Spackman M. A. ( 2013 ). Bonded Radii and the Contraction of the Electron Density of the Oxygen Atom by Bonded Interactions . The Journal of Physical Chemistry. A 117 , 1632 – 1640 . Google Scholar CrossRef Search ADS Ginibre C. , Wörner G. , Kronz A. ( 2002 ). Minor- and trace-element zoning in plagioclase: implications for magma chamber processes at Parinacota volcano, northern Chile . Contributions to Mineralogy and Petrology 143 , 300 – 315 . Google Scholar CrossRef Search ADS Grapes R. ( 1993 ). Barian mica and distribution of barium in metacherts and quartzo-feldspathic schists, southern Alps, New-Zealand . Mineralogical Magazine 57 , 265 – 272 . Google Scholar CrossRef Search ADS Green E. C. R. , Holland T. J. B. , Powell R. ( 2012 ). A thermodynamic model for silicate melt in CaO–MgO–Al2O3–SiO2 to 50 kbar and 1800°C . Journal of Metamorphic Geology 30 , 579 – 597 . Google Scholar CrossRef Search ADS Hammouda T. , Moine B. , Devidal J. , Vincent C. ( 2009 ). Trace element partitioning during partial melting of carbonated eclogites . Physics of the Earth and Planetary Interiors 174 , 60 – 69 . Google Scholar CrossRef Search ADS Hazen R. M. ( 1976 ). Effects of temperature and pressure on the cell dimension and X-ray temperature factors of periclase . American Mineralogist 61 , 266 – 271 . Hetherington C. J. , Gieré R. , Graeser S. ( 2003 ). Composition of Barium-rich white micas from the Berisal complex, Simplon region, Switzerland . The Canadian Mineralogist 41 , 1281 – 1292 . Google Scholar CrossRef Search ADS Hirschmann M. M. , Ghiorso M. S. ( 1994 ). Activities of nickel, cobalt, and manganese silicates in magmatic liquids and applications to olivine/liquid and to silicate/metal partitioning . Geochimica et Cosmochimica Acta 58 , 4109 – 4126 . Google Scholar CrossRef Search ADS Häkli T. A. , Wright T. L. ( 1967 ). The fractionation of nickel between olivine and augite as a geothermometer . Geochimica et Cosmochimica Acta 31 , 877 – 884 . Google Scholar CrossRef Search ADS Holland T. , Baker J. , Powell R. ( 1998 ). Mixing properties and activity-composition relationships of chlorites in the system MgO-FeO-Al2O3-SiO2-H2O . European Journal of Mineralogy 10 , 395 – 406 . Google Scholar CrossRef Search ADS Holland T. J. B. , Powell R. ( 2011 ). An improved and extended internally consistent thermodynamic dataset for phases of petrological interest, involving a new equation of state for solids . Journal of Metamorphic Geology 29 , 333 – 383 . Google Scholar CrossRef Search ADS Hugh-Jones D. A. , Angel R. J. ( 1994 ). A compressional study of MgSiO3 orthoenstatite up to 8.5 GPa . American Mineralogist 79 , 405 – 410 . Ibhi A. , Nachit H. , El Abia H. ( 2005 ). Titanium and barium incorporation into the phyllosilicate phases: the example of phlogopite-kinoshitalite solid solution . Journal de Physique IV France 123 , 331 – 335 . Google Scholar CrossRef Search ADS Israel S. , Saravanan R. , Srinivasan N. , Mohanlal S. ( 2003 ). An investigation on the bonding in MgO, CaO, SrO and BaO from the MEM electron density distributions . Journal of Physics and Chemistry of Solids 64 , 879 – 886 . Google Scholar CrossRef Search ADS Ivaldi G. , Catti M. , Ferraris G. ( 1988 ). Crystal structure at 25 and 700 degrees C of magnesiochloritoid from a high-pressure assemblage (Monte Rosa) . American Mineralogist 73 , 358 – 364 . Jackson C. R. , Shuster D. L. , Parman S. W. , Smye A. J. ( 2016 ). Noble gas diffusivity hindered by low energy sites in amphibole . Geochimica et Cosmochimica Acta 172 , 65 – 75 . Google Scholar CrossRef Search ADS Joswig W. , Fuess H. , Rothbauer R. , Takeuchi Y. , Mason S. A. ( 1980 ). A neutron diffraction study of a one-layer triclinic chlorite (penninite) . American Mineralogist 65 , 349 – 352 . Kohn M. J. ( 2003 ). Geochemical zoning in metamorphic minerals. In: Holland H. D. , Turekian K. K. (eds.) Treatise on Geochemistry . Oxford : Pergamon , pp. 229 – 261 . Google Scholar CrossRef Search ADS Kudoh Y. , Takeda H. ( 1986 ). Single crystal X-ray diffraction study on the bond compressibility of fayalite, Fe2SiO4 and rutile, TiO2 under high pressure . Physica B+C 139 , 333 – 336 . Google Scholar CrossRef Search ADS Kumagai Y. , Kawamoto T. , Yamamoto J. ( 2014 ). Evolution of carbon dioxide bearing saline fluids in the mantle wedge beneath the Northeast Japan arc . Contributions to Mineralogy and Petrology 168 , 1 – 13 . Google Scholar CrossRef Search ADS Lagache M. , Dujon S. ( 1987 ). Distribution of strontium between plagioclases and 1 molar aqueous chloride solutions at 600°C, 1.5 kbar and 750°C, 2 kbar . Bulletin de Minéralogie 110 , 551 – 561 . Levien L. , Prewitt C. T. ( 1981 ). High-pressure structural study of diopside . American Mineralogist 66 , 315 – 323 . Lewis G. V. , Catlow C. R. A. ( 1985 ). Potential models for ionic oxides . Journal of Physics C: Solid State Physics 18 , 1149 – 1161 . Google Scholar CrossRef Search ADS Libourel G. ( 1999 ). Systematics of calcium partitioning between olivine and silicate melt: implications for melt structure and calcium content of magmatic olivines . Contributions to Mineralogy and Petrology 136 , 63 – 80 . Google Scholar CrossRef Search ADS Lewis J. , Schwarzenbach D. , Flack H. D. ( 1982 ). Electric field gradients and charge density in corundum, α-Al2O3 . Acta Crystallographica Section A 38 , 733 – 739 . Google Scholar CrossRef Search ADS Massonne H. J. , Szpurka Z. ( 1997 ). Thermodynamic properties of white micas on the basis of high-pressure experiments in the systems K2O-MgO-Al2O3-SiO2-H2O and K2O-FeO-Al2O3-SiO2-H2O . Lithos 41 , 229 – 250 . Google Scholar CrossRef Search ADS Matzen A. K. , Baker M. B. , Beckett J. R. , Stolper E. M. ( 2013 ). The temperature and pressure dependence of nickel partitioning between olivine and silicate melt . Journal of Petrology 54 , 2521 – 2545 . Google Scholar CrossRef Search ADS McIntire W. ( 1963 ). Trace element partition coefficients—a review of theory and applications to geology . Geochimica et Cosmochimica Acta 27 , 1209 – 1264 . Google Scholar CrossRef Search ADS Mäder U. K. , Percival J. A. , Berman R. G. ( 1994 ). Thermobarometry of garnet–clinopyroxene–hornblende granulites from the Kapuskasing structural zone . Canadian Journal of Earth Sciences 31 , 1134 – 1145 . Google Scholar CrossRef Search ADS Meagher E. ( 1975 ). The crystal structures of pyrope and grossularite at elevated temperatures . American Mineralogist 60 , 218 – 228 . Mellini M. , Zanazzi P. F. ( 1987 ). Crystal structures of lizardite-1T and lizardite- 2H1 from Coli, Italy . American Mineralogist 72 , 943 – 948 . Miyake M. , Nakamura H. , Kojima H. , Marumo F. ( 1987 ). Cation ordering in Co-Mg olivine solid-solution series . American Mineralogist 72 , 594 – 598 . Morse S. ( 2015 ). Linear partitioning in binary solutions: a review with a novel partitioning array . American Mineralogist 100 , 1021 – 1032 . Google Scholar CrossRef Search ADS Mott N. F. , Littleton M. J. ( 1938 ). Conduction in polar crystals. I. Electrolytic conduction in solid salts . Trans. Faraday Soc 34 , 485 – 499 . Google Scholar CrossRef Search ADS Mysen B. O. ( 1976 ). Partitioning of samarium and nickel between olivine, orthopyroxene, and liquid; preliminary data at 20 kbar and 1025°C . Earth and Planetary Science Letters 31 , 1 – 7 . Google Scholar CrossRef Search ADS Mysen B. O. , Kushiro I. ( 1979 ). Pressure dependence of nickel partitioning between forsterite and aluminous silicate melts . Earth and Planetary Science Letters 42 , 383 – 388 . Google Scholar CrossRef Search ADS Okamura F. P. , Ghose S. , Ohashi H. ( 1974 ). Structure and crystal chemistry of calcium Tschermak’s Pyroxene, CaAlAlSiO6 . American Mineralogist 59 , 549 – 557 . Onuma N. , Higuchi H. , Wakita H. , Nagasawa H. ( 1968 ). Trace element partition between two pyroxenes and the host lava . Earth and Planetary Science Letters 5 , 47 – 51 . Google Scholar CrossRef Search ADS Orlando R. , Dovesi R. , Roetti C. , Saunders V. ( 1994 ). Convergence properties of the cluster model in the study of local perturbations in ionic systems. The case of bulk defects in MgO . Chemical Physics Letters 228 , 225 – 232 . Google Scholar CrossRef Search ADS Pabst S. , Zack T. , Savov I. P. , Ludwig T. , Rost D. , Tonarini S. , Vicenzi E. P. ( 2012 ). The fate of subducted oceanic slabs in the shallow mantle: insights from boron isotopes and light element composition of metasomatized blueschists from the Mariana forearc . Lithos 132–133 , 162 – 179 . Google Scholar CrossRef Search ADS Pacalo R. , Graham E. ( 1991 ). Pressure and temperature dependence of the elastic properties of synthetic MnO . Physics and Chemistry of Minerals 18 , 69 – 80 . Google Scholar CrossRef Search ADS Palin E. , Dove M. ( 2004 ). Investigation of Al/Si ordering in tetrahedral phyllosilicate sheets by Monte Carlo simulation . American Mineralogist 89 , 176 – 184 . Google Scholar CrossRef Search ADS Palin E. , Dove M. , Redfern S. , Bosenick A. , Sainz-Diaz C. , Warren M. ( 2001 ). Computational study of tetrahedral Al-Si ordering in muscovite . Physics and Chemistry of Minerals 28 , 534 – 544 . Google Scholar CrossRef Search ADS Parra T. , Vidal O. , Agard P. ( 2002 ). A thermodynamic model for Fe-Mg dioctahedral K white micas using data from phase-equilibrium experiments and natural pelitic assemblages . Contributions to Mineralogy and Petrology 143 , 706 – 732 . Google Scholar CrossRef Search ADS Pickles J. R. , Blundy J. D. , Brooker R. A. ( 2016 ). Trace element thermometry of garnet-clinopyroxene pairs . American Mineralogist 101 , 1438 – 1450 . Google Scholar CrossRef Search ADS Plunder A. , Agard P. , Dubacq B. , Chopin C. , Bellanger M. ( 2012 ). How continuous and precise is the record of P–T paths? Insights from combined thermobarometry and thermodynamic modelling into subduction dynamics (Schistes Lustrés, W. Alps) . Journal of Metamorphic Geology 30 , 323 – 346 . Google Scholar CrossRef Search ADS Purton J. , Allan N. , Blundy J. , Wasserman E. ( 1996 ). Isovalent trace element partitioning between minerals and melts: A computer simulation study . Geochimica et Cosmochimica Acta 60 , 4977 – 4987 . Google Scholar CrossRef Search ADS Putirka K. D. ( 2008 ). Thermometers and barometers for volcanic systems . Reviews in Mineralogy and Geochemistry 69 , 61 – 120 . Google Scholar CrossRef Search ADS Prewitt C. T. , Burnham C. W. ( 1966 ). The crystal structure of jadeite, NaAlSi2O6 . American Mineralogist 51 , 956 – 975 . Rajamani V. , Brown G. E. , Prewitt C. T. ( 1975 ). Cation ordering in Ni-Mg olivine . American Mineralogist 60 , 292 – 299 . Redhammer G. J. , Roth G. ( 2002 ). Single-crystal structure refinements and crystal chemistry of synthetic trioctahedral micas KM3(Al3+, Si4+)4O10(OH)2, where M = Ni2+, Mg2+, Co2+, Fe2+, or Al3+ . American Mineralogist 87 , 1464 – 1476 . Google Scholar CrossRef Search ADS Redfern S. A. T. , Artioli G. , Rinaldi R. , Henderson C. M. B. , Knight K. S. , Wood B. J. ( 2000 ). Octahedral cation ordering in olivine at high temperature. II: an in-situ neutron powder diffraction study on synthetic MgFeSiO4 (Fa50) . Physics and Chemistry of Minerals 27 , 630 – 637 . Google Scholar CrossRef Search ADS Richardson S. M. , Richardson J. W. ( 1982 ). Crystal structure of a pink muscovite from Archer’s Post, Kenya; implications for reverse pleochroism in dioctahedral micas . American Mineralogist 67 , 69 – 75 . Rinaldi R. , Artioli G. , Wilson C. C. , McIntyre G. ( 2000 ). Octahedral cation ordering in olivine at high temperature. I: inthinspacesitu neutron single-crystal diffraction studies on natural mantle olivines (Fa12 and Fa10) . Physics and Chemistry of Minerals 27 , 623 – 629 . Google Scholar CrossRef Search ADS Rodehorst U. , Geiger C. A. , Armbruster T. ( 2002 ). The crystal structures of grossular and spessartine between 100 and 600 K and the crystal chemistry of grossular-spessartine solid solutions . American Mineralogist 87 , 542 – 549 . Google Scholar CrossRef Search ADS Rubatto D. , Hermann J. ( 2003 ). Zircon formation during fluid circulation in eclogites (Monviso, Western Alps): implications for Zr and Hf budget in subduction zones . Geochimica et Cosmochimica Acta 67 , 2173 – 2187 . Google Scholar CrossRef Search ADS Sahoo D. ( 1984 ). Elastic continuum theories of lattice defects: a review . Bulletin of Materials Science 6 , 775 – 798 . Google Scholar CrossRef Search ADS Sainz-Diaz C. I. , Hernández-Laguna A. , Dove M. T. ( 2001a ). Modeling of dioctahedral 2: 1 phyllosilicates by means of transferable empirical potentials . Physics and Chemistry of Minerals 28 , 130 – 141 . Google Scholar CrossRef Search ADS Sainz-Diaz I. C. , Hernández-Laguna A. , Dove T. M. ( 2001b ). Theoretical modelling of cis-vacant and trans-vacant configurations in the octahedral sheet of illites and smectites . Physics and Chemistry of Minerals 28 , 322 – 331 . Google Scholar CrossRef Search ADS Salje E. K. H. ( 1995 ). Chemical mixing and structural phase transitions: the plateau effect and oscillatory zoning near surfaces and interfaces . European Journal of Mineralogy 7 , 791 – 806 . Google Scholar CrossRef Search ADS Sasaki S. , Fujino K. , Takéuchi Y. ( 1979 ). X-ray determination of electron-density distributions in oxides, MgO, MnO, CoO, and NiO, and atomic scattering factors of their constituent atoms . Proceedings of the Japan Academy, Series B 55 , 43 – 48 . Google Scholar CrossRef Search ADS Scambelluri M. , Pettke T. , Cannaò E. ( 2015 ). Fluid-related inclusions in Alpine high-pressure peridotite reveal trace element recycling during subduction-zone dehydration of serpentinized mantle (Cima di Gagnone, Swiss Alps) . Earth and Planetary Science Letters 429 , 45 – 59 . Google Scholar CrossRef Search ADS Schleid T. , Meyer G. ( 1989 ). Single crystals of rare earth oxides from reducing halide melts . Journal of the Less Common Metals 149 , 73 – 80 . Google Scholar CrossRef Search ADS Schreiber H. D. ( 1979 ). Experimental studies of nickel and chromium partitioning into olivine from synthetic basaltic melts. Proceedings of the Lunar and Planetary Science Conference 10th. pp. 509 – 516 . Shannon R. D. ( 1976 ). Revised effective ionic radii and systematic studies of interatomic distances in halides and chalcogenides . Acta Crystallographica Section A 32 , 751 – 767 . Google Scholar CrossRef Search ADS Sharp Z. D. , Hazen R. , Finger L. W. ( 1987 ). High-pressure crystal chemistry of monticellite CaMgSiO4 . American Mineralogist 72 , 748 – 755 . Shimizu N. ( 1975 ). Rare earth elements in garnets and clinopyroxenes from garnet lherzolite nodules in kimberlites . Earth and Planetary Science Letters 25 , 26 – 32 . Google Scholar CrossRef Search ADS Smye A. J. , Warren C. J. , Bickle M. J. ( 2013 ). The signature of devolatisation: Extraneous 40Ar systematics in high-pressure metamorphic rocks . Geochimica et Cosmochimica Acta 113 , 94 – 112 . Google Scholar CrossRef Search ADS Soga N. , Anderson O. L. ( 1966 ). High-temperature elastic properties of polycrystalline MgO and Al2O3 . Journal of the American Ceramic Society 49 , 355 – 359 . Google Scholar CrossRef Search ADS Spear F. ( 1993 ). Metamorphic Phase Equilibria and Pressure-Temperature-Time Paths . Mineralogical Society of America , 799 p. ISBN 0-939950-34-0. Thangadurai V. , Weppner W. ( 2005 ). Li6ALa2Ta2O12 (A = Sr, Ba): novel garnet-like oxides for fast lithium ion conduction . Advanced Functional Materials 15 , 107 – 112 . Google Scholar CrossRef Search ADS Tipper E. T. , Calmels D. , Gaillardet J. , Louvat P. , Capmas F. , Dubacq B. ( 2012 ). Positive correlation between Li and Mg isotope ratios in the river waters of the Mackenzie Basin challenges the interpretation of apparent isotopic fractionation during weathering . Earth and Planetary Science Letters 333-334 , 35 – 45 . Google Scholar CrossRef Search ADS Tuff J. , Gibson S. A. ( 2007 ). Trace-element partitioning between garnet, clinopyroxene and Fe-rich picritic melts at 3 to 7 GPa . Contributions to Mineralogy and Petrology 153 , 369 – 387 . Google Scholar CrossRef Search ADS van Westrenen W. , Draper D. S. ( 2007 ). Quantifying garnet-melt trace element partitioning using lattice-strain theory: new crystal-chemical and thermodynamic constraints . Contributions to Mineralogy and Petrology 154 , 717 – 730 . Google Scholar CrossRef Search ADS van Westrenen W. , Blundy J. , Wood B. ( 1999 ). Crystal-chemical controls on trace element partitioning between garnet and anhydrous silicate melt . American Mineralogist 84 , 838 – 847 . Google Scholar CrossRef Search ADS van Westrenen W. , Allan N. , Blundy J. , Purton J. , Wood B. ( 2000 ). Atomistic simulation of trace element incorporation into garnets - comparison with experimental garnet-melt partitioning data . Geochimica et Cosmochimica Acta 64 , 1629 – 1639 . Google Scholar CrossRef Search ADS van Westrenen W. , Allan N. L. , Blundy J. D. , Lavrentiev M. Y. , Lucas B. R. , Purton J. A. ( 2003 ). Trace element incorporation into pyrope-grossular solid solutions: an atomistic simulation study . Physics and Chemistry of Minerals 30 , 217 – 229 . Google Scholar CrossRef Search ADS Verlaguet A. , Brunet F. , Goffé B. , Murphy W. M. ( 2006 ). Experimental study and modeling of fluid reaction paths in the quartz - kyanite ±muscovite–water system at 0.7 GPa in the 350-550°C range: Implications for Al selective transfer during metamorphism . Geochimica et Cosmochimica Acta 70 , 1772 – 1788 . Google Scholar CrossRef Search ADS Vidal O. , Dubacq B. ( 2009 ). Thermodynamic modelling of clay dehydration, stability and compositional evolution with temperature, pressure and H2O activity . Geochimica et Cosmochimica Acta 73 , 6544 – 6564 . Google Scholar CrossRef Search ADS Vidal O. , Parra T. ( 2000 ). Exhumation paths of high-pressure metapelites obtained from local equilibria for chlorite–phengite assemblages . Geological Journal 35 , 139 – 161 . Google Scholar CrossRef Search ADS Vidal O. , De Andrade V. , Lewin E. , Munoz M. , Parra T. , Pascarelli S. ( 2006 ). P-T-deformation-Fe3+/Fe2+ mapping at the thin section scale and comparison with XANES mapping: application to a garnet-bearing metapelite from the Sambagawa metamorphic belt (Japan) . Journal of Metamorphic Geology 24 , 669 – 683 . Google Scholar CrossRef Search ADS Vieillard P. ( 1994a ). Prediction of enthalpy of formation based on refined crystal structures of multisite compounds: Part 2. Application to minerals belonging to the system Li2O-Na2O-K2O-BeO-MgO-CaO-MnO-FeO-Fe2O3-Al2O3-SiO2-H2O. Results and discussion . Geochimica et Cosmochimica Acta 58 , 4065 – 4092 , 4095–4107. Google Scholar CrossRef Search ADS Vieillard P. ( 1994b ). Prediction of enthalpy of formation based on refined crystal structures of multisite compounds: Part 1. Theories and examples . Geochimica et Cosmochimica Acta 58 , 4049 – 4063 . Google Scholar CrossRef Search ADS Vinograd V. L. , Sluiter M. H. ( 2006 ). Thermodynamics of mixing in pyrope-grossular, Mg3Al2Si3O12-Ca3Al2Si3O12, solid solution from lattice dynamics calculations and Monte Carlo simulations . American Mineralogist 91 , 1815 – 1830 . Google Scholar CrossRef Search ADS Watson E. B. ( 1979 ). Calcium content of forsterite coexisting with silicate liquid in the system Na2O-CaO-MgO-Al2O3-SiO2 . American Mineralogist 64 , 824 – 829 . Welch M. D. , Marshall W. G. ( 2001 ). High-pressure behavior of clinochlore . American Mineralogist 86 , 1380 – 1386 . Google Scholar CrossRef Search ADS Whitney D. L. , Evans B. W. ( 2010 ). Abbreviations for names of rock-forming minerals . American Mineralogist 95 , 185 – 187 . Google Scholar CrossRef Search ADS Wood B. J. ( 1974 ). The solubility of alumina in orthopyroxene coexisting with garnet . Contributions to Mineralogy and Petrology 46 , 1 – 15 . Google Scholar CrossRef Search ADS Wood B. J. , Blundy J. D. ( 1997 ). A predictive model for rare earth element partitioning between clinopyroxene and anhydrous silicate melt . Contributions to Mineralogy and Petrology 129 , 166 – 181 . Google Scholar CrossRef Search ADS Wood B. J. , Blundy J. D. ( 2001 ). The effect of cation charge on crystal-melt partitioning of trace elements . Earth and Planetary Science Letters 188 , 59 – 71 . Google Scholar CrossRef Search ADS Wood B. J. , Nicholls J. ( 1978 ). The thermodynamic properties of reciprocal solid solutions . Contributions to Mineralogy and Petrology 66 , 389 – 400 . Google Scholar CrossRef Search ADS Wyckoff R. ( 1963 ). Crystal Structures - Volume 1 . Interscience Publishers . Zack T. , Foley S. F. , Rivers T. ( 2002 ). Equilibrium and disequilibrium trace element partitioning in hydrous eclogites (Trescolmen, Central Alps) . Journal of Petrology 43 , 1947 – 1974 . Google Scholar CrossRef Search ADS Zanazzi P. , Montagnoli M. , Nazzareni S. , Comodi P. ( 2006 ). Structural effects of pressure on triclinic chlorite: A single-crystal study . American Mineralogist 91 , 1871 – 1878 . Google Scholar CrossRef Search ADS © The Author(s) 2018. Published by Oxford University Press. All rights reserved. For permissions, please e-mail: journals.permissions@oup.com This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/about_us/legal/notices)
Journal
Journal of Petrology – Oxford University Press
Published: Mar 7, 2018