A computer simulation study of natural silicate melts. Part II: High pressure properties

Geochimica et Cosmochimica Acta 71 (2007) 4538–4556 www.elsevier.com/locate/gca

A computer simulation study of natural silicate melts. Part II: High pressure properties Bertrand Guillot *, Nicolas Sator Laboratoire de Physique The´orique de la Matie`re Condense´e, Universite´ Pierre et Marie Curie (Paris 6), UMR CNRS 7600, 4 place Jussieu, 75252 Paris Cedex 05, France Received 20 November 2006; accepted in revised form 21 May 2007; available online 17 July 2007

Abstract The thermodynamic, structural and transport properties of natural silicate melts under pressure are investigated by molecular dynamics simulation with the help of a force field recently introduced by us [Guillot B. and Sator N. (2007) A computer simulation study of natural silicate melts. Part I: low pressure properties. Geochim. Cosmochim. Acta 71, 1249–1265]. It is shown that the simulation reproduces accurately the bulk moduli of a large variety of silicate liquids as evaluated from ultrasonic studies. The equations of state (EOS) of the simulated melts are in good agreement with the density data on mid-ocean ridge basalt, komatiite, peridotite and fayalite as obtained either by sink/float experiments or by shock-wave compression. From the structural point of view it is shown that the population of [5]Al and [6]Al species increases rapidly upon initial compression (0–50 kbar) whereas for Si these highly coordinated species are found in a significant abundance (>5%) only above 50 kbar for [5]Si and 100–150 kbar for [6]Si. This increase of the coordination of network formers is not the only response of the melt structure to the densification: there is also a large redistribution of the T–O–T (T = Si, Al) bond angles with the pressure and noticeably upon initial compression in rhyolitic and basaltic liquids. Furthermore, a detailed analysis of the population of bridging oxygens (BO) and nonbridging oxygens (NBO) points out that the polymerization of the melt generally increases when the pressure increases, the magnitude of this polymerization enhancement being all the more important that the melt is depolymerized at low pressure. The role of triclusters (threefold coordinated oxygens to network former cations) is particularly emphasized in acidic and basaltic liquids. The pressure-induced redistribution of the oxygen atoms through the melt structure is also stressed. Finally, the simulation predicts a nonmonotonic behavior of the diffusivity of network former ions when the pressure increases, a feature with depends on the melt composition. This could have a counterpart in the electrical conductivity at sufficiently high temperature when the viscosity of the liquid is low.  2007 Elsevier Ltd. All rights reserved.

1. INTRODUCTION A precise knowledge of the thermophysical and structural properties of natural silicate melts is a prerequisite for understanding many processes at work in the earth’s interior. Since melts are generated at depth, it is essential to know how their physical properties change with pressure and temperature. Thus, to estimate the buoyancy of melts generated at depth, to calculate the mineral-melt partition

*

Corresponding author. Fax: +33 (0)1 44 27 51 00. E-mail address: [email protected] (B. Guillot).

0016-7037/$ - see front matter  2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.gca.2007.05.029

coefficients of elements and volatiles or to evaluate the solubility of rare gases in melts at high pressure, an accurate evaluation of the equation of state (EOS) for liquid silicates is needed to make reliable predictions. However density measurements of melts under pressure remain somewhat scarce (Fujii and Kushiro, 1977; Rigden et al., 1984; Agee and Walker, 1988; Miller et al., 1991; Agee and Walker, 1993; Ohtani et al., 1993; Suzuki et al., 1995, 1998; Circone and Agee, 1996; Smith and Agee, 1997; Ohtani et al., 1998; Ohtani and Maeda, 2001; Chen et al., 2002; Suzuki and Ohtani, 2003) and mostly incomplete (in general only one or two state points are evaluated) which prevents the accurate determination of the EOS. Moreover in these studies,

Simulation study of natural silicate melts at high pressure

the experimental set up deals either with isostatic pressure (sink/float experiment) or with shock-wave compression (Hugoniot curve) and the link between the two methods is not straightforward; this is additionally obscured by experimental uncertainties. Correlatively, structural information on compression mechanisms in various silicate glasses quenched from melts at high pressure (up to 100 kbar) have been extracted from Raman, infrared, XANES and RMN spectroscopies (Xue et al., 1991, 1994; Yarger et al., 1995; Poe et al., 2001; Lee et al., 2003, 2004; Allwardt et al., 2005a; Lee, 2004, 2005; Lee et al., 2006). However, the data obtained from quenched melts give only a first order approximation of what really happens in the melt under pressure at superliquidus temperature (Allwardt et al., 2004, 2005b). So it can be useful to have a theoretical guide which gives unity to the thermodynamic and structural properties of silicate melts and is consistent with data of the literature. Our purpose is to tackle by computer simulation the evaluation of the thermophysical and structural properties of natural silicate melts under pressure. Our final goal is to provide a theoretical tool to investigate in the future some of the key processes taking place in melting regions at depth (e.g. solubility of volatiles in magmas and mineral-melt partitioning). We have recently proposed (Guillot and Sator, 2007) a simple atom–atom interaction potential to describe by molecular simulation the system involving the nine major oxides K2O–Na2O–CaO–MgO–FeO– Fe2O3–Al2O3–TiO2–SiO2 (KNCMFATS) present in common igneous rocks. Through use of the molecular dynamics method (MD) we have shown that this empirical force field is able to reproduce realistically a number of thermodynamic, structural and transport properties of a large variety of natural silicate melts from acidic to ultrabasic compositions at atmospheric pressure. Here we investigate the EOS (pressure versus density), the compression mechanisms (ionic coordinations, bond angle distributions, bonding and nonbonding oxygens) and the transport properties (diffusivities and electrical conductivity) under pressure of a selected set of molten rocks (namely, rhyolite, andesite, basalts, peridotite, komatiite, olivine and fayalite). 2. METHODOLOGY AND COMPUTATIONAL DETAILS In a previous study, we proposed (Guillot and Sator, 2007) a simple interionic potential to implement into a MD code to simulate natural silicate melts belonging to the KNCFMATS system. We recall briefly that the MD method consists of solving iteratively the equations of motion of an assembly of particles interacting through a given force field. The knowledge of the position and the velocity of each particle with time enables one to evaluate any macroscopic property by performing a time average of the corresponding quantity along the MD run (for further details see Allen and Tildesley, 1987). In this framework, the interaction energy between atoms i and j (where i, j = Si, Ti, Al, Fe3+, Fe2+, Mg, Ca, Na, K, O) is given by vðrij Þ ¼ zi zj =rij þ Bij erij =qij  C ij =r6ij

ð1Þ

4539

where rij is the distance between atoms i and j, zi is the effective charge associated with the atom i, and where Bij, qij and Cij are parameters describing repulsive and dispersive forces between the pair ij. To ensure the transferability of the interaction potential with the melt composition, the valence (charge) assigned to the atoms is kept fixed for all compositions (except for molten fayalite for which a specific set of parameters was defined to improve the description of its properties). The potential parameters (see Table 1) were obtained by constraining the MD simulations to reproduce at best the density at atmospheric pressure and the microscopic structure of a selected set of 11 natural silicates in the liquid phase (Guillot and Sator, 2007). It should be noted that for the system (CaO–MgO–Al2O3–SiO2) our force field becomes identical to that developed by Matsui (1994). This is important since Matsui has shown that his force field is transferable and reproduces quite well the symmetries, lattice parameters, cation–anion distances and bulk moduli of 27 crystals including oxides (periclase, lime, corundum, spinel, silica polymorphs), Mg meta and ortho-silicates, Al garnets and various Ca and Al bearing silicates (Matsui, 1996a). Moreover the structure and the EOS parameters of molten enstatite, wollastonite, diopside and anorthite deduced from the MD simulations compare well with experimental data (Matsui, 1996b). The present work therefore, represents an extension of the above studies to a wider compositional range in the KNCMFATS system. Here we evaluate the liquid properties under pressure of 10 compositions of common igneous rocks, namely a rhyolite (obsidian), an andesite (from Unzen volcano), a diopside–anorthite mixture (0.36:0.64), a mid-ocean ridge basalt (Kane zone of central Atlantic), a high-Ti lunar basalt (Apollo 14 black glass), a Fe-rich Martian basalt (Gusev crater), a peridotite (PHN1611), a komatiite (Belingwe), an olivine (San Carlos) and fayalite. The corresponding chemical compositions are reported in Table 2 and some details on these natural rocks are given in Appendix of our previous article (Guillot and Sator, 2007). The MD calculations were performed with the code DL_POLY_2 of the CCLRC Program library (Smith and Forester, 1996). The equations of motion are solved with the Verlet leap-frog algorithm (time step of 1 femtosecond) and we use a cubic simulation box with periodic boundary conditions in 3D. The long range coulombic interactions are accounted for by a Ewald sum (with aL = 5–7 where a is the width of the charge distribution on each ion and L is the size of the simulation box). A system size fixed to 1000 atoms to simulate different melts is a good compromise between accuracy and computational cost. As in our previous study, several test runs were performed with 8000 atoms to ensure that no significant size effects were detected for the thermodynamic properties and the liquid structure (except for the pffiffiffiffi expected 1= N reduction of the statistical noise). Long simulation runs (1 ns i.e. 106 time steps) in the microcanonical (NVE) and isobaric–isothermal (NPT) ensembles were crosschecked to evaluate the statistical uncertainties on PVT properties and very long simulation runs (10 ns i.e. 107 time steps) were performed when evaluating time dependent properties (self diffusion coefficients and electrical conductivity). Thus with N  1000 atoms and for all

4540

B. Guillot, N. Sator / Geochimica et Cosmochimica Acta 71 (2007) 4538–4556

Table 1 Potential parameters

O Si Ti Al Fe3+ Fe2+ Mg Ca Na K

z (e)

B (kJ mol1)

˚) q (A

C (A6 kJ mol1)

0.945 (1.20) 1.89 (2.40) 1.89 1.4175 1.4175 0.945 (1.20) 0.945 0.945 0.4725 0.4725

870570.0 (889916.0) 4853815.5 (5900530.0) 4836495.0 2753544.3 773840.0 1257488.6 (1450950.0) 3150507.4 15019679.1 11607587.5 220447.4

0.265 (0.265) 0.161 (0.161) 0.178 0.172 0.190 0.190 (0.190) 0.178 0.178 0.170 0.290

8210.17 (8210.17) 4467.07 (4467.07) 4467.07 3336.26 0.0 0.0 (0.0) 2632.22 4077.45 0.0 0.0

Values of B, q and C given in column 3, 4 and 5, respectively, correspond to oxygen–oxygen (first row) and cation–oxygen interactions, the cation–cation interactions being described only by coulombic repulsive forces. Values given in parenthesis for O, Si and Fe2+ are solely for fayalite.

Table 2 Chemical compositions (weight fraction) of the silicate melts simulated in this study Silicate

SiO2 (wt%)

TiO2 (wt%)

Al2O3 (wt%)

Rhyolite Andesite An0.36Di0.64 MORB Mars basalt LG14 Komatiite Peridotite Olivine Fayalite

74.51 56.65 50.34 50.59 47.68 34.00 46.73 45.10 40.68 29.49

0.10 (0) 1.01 (3) 0.0 (0) 1.52 (4) 0.54 (1) 16.40 (50) 0.31 (1)

13.25 17.41 15.38 15.11 10.96 4.60 6.30 2.80 0.01

(257) (203) (180) (185) (176) (136) (168) (159) (142) (143)

(54) (73) (66) (65) (48) (22) (27) (12) (0)

Fe2O3 (wt%)

FeO (wt%)

MgO (wt%)

0.32 (1) 4.63 (12) 0.0 (0) 1.15 (3) 3.09 (9)

1.28 (3) 3.53 (11) 0.0 (0) 8.39 (26) 15.82 (49) 24.50 (83) 10.76 (32) 10.40 (31) 8.76 (25) 70.51 (286)

0.08 4.30 10.80 7.77 12.62 13.30 28.42 38.40 50.52

CaO (wt%) (0) (23) (58) (42) (69) (79) (152) (203) (262)

0.75 7.38 23.48 11.87 7.96 6.90 6.29 3.40 0.06

(3) (28) (90) (47) (31) (30) (24) (13) (0)

Na2O (wt%)

K2O (wt%)

Total

4.15 (28) 3.23 (22) 0.0 (0) 2.94 (21) 2.68 (19) 0.23 (2) 0.85 (6)

5.64 (25) 1.56 (7) 0.0 (0) 0.13 (1) 0.06 (0) 0.16 (0) 0.13 (1)

100.08 99.70 100.00 99.47 101.41 100.09 99.79 100.10 100.03 100.00

(1000) (998) (1001) (1000) (1000) (1000) (1001) (1001) (1000) (1001)

In parenthesis are the numbers of cations of each species used in the simulation of the corresponding silicate melt. In the last column are indicated the total weight fraction and in parenthesis the total number of atoms (cations + oxygens) used in the simulation.

investigated compositions the statistical fluctuations are: ±5 kbar for the pressure, DT/T = ±1% and Dq/q = ±1%. Further details will be given in the text when necessary. 3. THERMODYNAMIC PROPERTIES OF THE SILICATE MELTS UNDER PRESSURE 3.1. Thermal expansion under pressure We have first evaluated the evolution of the molar volume with temperature under isobaric conditions. For each composition of the melt, several samples were prepared at high temperature (>2000 K) in the liquid phase for various pressures between zero pressure and 300 kbar. After equilibration, the samples were cooled down at constant pressure to low temperature (400 K) using a conservative cooling rate of 2.1011 K/s, this latter value being typical of cooling rates used in computer simulation studies of the literature (e.g. Vollmayr et al., 1996; Caprion and Schober, 2002; Guillot and Guissani, 2003). Thus, when the sample is cooled down, at some critical temperature (composition dependent) the simulated melt experiences a liquid–glass transition characterized by a kinetic arrest of the atoms and a change of slope of the molar volume (and of the enthalpy as well). As discussed in our previous

paper (Guillot and Sator, 2007), the accurate evaluation of the glass transition temperature by simulation is very costly in computer time and is beyond the scope of our study. However, in using different sizes of the system (N = 1000 and 8000 atoms) and by evaluating the self diffusion coefficients of network former ions at different temperatures along various isobars, we have carefully checked that the evolution with temperature of the thermodynamic properties of the investigated melts describes the fully relaxed liquid phase as long as T > 1300 K (the latter value being conservative for most of the simulated samples). Below this temperature and depending on the composition, the thermodynamic properties become sensitive to the size of the sample and to the cooling rate, meaning that the melt is entering into the glassy phase. For illustration it is shown in Fig. 1 the temperature evolution of the molar volume of our simulated MORB at various pressures. It is clear that within the statistical noise of the simulation (±1% on V), the molar volume of the liquid at isobaric conditions varies linearly with the temperature. This feature occurs with all melts investigated here (not shown). This is in accordance with recent measurements made at atmospheric pressure for calcium aluminosilicate and natural melts at superliquidus temperatures (Potuzak et al., 2006; Potuzak and Dingwell, 2006). We have reported

Simulation study of natural silicate melts at high pressure

4541

predicted by the EOS of Ghiorso and the one obtained by simulation do not deviate from each other by more than 2%. Above 100 kbar the deviation with the EOS of Ghiorso increases (e.g. 5% at 200 kbar), likely because the latter one is less constrained in this pressure range (see also the discussion in Section 3.2). More generally, since high pressure data on the expansivity of natural melts is scarce, it is noteworthy that Suzuki et al. (1998) for a FeO-rich peridotite evaluated the expansivity around 9.6 ± 3.1 · 104 cm3 mol1 K1 at 81 kbar and 2253 K, and around 8.6 ± 2.6 · 104 cm3 mol1 K1 at 151 kbar and 2603 K, as compared with 9.4 · 104 cm3 mol1 K1 at 100 kbar and 7.2 · 104 cm3 mol1 K1 at 200 kbar for our molten peridotite (PHN1611) in the same temperature range. 3.2. Bulk modulus and density of melt at high pressure

Fig. 1. Molar volume as function of temperature for molten MORB at various pressures. The wiggling lines are the simulation data, the triangles are the experimental data for Kilauea basalt (Nelson and Carmichael, 1979), the dots are for Slapany basalt (Potuzak and Dingwell, 2006) and the dotted lines are the predictions for MORB given by the multicomponent mixing model of Ghiorso and Kress (2004).

in Fig. 1 the 1-bar measurements of Nelson and Carmichael (1979) for a Kilauea basalt and that of Potuzak and Dingwell (2006) for a Slapany basalt whose the compositions are close to our MORB. Although the overall agreement between simulated and experimental data for the molar volume is excellent, the Slapany basalt exhibits a molar expansivity of 13.1 ± 2.0 · 104 cm3 mol1 K1 when that of Kilauea basalt amounts to 26.9 · 104 cm3 mol1 K1 (with a large uncertainty as high as 50%) and 25.5 ± 3.5 · 104 cm3 mol1 K1 for our simulated MORB. In this context, it is useful to compare our results at zero pressure with the values of the expansivity deduced from the multicomponent mixing model of Ghiorso and Kress (2004) which is based upon a large compilation of data on silicate liquids at atmospheric pressure (see the dotted line for P = 0 in Fig. 1). The model of Ghiorso and Kress predicts a slightly larger molar volume for MORB (by 2% at 1673 K) and a slightly smaller molar expansivity (19.0 · 104 cm3 mol1 K1 instead of 25.5 ± 3.5 · 104 cm3 mol1 K1). However, when applied to Slapany basalt, the model of Ghiorso and Kress underestimates its molar volume by 2% at 1673 K and overestimates the molar expansivity by about 50% (19.0 instead of 13.1 · 104 cm3 mol1 K1). Under pressure our simulation data exhibit (for all composition) a net decrease of (oV/oT)P with increasing pressure (see Fig. 1) whereas the EOS of Ghiorso (2004) predicts a molar expansivity essentially invariant with pressure. Nevertheless factually, over a large P–T domain (1473–2000 K and 0–100 kbar, see Fig. 1) the molar volume

We have performed a series of MD simulations in the isothermal–isobaric ensemble to calculate the evolution of the density with the pressure along several isotherms in the liquid range. In practice, each sample of a given composition was equilibrated around zero pressure (in fact P = 0 ± 5 kbar) at the temperature of investigation (generally between 1673 and 2673 K with a temperature fluctuation DT/T  1%) and next was progressively compressed at a rate of 0.1 kbar/ps up to very high pressure (600 kbar or more). By performing independent runs in the microcanonical (NVE) ensemble, we have carefully checked that this compression rate was sufficiently low for the sample to be in quasi equilibrium all along the isotherm. The data points (PVT) collected at the end of the simulations have been fitted with the help of equations of the literature, namely the commonly used Birch–Murnaghan equation (Birch, 1978), and a related equation due to Vinet (Vinet et al., 1989). These equations are     P BM ¼ ð3=2ÞK T X 7  X 5 f1  ð3=4Þ 4  K 0T ðX 2  1Þg ð2Þ ð1  X Þ ð3=2ÞðK 0 1Þð1X Þ T P Vinet ¼ 3K T e ð3Þ X2 where the constants KT and K 0T are the isothermal bulk modulus at 1 bar and its pressure derivative, X = (V0(P  0, T)/V(P, T))1/3 with V0(P  0, T) being the molar volume at one bar (i.e. around zero pressure for the simulation). To illustrate the ability of these two Eqs. (2 and 3) to fit the simulation data, it is shown in Fig. 2 the evolution of pressure versus density for our simulated peridotite along the 2673 K isotherm. The Birch–Murnaghan equation (BMEOS) provides an excellent description up to 350 kbar but it overestimates the compressibility of the simulated melt at higher pressures. The Vinet equation leads to a fit virtually undistinguishable from that given by BMEOS. However none of these equations are able to fit the simulation data over the entire pressure range investigated here (0–1500 kbar) and this is observed with all compositions (not shown). The conclusion is that our simulation data can be fitted accurately by the BMEOS and the Vinet EOS only in the pressure range 0–350 kbar. Values of the parameters for BMEOS obtained from the fits of our simulation data for the investigated melts are given in Table 3.

4542

B. Guillot, N. Sator / Geochimica et Cosmochimica Acta 71 (2007) 4538–4556 Table 3 Birch–Murnaghan equation parameters for simulated melts along various isotherms

Fig. 2. Pressure versus density for molten peridotite at 2673 K. The circles are the simulation data, the continuous curve represents the best fit obtained with the Birch–Murnaghan equation and the dashdotted curve that obtained with the Vinet equation (they are virtually undistinguishable from each other). The inset is an enlarged view of the 0–300 kbar range. The statistical uncertainties are of the order of the size of the circles.

A quick look at the latter indicates that for all compositions, the bulk modulus KT at zero pressure decreases when the temperature increases, while its pressure derivative K 0T is mostly invariant or increases somewhat. Notice that the value of KT resulting from the fitting procedure does not always provide a very accurate evaluation of the bulk modulus of the simulated melt at 1 bar. As a matter of fact, in our calculations the statistical uncertainty on the pressure (DP  ±5 kbar) is greater around zero pressure that at high pressure. In order to compare our results with the literature, in Table 4 are reported literature values of densities and bulk moduli (isothermal or adiabatic) at atmospheric pressure measured for analog melts or obtained from the multicomponent mixing models of Lange and Carmichael (1987) and Ghiorso and Kress (2004). As far as the density is concerned, the deviation between simulated values and measured (or predicted) ones is at most 2% and generally 1% or less. As for the bulk moduli, the agreement between our results and that obtained from ultrasonic studies is quite impressive considering the various sources of uncertainties present together in the simulations and in the experimental data. This excellent agreement over a large composition range is an important test of the validity of our calculations. An interesting feature exhibited by the simulation data is that for all compositions the bulk modulus decreases when the temperature increases. This is in accordance with the general trend shown by the literature data for silicate melts of various compositions (Rivers and Carmichael, 1987; Webb and Dingwell, 1994; Webb and Courtial, 1996). For instance over a limited temperature range (100–150 K) above the liquidus the measurements of Rivers and Carmichael (1987) lead to

T (K)

q0 (g/cm3)

KT (kbar)

K 0T

Rhyolite

1673 2073 2473

2.28 2.27 2.20

125 105 90

7.0 6.9 7.5

Andesite

2273

2.47

135

7.1

MORB

1673 2073 2473

2.69 2.59 2.48

175 155 135

7.2 7.2 7.2

BMars

1673 2073 2473

2.76 2.64 2.52

166 130 116

7.5 7.7 7.7

LG14

1673 2353

3.05 2.76

225 135

7.3 8.0

Komatiite

2273 2673 3073

2.58 2.44 2.32

160 115 95

7.0 8.0 8.2

Peridotite

2273 2673 3073

2.61 2.48 2.33

165 129 95

7.2 7.9 8.4

Olivine

2273 2673 3073

2.61 2.46 2.295

175 135 89

7.2 7.6 8.7

Fayalite

2273 2673 3073

3.45 3.30 3.10

155 130 87

6.4 6.9 8.2

These parameters fit accurately the simulation results only up to P  350 kbar.

(oKT/oT)0 = 20 · 103 kbar/K for a rhyolitic composition, 30.103 kbar/K for Kilauea basalt and 80.103 kbar/K for fayalite as compared with 50.103 kbar/K, 50.103 kbar/K and 63.103 kbar/K for corresponding simulated melts. This general tendency of the bulk modulus to decrease upon heating suggests that it is inaccurate to neglect the temperature dependence of the bulk modulus when investigating the compressibility curve of a melt well above the liquidus temperature at atmospheric pressure. Moreover some caution must be taken with the use of the multicomponent mixing models which sometimes predict erroneously an increase of KT with temperature (e.g. basalts with the model of Ghiorso and Kress, 2004). High pressure data on silicate melts, natural or otherwise, are scarce and are obtained either from sink/float experiments under isostatic pressure or from shock-wave compression studies. Since these two techniques imply different thermodynamic pathways, we will discuss their results separately. In using the sink/float method, Agee (1998) measured the density of a MORB compressed around 58.5 kbar at liquidus temperature (1950 K, see Yasuda et al., 1994) whereas Ohtani and Maeda (2001) investigated a model MORB compressed to about 150 kbar at 2473 and 2773 K. These results are presented in Fig. 3

Simulation study of natural silicate melts at high pressure

4543

Table 4 Comparison between the bulk moduli (isothermal or adiabatic) determined in the present study by applying the Birch–Murnaghan EOS to high pressure simulation data and one-atmosphere experimental data for analog melts T (K)

q0 (g/cm3)

K0 (kbar)

Rhyolite

1673 1673

2.28 2.30a, 2.31b

125 132c, 130b

This study (KT) Literature data

An0.36Di0.64

1673 1673

2.70 2.63d, 2.62b

235 230e, 225b

This study (KS) Literature data

MORB

1673 1673

2.69 2.67f, 2.65b

175 174g, 169–181h, 190b

This study (KT) Literature data

BMars

1673 1673

2.76 2.74i, 2.72b

200 202j, 200b

This study (KS) Literature data

LG14

1673 1673

3.05 3.07b, 3.09a

225 187b, 192i

This study (KT) Literature data

Komatiite

1823 1823

2.73 2.68k, 2.69b

210 230l, 202b

This study (KS) Literature data

Peridotite

1873 1873

2.75 2.74k, 2.75b

260 260m, 239b

This study (KS) Literature data

Fayalite

1573 1573

3.75 3.75n, 3.77b

220 220o, 211b

This study (KS) Literature data

The latter values are obtained from ultrasonic studies or deduced from multicomponent mixing models. a Lange and Carmichael (1987). b Ghiorso and Kress (2004). c NZC-7 in Rivers and Carmichael (1987). d Knoche et al. (1992). e Secco et al. (1991a). f CRB in Murase and McBirney (1973). g Kil-2 in Rivers and Carmichael (1987). h CMS4 and CMS5 in Manghnani et al. (1986). i Kress and Carmichael (1991). j G26 in Secco et al. (1991b). k Courtial et al. (1997). l ZV-10 in Secco et al. (1991b). m For 38% of MgO in Secco et al. (1991b). n Shiraishi et al. (1978). o Rivers and Carmichael (1987).

together with our simulation results for MORB along two isotherms (2073 and 2473 K). Considering the slight variation in composition between experimental and simulated samples and the experimental uncertainties generated by the sink/float experiments, the agreement is quite good. However above 120 kbar our results deviate significantly (upwards) from the compressibility curve proposed by Agee (1998) whereas a similar situation occurs above 200 kbar with the compressibility curve proposed by Ohtani and Maeda (2001). This feature can be rationalized in the following manner. Agee (1998) and Ohtani and Maeda (2001) neglect the temperature dependence of the isothermal bulk modulus and use for KT the value calculated at 1673 K by the partial molar compressibility of liquid oxides of Lange and Carmichael (1987). The KT values of 197 kbar (Agee, 1998) and 187 kbar (Ohtani and Maeda, 2001) are, therefore, too large to describe the isothermal compressibility curve at 2073 and 2473 K this has the consequence of an underestimation of the value of K 0T in the fitting procedure. In support of our suggestion, note that the bulk modulus of Kilauea basalt is 175 kbar at 1673 K, 161 kbar at 2073 K

and 149 kbar at 2473 K when using the ultrasonic data of Rivers and Carmichael (1987) and their estimated value for (oKT/oT)0 (30 · 103 kbar/K). At the same temperature, the values of KT for our simulated MORB are in fair agreement with the above values since we obtain 175, 155 and 135 kbar, respectively. This is why our simulation data suggest for K0 a value around 7 instead of 4–5. For completeness we have also reported in Fig. 3 the compressibility curves predicted by the EOS of Ghiorso (2004) which is based on exhaustive data compilation. This EOS describes a MORB melt much less compressible above 100 kbar than what is observed by Ohtani and Maeda (2001) or evaluated by simulation. Agee and Walker (1988, 1993) have reported density measurements for molten komatiite compressed around 60 and 90 kbar at 2173 K. The data are shown in Fig. 4 with our simulation results for the same ultrabasic composition compressed isothermally at 2273 K. Since measurements were made at 2173 K (or rescaled for this temperature) we have also evaluated by simulation the density at this very temperature for a pressure of 85 kbar.

4544

B. Guillot, N. Sator / Geochimica et Cosmochimica Acta 71 (2007) 4538–4556

Fig. 3. Isothermal compression curves for simulated MORB at 2073 and 2473 K. The circles are the simulation data and the thin lines are the fits obtained with the BMEOS (the associated parameters are given in Table 3). The black square represents the density measurement of Agee (1998) at 1950 K and 58.5 kbar and the dashed line is the BMEOS proposed by this author (with q0 = 2.60 g/cm3, KT = 193 kbar and K0 = 4.4). The black triangle is the high pressure data of Ohtani and Maeda (2001) at 2473 K and 150 kbar and the dotted curve is the associated BMEOS (with q0 = 2.52 g/cm3, KT = 187 kbar and K0 = 5.0). The two bold lines are the isothermal compression curves (2073 and 2473 K) predicted by the EOS of Ghiorso (2004).

The resulting value almost coincides with the experimental data (see Fig. 4). Moreover, at zero pressure, the density of the simulated komatiite is only 1.2% smaller than that reported by Courtial et al. (1997) and the parameters of the BMEOS associated with the simulated points at 2173 K (namely KT = 171 kbar and K0 = 6.8) fit in the bracketing proposed by Courtial et al. (KT = 145  180 kbar and K0 = 8.2  5.6) in applying the BMEOS to high pressure data of Agee and Walker and their own one-atmosphere experimental results. So, one may conclude that the simulation reproduces satisfactorily the compression curve of komatiite melt from ambient to 90 kbar. As for the EOS of Ghiorso (2004) it tends to overestimate somewhat the compressibility of the melt above 60 kbar (see Fig. 4). For peridotite melt, the available data (Agee and Walker, 1993; Suzuki et al., 1995; Suzuki and Ohtani, 2003) cover a rather broad temperature–pressure range (from 80 kbar at 2273 K to 205 kbar at 2633 K). These data are shown in Fig. 5 with our own simulation results for the two isotherms 2273 and 2673 K. Although the overall agreement is satisfying, the simulation underestimates the experimental pressure by 10% and overestimates the density by 1.5% over the 80–205 kbar range. At zero pressure and 2273 K the density of the simulated melt amounts to 2.61 g/cm3 instead of 2.66 g/cm3 as extrapolated from the data of Courtial et al. (1997). Finally, here again the simulation is in better agreement with high pressure data than does the EOS of Ghiorso (see Fig. 5).

Fig. 4. Isothermal compression curve for molten komatiite. The circles are the simulation data at 2273 K and the thin curve is the associated BMEOS (for the parameters, see Table 3). The black squares are the data of Agee and Walker (1988, 1993) at 2173 K and the empty squares at P  0 and 85 kbar are our simulation data at the same temperature. The full dot is the density value at atmospheric pressure obtained from the data of Courtial et al. (1997) after a slight linear extrapolation to 2173 K. The dotted curve is the fit of our simulation data at 2173 K by the BMEOS. The bold line is the isothermal compression curve (2173 K) predicted by the EOS of Ghiorso (2004).

Fig. 5. Isothermal compression curves for molten peridotite at 2273 and 2673 K. The circles are the simulation data at 2273 and 2673 K and the thin curves are the associated BMEOS (for the parameters, see Table 3). The black squares are the data of Suzuki and Ohtani (2003) at (205 kbar, 2633 K) and (131 kbar, 2303 K), respectively, and the black triangles are the data points of Agee and Walker (1993) at (80 kbar, 2273 K) and (100 kbar, 2273 K). The full dot is the density value at atmospheric pressure obtained from the data of Courtial et al. (1997) after a linear extrapolation to 2273 K. The two bold lines are the isothermal compression curves predicted by the EOS of Ghiorso (2004).

Simulation study of natural silicate melts at high pressure

4545

3.3. Shock-wave equation of state Over the past 20 years the EOS of a number of silicate liquids have been evaluated by shock-wave experiments on targets preheated to liquidus temperature (Asimow and Ahrens, 2005). In these experiments, the initial sample density is inferred from measured temperature (by thermocouple or pyrometry) and known (or estimated) low density and thermal expansion data. Measurement of the shockwave velocity and the projectile velocity that are generated in the melt sample allows determination of the pressure– density relation along the principal Hugoniot by application of the Rankine–Hugoniot conservation equations. Thus the pressure range 50–450 kbar was achieved by this technique in MORB, komatiite, fayalite and the diopside– anorthite eutectic (model basalt). From the simulation point of view, the Hugoniot curve can be evaluated from the expression (for a review see Brennan and Rice, 2003), U H ¼ U 0 þ ðP H þ P 0 ÞðV 0  V H Þ=2

ð4Þ

which expresses the conservation of mass, momentum and energy during an adiabatic compression of the shocked sample. In Eq. (4), UH is the internal energy along the Hugoniot, U0 the internal energy of the sample in the initial state (P0, V0, T0) and PH, VH the pressure and molar volume in the Hugoniot state. The determination of the Hugoniot consists in finding the set of state points (PH, VH, TH) which satisfy the Eq. (4). In practice, for a given pressure P (>P0) one evaluates by simulation in the isothermal–isobaric ensemble the internal energy U as function of the temperature T (>T0). The intersection of the curve U(T) with the function, U0 + (P + P0)(V0  V)/2, then defines the temperature T at which the Eq. (4) is verified for a given pressure. Repeating this procedure for other P’s allows one to define the Hugoniot curve. (Each state point was evaluated with the following statistical uncertainties: ±0.02 g/cm3 for the density, ±4% for the temperature and ±5% for the pressure.) It should be noted that a shock compression although adiabatic, differs substantially from an isentropic compression. The latter is reversible while a shock compression is not. Hence for high shock waves the pressure along the isentrope is below that on the Hugoniot whereas for weak shock waves the Hugoniot is very nearly isentropic and the pressures are almost identical to each other. So, for the sake of completeness we have also evaluated the isentrope by performing a series of calculations in the microcanonical ensemble (NVE). In starting from an initial state (N0, V0, E0) corresponding to a temperature T0 and a pressure P0, the progressive reduction of the volume (V < V0) at constant energy E0 (adiabatic compression) leads to the corresponding isentropic compression curve, (P0, V0, T0)E0 fi (P, V, T)E0. Fig. 6 presents the Hugoniot data of Chen et al. (2002) for molten fayalite initially at 1573 K and 1 bar, as well as our simulation results obtained with the same initial conditions. The agreement between experiment and simulation is excellent up to 250 kbar. Above 250 kbar, the experimental Hugoniot becomes erratic (see also Fig. 10 in Chen

Fig. 6. Pressure–density Hugoniot for molten fayalite initially at 1573 K. The circles are the simulation results for the Hugoniot (the dotted curve is a guideline for the eyes) and the full dots represent the Hugoniot data of Chen et al. (2002). The triangles are the simulation results for the isentrope (the dashed curve is a guideline for points above 250 kbar) and the thin curve is the BMEOS associated with this isentrope in the 0–250 kbar range (for the parameters, see Table 5). Notice that the BMEOS becomes unreliable above 230 kbar. The inset shows the temperature– pressure relationship along the simulated Hugoniot (circles) and the simulated isentrope (triangles), the dotted and dashed curves being just guidelines for the eyes.

et al., 2002) and it is difficult to make a constructive comparison. Interestingly enough the simulated isentrope and the simulated Hugoniot are very close to each other up to 300 kbar, the Hugoniot separates rapidly from the isentrope at higher pressure. Correspondingly, the temperature along the Hugoniot increases steadily with the pressure (see the inset in Fig. 6) when it tends to level off along the isentrope, as expected theoretically (Rigden et al., 1988). It is useful to fit our results for the isentrope by a BMEOS. Values of KS and K 0S are 220 kbar and 7.8, respectively, as compared with 259 kbar and 5.36 deduced by Chen et al. (see also Table 5). It is meaningful that our value of KS at 1573 K and zero pressure is equal to the ultrasonic measurement of Rivers and Carmichael (1987). A similar investigation was carried out on an initially molten komatiite at 1823 K. Results of the simulations and Hugoniot data (Miller et al., 1991) on a synthetic komatiite are shown in Fig. 7. The calculated Hugoniot is very close to the experimental data up to approximately 100 kbar but deviates upwards at higher compression rate. The simulated isentrope is virtually undistinguishable from the simulated Hugoniot up to 120–130 kbar and it deviates significantly from the isentrope deduced by Miller et al. (1991) only above 100 kbar (values of KS and K 0S are collected in Table 5). As far as the experimental data are reliable above 150 kbar, one concludes that the simulated komatiitic melt is slightly less compressible than the experimentally investigated komatiite at high compression rate. This difference may be ascribed to a small variation in com-

4546

B. Guillot, N. Sator / Geochimica et Cosmochimica Acta 71 (2007) 4538–4556

Table 5 Birch–Murnaghan equation parameters for isentropic compression of various silicate melts T0 (K)

q0 (g/cm3)

KS (kbar)

K 0S

Ref.

An0.36Di0.64

1673 1673 1673

2.70 2.615 2.61

235 226 242

7.7 4.15 4.85

This study Rigden et al. (1984) Rigden et al. (1988)

MORB

1673 1673 1673

2.69 2.66 2.719

225 112 244

7.0 5.9 4.87

This study Rowan et al. (1991) Chen et al. (2002)

BMars

1673

2.76

200

7.6

This study

Komatiite

1823 1823

2.73 2.745

210 270

8.2 4.9

This study Miller et al. (1991)

Peridotite

1873

2.75

260

7.0

This study

Olivine

1873

2.75

270

7.2

This study

Fayalite

1573 1573

3.75 3.75

220 259

7.8 5.36

This study Chen et al. (2002)

Fig. 7. As in Fig. 6 but for molten komatiite initially at 1823 K. The circles are our simulation results for the Hugoniot (the dotted curve is just a guideline) and the triangles are those for the isentrope. The full dots are the Hugoniot data of Miller et al. (1991). The thin curve is the BMEOS associated with the simulated isentrope (for the parameters, see Table 5), and the dash-dotted line is the BMEOS proposed by Miller et al. for the isentrope. Notice that the BMEOS associated with the simulated isentrope becomes unreliable above 350 kbar (compare the thin curve with the dashed curve which is a guideline for the simulation points above 350 kbar). The inset shows the temperature–pressure relationship along the simulated Hugoniot (circles) and the simulated isentrope (triangles), the dotted and dashed curves being just guidelines for the eyes.

position between the two komatiites; the one investigated by Miller et al. (1991) being slightly richer in Al2O3 and poorer in MgO and FeO than the simulated one. Concerning molten basalt, Rigden et al. (1984, 1988) have investigated the anorthite–diopside eutectic to model a MORB. The Hugoniot data for this Fe-free model basalt are shown in Fig. 8 with our simulation results for a true MORB composition (Kane zone of central Atlantic) and

Fig. 8. As in Fig. 6 but for molten MORB initially at 1673 K. The circles are our simulation results for the Hugoniot of MORB (the dotted curve is just a guideline), the triangles are those for the isentrope and the thin curve is the BMEOS associated with this isentrope (for parameters, see Table 5). Notice that the latter BMEOS becomes unreliable above 270 kbar (compare the thin curve with the dashed curve which is a guideline for the simulation points above 270 kbar). The empty squares are our simulation results for the isentrope of An0.36Di0.64 and the thin curve is the BMEOS associated with this isentrope (for parameters, see Table 5). Notice that the latter BMEOS becomes unreliable above 270 kbar (compare the thin curve with the dashed curve which is a guideline for the simulation points). The full dots are the Hugoniot data of Rigden et al. (1988) for the anorthite–diopside eutectic. The inset shows for MORB, the temperature–pressure relationship along the simulated Hugoniot (circles) and the simulated isentrope (triangles), the dotted and dashed curves being just guidelines for the eyes.

the same initial conditions (1673 K, P  0). The agreement between the two Hugoniots is satisfactory if one considers the slight shift in density at 1 bar between simulated and real melts (2.70 g/cm3 for simulated An0.36Di0.64 and 2.69 g/cm3 for simulated MORB instead of 2.63 and

Simulation study of natural silicate melts at high pressure

2.67 g/cm3 as observed experimentally, see Table 4). To check the sensitivity of the simulation to a change in composition, we have calculated the isentrope for the anorthite– diopside eutectic and found that it matches fairly well the Hugoniot data of Rigden et al. (1988) in spite of the slight density shift (3%) at low pressure already mentioned (see Fig. 8). Hence, we are confident in our prediction for real MORB. To be complete, it is noteworthy that Rowan et al. (1991) have reported some shock-wave data on a real MORB (see Table 5), however it has been recognized that their measurements are in strong disagreement with static compression data (Agee, 1998) and with the shock-wave compression study of Rigden et al. (1988) presently discussed. Through another way, in an attempt to better predict the compressibility of an idealized molten basalt, Chen et al. (2002) have estimated the EOS of an hypothetical mixture (An0.36Di0.64)Fa0.15 whose composition is close to that of an average MORB. However, the resulting BMEOS KS = 244 kbar and parameters (q0 = 2.719 g/cm3, 0 K S ¼ 4:87) lead to an isentrope which describes a melt much more compressible above 70 kbar than the one obtained in the present study (q0 = 2.69 g/cm3, KS = 225 kbar and K 0S ¼ 7:0). Finally, for the sake of heuristicity, we have also evaluated the isentropes of a Fe-rich Martian basalt, an olivine melt and a peridotite melt. The isentrope of the Martian basalt shows a compressibility curve very close to that of a MORB except above 200 kbar where its compressibility is slightly larger (not shown). As for olivine and peridotite, their isentropic compressibility is found to be close to each other (and to that of komatiite as well) with the sequence, Ol < Pe < Ko, the olivine melt being the less compressible one (all values of KS and K 0S are given in Table 5). This is in accordance with the trend exhibited by the adiabatic bulk modulus with MgO content in magnesiosilicates (see Secco et al., 1991b). 4. STRUCTURAL CHANGES IN HIGH PRESSURE LIQUIDS In the last two decades a number of experimental studies using mostly NMR spectroscopy (but also Raman, infrared and XANES) have been devoted to the investigation of glasses quenched from high pressure melts (Stebbins and McMillan, 1989; Stebbins and Sykes, 1990; Xue et al., 1991, 1994; Yarger et al., 1995; Poe et al., 2001; Lee et al., 2003, 2004; Lee, 2004, 2005; Allwardt et al., 2004, 2005a; Lee et al., 2006). Alkali silicates (e.g. Na2Si2O5, Na2Si3O7, Na2Si4O9 and K2Si4O9) and aluminosilicate glasses (e.g. NaAlSi3O8) have been investigated as analogs of natural melts (granitic to basaltic composition). An important finding of these studies is that the pressure-induced structural changes taking place in silicate melts and glasses are complex function of composition. For instance, the abundance of highly coordinated species [5,6]Si and [5,6]Al increases all the more rapidly with increasing pressure such that the melt is more depolymerised at atmospheric pressure. Moreover, the polymerization of the melt tends to increase with pressure as a result of the conversion of a number of nonbridging oxygens (NBO) into bridging oxy-

4547

gens (BO), the magnitude of this conversion depends on the degree of polymerization of the melt at low pressure (the lower this degree the greater the conversion). Recent progress on determining the structure of melts have been achieved thanks to 2D triple quantum magic-angle spinning (3QMAS) NMR. It is now possible to identify oxygens linking 4-coordinated Si and Al to 5- and 6-coordinated Al and Si (e.g. [4]Si–O–[5,6]Al and [4]Al–O–[5,6]Si) as well as the oxygens linking a network modifier cation to a highly coordinated network former cation (e.g. Na–O–[5,6]Si). Furthermore, the chemical shift of the BO and NBO peaks in NMR spectra of glasses quenched from high pressure melts indicates that the T–O bond lengths (with T = Si, Al) seems to increase under pressure. In this context computer simulations are of particular interest as they lead directly to the atomic configurations. We have shown previously (Guillot and Sator, 2007) that the main structural parameters, average cation–oxygen distances and cation coordination numbers, as deduced from diffraction studies (Brown et al., 1995) on various silicate melts around atmospheric pressure, are well reproduced by our simulations. So, for illustration, we have evaluated by MD simulations the evolution of the structure of the melt with the pressure for three compositions, namely rhyolite, MORB and peridotite. 4.1. Pressure-induced coordination changes in magmatic liquids In Fig. 9 is presented a population analysis of n-coordinated Al and Si network former cations in the three melts at the occurrence 2273 K and P = 0, 50, 100, 150 and 200 kbar. In our calculation a cation is considered as ncoordinated if there are n oxygen atoms within a sphere of radius rcut, where rcut is the location of the first minimum in the corresponding pair correlation function (in practice ˚ and rcut(Al–O) = 2.60 A ˚ ). In accorrcut(Si–O) = 2.30 A dance with the findings of the experimental literature discussed above, we found that the fraction of 5- and 6coordinated Al and Si increases continuously with the pressure and that the magnitude of this increase depends on the degree of polymerization of the melt at low pressure (more depolymerized the melt higher the population increase of [5,6] Al and [5,6]Si at isobaric conditions). Correlatively the mean T–O bond length increases with coordination number ˚ for [4]Si at low pressure (1.83 A ˚ for [4]Al) to from 1.67 A [5] [5] ˚ ˚ ˚ for [6]Si 1.77 A for Si (1.89 A for Al) and 1.84 A ˚ for [6]Al) at 200 kbar, values which depend very lit(1.97 A tle on composition. Nevertheless a significant population of 6-coordinated Si is delayed at a higher pressure with respect to [6]Al. For instance in rhyolite at 50 kbar only 0.3% of [6]Si is found as compared with 7% of [6]Al. Interestingly, the fraction of [5]Al goes through a maximum value when the pressure increases (around 150 kbar in rhyolite and 100 kbar in MORB and in peridotite), a feature observed experimentally by Yarger et al. (1995) in Na3AlSi7O17. Moreover, the highly coordinated species [5,6]Al and [5,6]Si are generally found in a higher proportion in our simulations than it is observed by NMR in aluminosilicate glasses. Thus in a fully polymerized glass (NaAlSi3O8) quenched

4548

B. Guillot, N. Sator / Geochimica et Cosmochimica Acta 71 (2007) 4538–4556

Fig. 9. Population (in %) of n-coordinated Si and Al in simulated melts as function of pressure. The sticks are displaced on the right side when the pressure increases by increment of 50 kbar (P = 0, 50, 100, 150 and 200 kbar) at fixed temperature (2273 K).

from a high pressure melt, Lee et al. (2004) estimated the fraction of [5]Al to about 7–8% at 80 kbar when Allwardt et al. (2005b) estimated this fraction to about 15% in the same glass at 100 kbar, as compared with our evaluation 47% in rhyolite at 100 kbar. For basaltic like compositions, Lee et al. (2004) estimated the total fraction [5] Al + [6]Al about 37% in (Na2O)0.75(Al2O3)0.253SiO2 at 80 kbar when Allwardt et al. (2005b) obtained 62% (with 39% of [5]Al and 23% of [6]Al) in Na3AlSi7O17 at 100 kbar and we obtain 49% [5]Al + 21.5% [6]Al in MORB at 100 kbar. As frequently stressed in the literature (e.g. Allwardt et al., 2004, 2005b), it is likely that the evaluation of highly coordinated species in quenched glasses underestimate the true values reached in the liquid phase. 4.2. Compression mechanisms in melts The difference between Al and Si to a pressure increase is well illustrated in Fig. 10, where the average coordination numbers of Al and Si are reported as function of pressure for the three investigated melts (the coordination numbers are calculated by integrating the first peak of the corresponding pair correlation function). The average coordination number of Al increases much more rapidly with pressure than the coordination number of Si does. Furthermore, the average coordination number of Si seems to increase significantly only above 80 kbar (Nc = 4.0 at P = 0

Fig. 10. Pressure dependence of the average coordination numbers for Si, Al, and Mg in rhyolite (full circles), MORB (empty squares) and peridotite (full triangles) at 2273 K.

and 4.1 at 75 kbar), which suggests that below this pressure another compression mechanism is responsible for the densification. This is apparent in the bond angle distributions

Simulation study of natural silicate melts at high pressure

T–O–T and T–O–T0 (with T, T0 = Si, Al) presented in Fig. 11. At low pressure the Si–O–Si intertetrahedral angle exhibits a broad distribution centered about 142 in rhyolite (as compared with 144 in silica), 133 in MORB and 135 in peridotite melt. Upon initial compression (e.g. 50 kbar), the maximum of the distribution shifts towards a smaller angle, the shift being less pronounced when the melt is more depolymerized. This is in agreement with the pressure-induced folding mechanism of corner-shared silicate tetrahedra proposed by Stolper and Ahrens (1987). Above 100 kbar a marked shouldering appears in the bond angle distribution around 95 which corresponds to a rising population of [5]Si and [6]Si species. In contrast, the Al–O–Al bond angle distribution at low pressure is broad and rather flat between 90 and 120 meaning that adjacent AlO4 tetrahedra prefer to be in a folded configuration. When the pressure increases, a genuine peak develops rapidly around 90 which is the signature of the presence of [5]Al and [6]Al species. As for the Si–O–Al bond angle distribution, it presents some similarities with the Si–O–Si and Al–O–Al distributions. Thus upon initial compression the angle between adjacent SiO4 and AlO4 tetrahedra narrows rapidly while the coordination of Al cations increases by attracting nonbonding oxygens located in the neighbourhood. The peridotite case is particular because the structure of the melt is a highly distorted network composed essentially

4549

of interconnected SiO4 tetrahedra (Al plays a minor role due to its low concentration) and MgOn = 4, 5, 6 polyhedra (Kohara et al., 2004). This is illustrated in Fig. 11 where the Si–O–Mg bond angles present a broad distribution at low pressure. With increasing pressure, this distribution narrows drastically and peaks at about 90 as a consequence of the progressive folding of the intertetrahedral angle between adjacent SiO4 units and between a SiO4 unit and its connected MgOn polyhedron. In contrast the Mg– O–Mg bond angle distribution evolves only marginally with pressure, because the interconnected MgOn polyhedra are in a compact configuration from the low pressure. In fact, its evolution is essentially governed by the increase of the coordination number of Mg with pressure (from 4.9 around zero pressure to 6.8 at 200 kbar, see Fig. 10). As for the Si– O–Si bond angle distribution, its evolution with pressure is similar to that exhibited by MORB (see Fig. 11) which expresses the folding of the SiO4 tetrahedra upon initial compression followed by the occurring of [5]Si and [6]Si species at higher pressure (see the prominent feature near 90 in Fig. 11). 4.3. Degree of polymerization in high pressure melts The above findings are the expression of the pressure-induced redistribution of oxygens through the melt. This

Fig. 11. Bond angle distribution (T–O–T with T = Si, Al or Mg) in rhyolite, MORB and peridotite for P = 0, 50, 100, 150 and 200 kbar. The arrows indicate the evolution of the distribution when the pressure increases.

4550

B. Guillot, N. Sator / Geochimica et Cosmochimica Acta 71 (2007) 4538–4556

redistribution depends on the degree of polymerization of the melt at ambient pressure. As shown in Fig. 12, where a population analysis of oxygen atoms is presented, this redistribution is much more important in a highly polymerized melt as rhyolite than in a depolymerized melt as peridotite. For the highly compressible rhyolite, the distribution of the oxygens is asymmetric and shifts rapidly (upon initial compression) towards larger numbers of neighbours with an average coordination number which increases from 7.0 at low pressure to 12.5 at 200 kbar. For MORB the distribution is more symmetrical and the average coordination number increases from 9.3 to 12.2 at 200 kbar. For the weakly compressible peridotite melt, the distribution is symmetrical and shifts very little with the pressure: the coordination number only varies from 11.6 to 12.0 between 0 and 200 kbar. Correspondingly, for all composition, the average O–O distance decreases from ˚ at low pressure to 2.57 A ˚ at 200 kbar. 2.66 A Key information is given by the evolution with pressure and composition of the population of bridging oxygens (BO) and nonbridging oxygens (NBO). A detailed population analysis of BOs and NBOs is given in Fig. 13 as function of pressure and composition. The BO atoms are

Fig. 12. Population analysis of oxygen atoms in the first shell around a central oxygen as function of the pressure. In our calculations the extension of the first shell is defined by the location of the first minimum in the oxygen–oxygen pair distribution function. In practice, when the pressure is increasing from 0 to ˚ for 200 kbar, this distance evolves in the range 3.50–3.60 A ˚ for MORB and 4.09–3.60 A ˚ for peridotite. rhyolite, 3.90–3.70 A The dotted curves are guidelines to better see the envelope of the populations for a given pressure (P = 0, 50, 100, 150 and 200 kbar). Notice that the population distribution shifts steadily to the right when pressure increases.

defined as oxygen atoms that are connected to two network former cations (Si or Al) whereas NBO atoms are oxygen atoms connected to only one network former cation. In our calculations two atoms (e.g. Si–O) are connected if the distance between the two atoms is less than the location of the first minimum in the corresponding pair distribution ˚ for Si– function (in practice these cutoff distances are 2.30 A ˚ O and 2.60 A for Al–O). For the BOs we have distinguished the Si–O–Si, Si–O–Al and Al–O–Al bonds and also the triclusters O–(Si)3, Al–O–(Si)2, (Al)2–O–Si and O–(Al)3 where an oxygen is threefold coordinated to network former cations. For the NBOs we have distinguished the oxygens linked to one network former cation (Si or Al) from the free oxygens only linked to network modifier cations. For rhyolite one notices that the total number of BOs (the sum of all contributions) is high, as expected, and barely increases with the pressure (93% at zero pressure and 96% at 200 kbar). However the respective populations of BOs are modified by the pressure. Thus the number of triclusters (especially Al–O–(Si)2 and O–(Si)3) increases progressively from 5% at low pressure to 33% at 200 kbar. This rising population of triclusters is formed at the expense of Si–O–Si and Al–O–Si clusters and is favoured by the emergence of highly coordinated [5,6]Si and [5,6]Al species. Correlatively, the pressure induces the transformation of [4] Si–O–[4]Si and [4]Al–O–[4]Si into [4,5]Si–O–[4]Si and [5,6] Al–O–[4,5]Si. In MORB, the total number of BOs increases from 65% at zero pressure to 79% at 200 kbar. Most of this increase is due, here also, to a rising population of triclusters (from 3.2% to 18.8% at 200 kbar). For the sake of comparison with experimental data, notice that for Na2O–3SiO2, a partially depolymerised melt used as a model for basaltic melts, Lee et al. (2003) measured by NMR a decrease of NBO fraction from 28% at 1 bar to 22% at 100 kbar as compared with a decrease from 35% to 27.4% in our MORB at the same pressures. An excellent agreement if one considers the large difference in composition between the two melts. In the peridotitic melt, the total number of NBOs decreases from 72% at zero pressure to 55% at 200 kbar (increasing polymerization). More precisely, as indicated by Fig. 13, when the pressure increases the abundances of Mg–O–(Si, Al) bonds and the Mg–O–Mg bonds diminish to the benefit of the Si–O–(Si, Al) bonds. This means that a fraction of the oxygens connecting SiO4–MgOn and MgOn–MgOn polyhedra in the low pressure melt tends to migrate preferentially towards SiO4 units to form an interconnected network of SiOn = 4, 5, 6 polyhedra at high pressure. Here the role of triclusters is marginal even if their abundance increases with pressure (from 0.15% at zero pressure to 3.5% at 200 kbar). To be complete, it is noteworthy that the abundance of free oxygens is negligible in rhyolite, very low in MORB (1%) and significant in peridotite (10%) where they correspond essentially to Mg–O–Mg bonds. 5. DIFFUSIVITY AND ELECTRICAL CONDUCTIVITY In order to get some insight about the link between viscous flow, transport properties and melt structure, many

Simulation study of natural silicate melts at high pressure

4551

Fig. 13. Population of bridging oxygens (BO) and nonbridging oxygens (NBO) as function of the pressure. The symbols are the following. For BOs: Si–O–Si (filled circles), Al–O–Si (empty circles), Al–O–Al (crosses), Al–O–(Si)2 (full squares), (Al)2–O–Si (empty squares), O–(Si)3 (full triangles), O–(Al)3 (empty triangles). For NBOs: O–Si (full pentagons), O–Al (empty pentagons), O linked only to network modifying cations (stars). Notice that in peridotite the population of free oxygens (stars) is represented essentially by Mg–O–Mg bonds.

experimental studies have been devoted to the diffusion of Si and O atoms into silicate melts. Thus in polymerized melts like jadeite (Shimizu and Kushiro, 1984), aluminosilicate liquids (Poe et al., 1997) and Na2Si4O9 (Rubie et al., 1993), it was found that both O and Si self diffusion coefficients increase on initial compression and reach a maximum value before decreasing at higher pressure. In depolymerized melts like diopside, Reid et al. (2001) have observed in the liquid at 2273 K a decrease of the Si and O self diffusions with pressure up to 110 kbar after which an increase is observed up to 150 kbar, the highest pressure investigated. For compositions close to those of natural silicates, Tinker and Lesher (2001) observed that Si and O exhibit a diffusivity maximum in dacitic liquid at about 40 kbar whereas for a basaltic composition (the diopside–anorthite eutectic), Tinker et al. (2003) noticed a continuous increase of the diffusivities in the pressure range 10–40 kbar. We have evaluated by MD simulations the self diffusion coefficients of all ionic species in rhyolite, andesite, MORB, peridotite and komatiite melt at 2273 K in the pressure range 0–200 kbar. The results are reported in Fig. 14 for network formers and in Fig. 15 for structure modifiers.

The self diffusion coefficient was computed from the MD run via the Einstein relation (Kubo, 1966) Ds ¼ lim

t!1

Ns 1 X hð~ ri ðtÞ ~ ri ð0ÞÞ2 i N s i¼1 6t

ð5Þ

where ri is the position of the ith ion of species s, and where the bracket expresses an ensemble average taken over many origin times. Long simulation runs were performed (10 ns long and 103 origin times) to reduce the uncertainties on D value to around ±5% for the most abundant ions (Si and O) and around ±40% for the less abundant cations (Fe3+, Na and K). A quick look at Figs. 14 and 15 shows that network formers and network modifiers behave differently with pressure. When the diffusivity of network modifiers is always decreasing when the pressure raises, the diffusivity of network formers exhibits a more complex behavior with pressure. The diffusivity of Si, O, Al and Fe3+ in rhyolite exhibits a maximum at about 100–150 kbar, at about 50 kbar in MORB, and at about 25 kbar in molten peridotite. In andesite two maxima are visible (at 25 and 100 kbar) whereas in molten komatiite no maximum

4552

B. Guillot, N. Sator / Geochimica et Cosmochimica Acta 71 (2007) 4538–4556

Fig. 14. Pressure dependence of the self diffusion coefficients of network former ions in rhyolite (full dots), andesite (empty circles), MORB (squares), komatiite (empty triangles) and peridotite (full triangles). The temperature of the simulated melts is 2273 K.

appears. Notice that the values for Al and Fe3+ are calculated with a greater uncertainty than for Si and O, a feature which tends to blur somewhat the similarity of behavior between all these ions for a given composition. The maximum of diffusivity for Si and O observed in dacitic liquid around 40–50 kbar (Tinker and Lesher, 2001) is quite compatible with our results for acidic melts if one considers that the dacite composition is intermediate between that of rhyolite and andesite. In the same way, the weak increase of Si and O self diffusion coefficients observed in diopside– anorthite melt between 10 and 40 kbar (Tinker et al., 2003) parallels the barely visible variation of the corresponding diffusivities calculated for MORB in the same pressure range. For molten peridotite, the simulation data show a continuous decrease of the diffusivities with increasing pressure (when P > 25 kbar) as it is observed by Reid et al. (2001) in diopside liquid (considered as a peridotitic melt analog) between 30 and 110 kbar. Nevertheless the diffusivity increase observed between 110 and 150 kbar is not predicted for MORB, a discrepancy may be due to differences in composition between the two melts. The pressure enhancement of ion mobilities in highly polymerized liquid silicates is well established by simulation studies (Angell et al., 1982; Nevins and Spera, 1998; Bryce et al., 1999). For strong glass formers like silica it has been argued that upon compression those melts undergo a strong-to-fragile transition (Angell, 1995; Guillot and Guissani, 1997; Horbach and Kob, 1999; Saika-Voivod et al., 2004). In a polymerized melt, upon initial compression the folding of the intertetrahedral angles

(see T–O–T angles in Fig. 11) modifies the potential energy landscape. The barrier height to diffusion decreases between oxygens belonging to adjacent tetrahedral and the diffusivity of oxygen increases. However this process becomes less effective upon further compression because the barrier height to diffusion rises up due to the packing. Then a normal behavior is restored, the diffusivity decreases upon further compression. This scenario likely applies to acidic melts like rhyolite, dacite and andesite whereas for depolymerized melts the role of nonbridging oxygen is predominant and a normal behavior of the diffusivity is expected with pressure (for a semi-quantitative theory on anomalous pressure dependence O2 diffusivity, see Lee et al., 2004). In this context one can wonder whether the diffusivity anomaly discussed above can be observed by investigating the electrical conductivity of molten rocks. So we have evaluated the electrical conductivity of our silicate liquids at 2273 K, a temperature at which the viscosity of all simulated compositions is lower than 1.PaS (according to Eyring equation). This latter assumption is likely fulfilled by real melts except rhyolite for which the viscosity is higher, about 100 PaS at ambient pressure (Neuville et al., 1993). Keeping in mind this information, the electrical conductivity was evaluated from the mean square displacement of the ions (Kubo, 1966), D P 2 E  N zi ð~  r ðtÞ ~ r ð0ÞÞ i i i¼1 1 r ¼ lim t!1 kTV 6t

ð6Þ

Simulation study of natural silicate melts at high pressure

4553

Fig. 15. As in Fig. 14 but for network modifying cations.

where k is the Boltzmann constant, V the volume of the sample, N the total number of ions, and zi the charge carried by the ion i. In Fig. 16 the evolution of the conductivity with pressure for four different compositions is shown. Despite the large error bars associated with the computed values, it appears that the pressure dependence of the conductivities of rhyolite, andesite and MORB shares some similarity with the pressure behavior exhibited by the diffusion coefficients of network former ions in these melts whereas that of peridotite is also reminiscent of the diffusivity of Mg2+ (compare Fig. 16 with Figs. 14 and 15). As a matter of fact, for rhyolite, andesite and MORB the conductivity increases upon initial compression, reaches an inconspicuous maximum about or above 100 kbar and tends to saturate at higher pressure whereas for peridotite the conductivity decreases upon initial compression. As far as the magnitude of the conductivity is concerned, the calculated values may seem very high (10–150 S/m). In fact a simple high temperature extrapolation up to superliquidus temperatures of existing conductivity data for various silicate melts at low pressure (Presnall et al., 1972; Rai and Manghnani, 1978; Tyburczy and Waff, 1983; Sakamoto et al., 2002; Maumus et al., 2005) suggests the same order of magnitude for the conductivity of real liquids. 6. CONCLUSION From this simulation study of natural silicate melts under pressure, we will emphasize the following main points. The equations of state of our simulated liquids are in very

good agreement with the available high pressure data on MORB, komatiite, peridotite and fayalite as obtained either by sink/float experiment or by shock-wave compression. Moreover, the ability of the simulation to reproduce accurately the bulk moduli of a large variety of silicate liquids as evaluated by ultrasonic studies has to be emphasized. Our calculations also suggest that the temperature dependence of the bulk modulus cannot be neglected when investigating the compressibility curve of a melt at superliquidus temperature. Besides, the good reproduction of the available Hugoniot data is an important test in this respect. A consequence of taking into account the decrease of the bulk modulus with increasing temperature is that a value of K0 around 7–8 is required to reproduce at best by a BMEOS, the high temperature–high pressure data of silicate liquids. As far as the structure of the liquid silicates under pressure is concerned, we have stressed that the average coordination number of Al increases rapidly upon initial compression whereas for Si it varies very little in the 0–50 kbar range and more noticeably only above 80 kbar. At low and moderate pressure, this is the reduction of the intertetrahedral angle (T–O–T) which is the main structural response of the melt to the pressure, while at higher pressure it is the occurrence of penta- and hexacoordinated (Al, Si) species. Our results show that the degree of polymerization of the melt increases when applying pressure and that this tendency is all the more marked that the melt is more depolymerized at low pressure. These findings are in accordance with current ideas based on NMR data. Furthermore, in high pressure melts

4554

B. Guillot, N. Sator / Geochimica et Cosmochimica Acta 71 (2007) 4538–4556

Fig. 16. Pressure dependence of the electrical conductivity in silicate melts simulated at 2273 K. Notice the large error bars associated with the calculated points.

an important contribution of penta- and hexa-coordinated network former cations is to promote a significant fraction of triclusters. Although the presence of triclusters is frequently mentioned in the literature (Toplis et al., 1997; Poe et al., 2001; Tossell and Cohen, 2001; Stebbins et al., 2001; Lee et al., 2004; Allwardt et al., 2005b; Benoit et al., 2005), the present study suggests that their role could be quite significant in granitic and basaltic melts under pressure. Finally, it was shown that the diffusivies of network former ions exhibit a complex behaviour with pressure which depends on the degree of polymerization of the melt (the behaviour is nonmonotonous in acidic and basaltic liquids and tends to be simply decreasing in ultrabasic liquids). In contrast, the diffusivity of network modifiers is always decreasing when the pressure rises. Although the calculated electrical conductivity and its pressure dependence seem to share some similarities with the diffusivity behaviour, the relationship is not straightforward. ACKNOWLEDGMENT We thank two anonymous reviewers for critical comments and constructive suggestions. REFERENCES Agee C. B. (1998) Crystal-liquid density inversions in terrestrial and lunar magmas. Phys. Earth Planet. Inter. 107, 63–74.

Agee C. B. and Walker D. (1988) Static compression and olivine flotation of ultrabasic silicate liquid. J. Geophys. Res. 93, 3437–3449. Agee C. B. and Walker D. (1993) Olivine flotation in mantle melt. Earth Planet. Sci. Lett. 114, 315–324. Allen M. P. and Tildesley D. J. (1987) Computer Simulations of Liquids. Oxford University Press. Allwardt J. R., Schmidt B. C. and Stebbins J. F. (2004) Structural mechanisms of compression and decompression in highpressure K2Si4O9 glasses: an investigation utilizing Raman and NMR spectroscopy of glasses and crystalline materials. Chem. Geol. 213, 137–151. Allwardt J. R., Stebbins J. F., Schmidt B. C., Frost D. J., Withers A. C. and Hirschmann M. M. (2005a) Aluminum coordination and the densification of high pressure aluminosilicate glasses. Am. Mineral. 90, 1218–1222. Allwardt J. R., Poe B. T. and Stebbins J. F. (2005b) The effect of fictive temperature on Al coordination in high-pressure (10 GPa) sodium aluminosilicate glasses. Am. Mineral. 90, 1453–1457. Angell C. A. (1995) Simulation of glasses and glass-forming liquids after two decades: some perspectives. Comp. Mat. Sci. 4, 285–291. Angell C. A., Cheeseman P. A. and Tamaddon S. (1982) Pressure enhancement of ion mobilities in liquid silicates from computer simulation studies to 800 kilobars. Science 218, 885–887. Asimow P. D. and Ahrens T. J. (2005) Silicate liquid equations of state from molten shock experiments. AGU Abstract #MR12A-03. Benoit M., Profeta M., Mauri F., Pickard C. J. and Tuckerman M. E. (2005) First-principles calculation of the 17O NMR parameters of a calcium aluminosilicate glass. J. Phys. Chem. B 109, 6052–6060.

Simulation study of natural silicate melts at high pressure Birch F. (1978) Finite strain isotherm and velocities for singlecrystal and polycrystalline NaCl at high pressures and 300 K. J. Geophys. Res. 95, 1257–1268. Brennan J. K. and Rice B. M. (2003) Efficient determination of Hugoniot states using classical molecular simulation techniques. Mol. Phys. 101, 3309–3322. Brown, Jr., G. E., Farges F. and Calas G. (1995) X-ray scattering and X-ray spectroscopy studies of silicate melts. Rev. Mineral. 32, 317–410. Bryce J. G., Spera F. J. and Stein D. J. (1999) Pressure dependence of self-diffusion in the NaAlO2–SiO2 system: compositional effects and mechanisms. Am. Mineral. 84, 345–356. Caprion D. and Schober H. R. (2002) Influence of the quench rate and the pressure on the glass transition temperature in selenium. J. Chem. Phys. 117, 2814–2818. Chen G. Q., Ahrens T. J. and Stolper E. M. (2002) Shock-wave equation of state of molten and solid fayalite. Phys. Earth Planet. Inter. 134, 35–52. Circone S. and Agee C. B. (1996) Compressibility of molten high-Ti mare glass: evidence for crystal-liquid density inversions in the lunar mantle. Geochim. Cosmochim. Acta 60, 2709–2720. Courtial Ph., Ohtani E. and Dingwell D. R. (1997) High-temperature densities of some mantle melts. Geochim. Cosmochim. Acta 61, 3111–3119. Fujii T. and Kushiro I. (1977) Density, viscosity and compressibility of basaltic liquid at high pressures. Carnegie Inst. Wash Year book 76, 419–424. Ghiorso M. (2004) An equation of state for silicate melts. IV. Calibration of a multicomponent mixing model to 40 GPa. Am. J. Sci. 304, 811–838. Ghiorso M. and Kress V. C. (2004) An equation of state for silicate melts. II. Calibration of volumetric properties at 105 Pa. Am. J. Sci. 304, 679–751. Guillot B. and Guissani Y. (1997) Boson peak and high frequency modes in amorphous silica. Phys. Rev. Lett. 78, 2401–2404. Guillot B. and Guissani Y. (2003) Polyamorphism in low temperature water: a simulation study. J. Chem. Phys. 119, 11740–11752. Guillot B. and Sator N. (2007) A computer simulation study of natural silicate melts. Part I: low pressure properties. Geochim. Cosmochim. Acta 71, 1249–1265. Horbach J. and Kob W. (1999) Static and dynamics properties of a viscous silica melt. Phys. Rev. B 60, 3169–3181. Knoche R., Dingwell D. B. and Webb S. L. (1992) Temperaturedependent thermal expansivities of silicate melts: the system anorthite–diopside. Geochim. Cosmochim. Acta 56, 689–699. Kohara S., Suzuya K., Takeuchi K., Loong C.-K., Grimsditch M., Weber J. K. R., Tangeman J. A. and Keys T. S. (2004) Glass formation at the limit of insufficient network formers. Science 303, 1649–1652. Kress V. C. and Carmichael I. S. E. (1991) The compressibility of silicate liquids containing Fe2O3 and the effect of composition, temperature, oxygen fugacity and pressure on their redox states. Contrib. Mineral. Petrol. 108, 82–92. Kubo R. (1966) The fluctuation-dissipation theorem. Rep. Prog. Phys. 29, 255–284. Lange R. A. and Carmichael I. S. E. (1987) Densities of Na2O– K2O–CaO–MgO–FeO–Fe2O3–Al2O3–TiO2–SiO2 liquids: new measurements and derived partial molar properties. Geochim. Cosmochim. Acta 51, 2931–2946. Lee S. K., Fei Y., Cody G. D. and Mysen B. O. (2003) Order and disorder in sodium silicate glasses and melts at 10 GPa. Geophys. Res. Lett. 30, 1845. Lee S. K. (2004) Structure of silicate glasses and melts at high pressure: quantum chemical calculations and solid-state NMR. J. Phys. Chem. B 108, 5889–5900.

4555

Lee S. K. (2005) Microscopic origins of macroscopic properties of silicate melts and glasses at ambient and high pressure: implications for melt generation and dynamics. Geochim. Cosmochim. Acta 69, 3695–3710. Lee S. K., Cody G. D., Fei Y. and Mysen B. O. (2004) Nature of polymerization and properties of silicate melts and glasses at high pressure. Geochim. Cosmochim. Acta 68, 4189–4200. Lee S. K., Cody G. D., Fei Y. and Mysen B. O. (2006) The effect of Na/Si on the structure of sodium silicate and aluminosilicate glasses quenched from melts at high pressure: a multi-nuclear (Al-27, Na-23, O-17) 1D and 2D solid-state NMR study. Chem. Geol. 229, 162–172. Manghnani M. H., Sato H. and Rai C. S. (1986) Ultrasonic velocity and attenuation measurements on basaltic melts to 1500 C: role of composition and structure in the viscoelastic properties. J. Geophys. Res. 91, 9333–9342. Matsui M. (1994) A transferable interatomic potential model for crystals and melts in the system CaO–MgO–Al2O3–SiO2. Mineral. Mag. 58A, 571–572. Matsui M. (1996a) Molecular dynamics study of the structures and bulk moduli of crystals in the system CaO–MgO–Al2O3–SiO2. Phys. Chem. Mineral. 23, 345–353. Matsui M. (1996b) Molecular dynamics simulation of structures, bulk moduli, and volume thermal expansivities of silicate liquids in the system CaO–MgO–Al2O3–SiO2. Geophys. Res. Lett. 23, 395–398. Maumus J., Bagdassarov N. and Schmeling H. (2005) Electrical conductivity and partial melting of mafic rocks under pressure. Geochim. Cosmochim. Acta 69, 4703–4718. Miller G. H., Stolper E. M. and Ahrens T. J. (1991) The equation of state of a molten komatiite 1. Shock wave compression to 36 GPa. J. Geophys. Res. 96, 11,831–11,848. Murase T. and McBirney A. R. (1973) Properties of some common igneous rocks and their melts at high temperatures. Geol. Soc. Am. Bull. 84, 3563–3592. Nelson S. A. and Carmichael I. S. E. (1979) Partial molar volumes of oxide components in silicate liquids. Contrib. Mineral. Petrol. 71, 117–124. Neuville D. R., Courtial D., Dingwell D. B. and Richet P. (1993) Thermodynamic and rheological properties of rhyolite and andesite melts. Contrib. Mineral. Petrol. 113, 572–581. Nevins D. and Spera F. J. (1998) Molecular dynamics simulations of molten CaAl2Si2O8: dependence of structure and properties on pressure. Am. Mineral. 83, 1220–1230. Ohtani E. and Maeda M. (2001) Density of basaltic melt at high pressure and stability of the melt at the base of the lower mantle. Earth Planet. Sci. Lett. 193, 69–75. Ohtani E., Suzuki A. and Kato T. (1993) Flotation of olivine in the peridotite melt at high pressure. Proc. Jpn. Acad. Ser. B 69, 23–28. Ohtani E., Suzuki A. and Kato T. (1998) Flotation of olivine and diamond in mantle melt at high pressure: implications for fractionation in the deep mantle and ultradeep origin of diamond. Geophys. Mono. 101, 227–239. Poe B. T., McMillan P. F., Rubie D. C., Chakraborty S., Yarger J. and Diefenbacher J. (1997) Silicon and oxygen self-diffusivities in silicate liquids measured to 15 GPa and 2800 K. Science 276, 1245–1248. Poe B. T., Romano C., Zotov N., Cibin G. and Marcelli A. (2001) Compression mechanisms in aluminosilicate melts: Raman and XANES spectroscopy of glasses quenched from pressures up to 10 GPa. Chem. Geol. 174, 21–31. Potuzak M. and Dingwell D. B. (2006) Temperature-dependent thermal expansivities of multicomponent natural melts between 993 and 1803 K. Chem. Geol. 229, 10–27.

4556

B. Guillot, N. Sator / Geochimica et Cosmochimica Acta 71 (2007) 4538–4556

Potuzak M., Solvang M. and Dingwell D. B. (2006) Temperature independent thermal expansivities of calcium aluminosilicate melts between 1150 and 1973 K in the system anorthite– wollastonite–gehlenite (An–Wo–Geh): a density model. Geochim. Cosmochim. Acta 70, 3059–3074. Presnall D. C., Simmons C. L. and Porath H. (1972) Changes in electrical conductivity of a synthetic basalt during melting. J. Geophys. Res. 77, 5665–5672. Rai C. S. and Manghnani M. H. (1978) Electrical conductivity of ultramafic rocks to 1820 Kelvin. Phys. Earth Planet. Int. 17, 6–13. Reid J. E., Poe B. T., Rubie D. C., Zotov N. and Wiedenbeck M. (2001) The self-diffusion of silicon and oxygen in diopside (CaMgSi2O6) liquid up to 15 GPa. Chem. Geol. 174, 77–86. Rigden S. M., Ahrens T. J. and Stolper E. M. (1984) Densities of liquid silicates at high pressures. Science 226, 1071–1074. Rigden S. M., Ahrens T. J. and Stolper E. M. (1988) Shock compression of molten silicate: results for a model basaltic composition. J. Geophys. Res. 93, 367–382. Rivers M. L. and Carmichael I. S. E. (1987) Ultrasonic studies of silicate melts. J. Geophys. Res. 92, 9247–9270. Rowan L. R., Ahrens T. J. and Stolper E. M. (1991) Equation of state of molten basalt. EOS 72, 548–549. Rubie D. C., Ross, II, C. R., Carroll M. R. and Elphick S. C. (1993) Oxygen self-diffusion in Na2Si4O9 liquid up to 10 GPa and estimation of high-pressure melt viscosities. Am. Mineral. 78, 574–582. Saika-Voivod I., Sciortino F. and Poole P. H. (2004) Free energy and configurational entropy of liquid silica: fragile-to-strong crossover and polyamorphism. Phys. Rev. E 69, 0415031-13. Sakamoto D., Yoshiasa A., Yamanaka T., Ohtaka O. and Ota K. (2002) Electric conductivity of olivine under pressure investigated using impedance spectroscopy. J. Phys. Condens. Matter 14, 11375–11379. Secco R. A., Manghnani M. H. and Liu T. (1991a) The bulk modulus—attenuation—viscosity systematics of diopside–anorthite melts. Geophys. Res. Lett. 18, 93–96. Secco R. A., Manghnani M. H. and Liu T. (1991b) Velocities and compressibilities of komatiitic melts. Geophys. Res. Lett. 18, 1397–1400. Shimizu N. and Kushiro I. (1984) Diffusivity of oxygen in jadeite and diopside melts at high pressures. Geochim. Cosmochim. Acta 48, 1295–1303. Shiraishi Y., Ykeda K., Tamura A. and Saitoˆ T. (1978) On the viscosity and density of the molten FeO–SiO2 system. Trans Jpn. Inst. Metal. 19, 264–274. Smith J. R. and Agee C. B. (1997) Compressibility of molten ‘‘green glass’’ and crystal-liquid density crossovers in low-Ti lunar magma. Geochim. Cosmochim. Acta 61, 2139–2145. Smith W. and Forester T. (1996) DL_POLY_2.0: a generalpurpose parallel molecular dynamics simulation package. J. Mol. Graph. 14, 136–141. Stebbins J. F. and McMillan P. (1989) Five- and six-coordinated Si in K2Si4O9 glass quenched from 19 GPa and 1200 C. Am. Mineral. 74, 965–968. Stebbins J. F. and Sykes D. (1990) The structure of NaAlSi3O8 liquid at high pressure: new constraints from NMR spectroscopy. Am. Mineral. 75, 943–946.

Stebbins J. F., Oglesby J. V. and Kroeker S. (2001) Oxygen triclusters in crystalline CaAl4O7 (grossite) and uncalcium aluminosilicate glasses: 17O NMR. Am. Mineral. 86, 1307–1311. Stolper E. M. and Ahrens T. J. (1987) On the nature of pressureinduced coordination changes in silicate melts and glasses. Geophys. Res. Lett. 12, 1231–1233. Suzuki A., Ohtani E. and Kato T. (1995) Flotation of diamond in mantle melt at high pressure. Science 269, 216–218. Suzuki A., Ohtani E. and Kato T. (1998) Density and thermal expansion of a peridotite melt at high pressure. Phys. Earth Planet. Inter. 107, 53–61. Suzuki A. and Ohtani E. (2003) Density of peridotite melts at high pressure. Phys. Chem. Mineral. 30, 449–456. Tinker D. and Lesher C. E. (2001) Self diffusion of Si and O in dacitic liquid at high pressures. Am. Mineral. 86, 1–13. Tinker D., Lesher C. E. and Hutcheon I. D. (2003) Self diffusion of Si and O in diopside–anorthite melt at high pressures. Geochim. Cosmochim. Acta 67, 133–142. Toplis M. J., Dingwell D. B. and Lenci T. (1997) Peraluminous viscosity maxima in Na2O–Al2O3–SiO2 liquids: the role of triclusters in tectosilicate melts. Geochim. Cosmochim. Acta 61, 2605–2612. Tossell J. A. and Cohen R. E. (2001) Calculation of the electric field gradients at ‘tricluster’-like O atoms in the polymorphs of Al2SiO5 and in aluminosilicate molecules: models for tricluster O atoms in glasses. J. Non-Cryst. Solids 286, 187–199. Tyburczy J. A. and Waff H. S. (1983) Electrical conductivity of molten basalt and andesite to 25 kilobars pressure: geophysical significance and implications for charge transport and melt structure. J. Geophys. Res. 88, 2413–2430. Vinet P., Rose J. H., Ferrante J. and Smith J. R. (1989) Universal features of the equation of state of solids. J. Phys. Condens. Matter 1, 1941–1963. Vollmayr K., Kob W. and Binder K. (1996) Cooling-rate effects in amorphous silica: a computer-simulation study. Phys. Rev. B 54, 15,808–15,827. Webb S. L. and Dingwell D. B. (1994) Compressibility of titanosilicate melts. Contrib. Mineral. Petrol. 118, 157–168. Webb S. L. and Courtial Ph. (1996) Compressibility of melts in the CaO–Al2O3–SiO2 system. Geochim. Cosmochim. Acta 60, 75– 86. Xue X., Stebbins J. F., Kanzaki M., McMillan P. F. and Poe B. (1991) Pressure-induced silicon coordination and tetrahedral structural changes in alkali silicate melts up to 12 GPa: NMR, Raman and infrared spectroscopy. Am. Mineral. 76, 8–26. Xue X., Stebbins J. F. and Kanzaki M. (1994) Correlations between 17O NMR parameters and local structure around oxygen in high-pressure silicates: implications for the structure of silicate melts at high pressure. Am. Mineral. 79, 31–42. Yarger J. L., Smith K. H., Nieman R. A., Diefenbacher J., Wolf G. H., Poe B. T. and McMillan P. F. (1995) Al coordination changes in high-pressure aluminosilicate liquids. Science 270, 1964–1967. Yasuda A., Fujii T. and Kurita K. (1994) Melting phase relations of an anhydrous mid-ocean ridge basalt from 3 to 20 GPa: implications for the behavior of subducted oceanic crust in the mantle. J. Geophys. Res. 99, 9401–9414. Associate editor: Bjorn Mysen