Viscosity and viscosity anomalies of model silicates and magmas: a numerical investigation

0.26 MB·English

by  M. Bauchy

Viscosityandviscosityanomaliesofmodel silicatesandmagmas: anumerical investigation
M. Bauchy, B. Guillot, M. Micoulaut∗, N. Sator
Laboratoire de Physique Théorique de la Matière Condensée,
Université Pierre et Marie Curie, 4 Place Jussieu, F-75252 Paris Cedex 05, France
(Dated:September17,2012)
We present results for transport properties (diffusion and viscosity) using computer simulations. Focus is
madeonadensifiedbinarysodiumdisilicate2SiO -Na O(NS2)liquidandonmulticomponentmagmaticliq-
2 2
uids (MORB, basalt). In the NS2 liquid, results show that a certain number of anomalies appear when the
2 systemisdensified:theusualdiffusivitymaxima/minimaisfoundforthenetwork-formingions(Si,O)whereas
1 thesodium atomdisplaysthreedistinctrégimesfor diffusion. Someofthesefeaturescanbecorrelated with
0 theobtained viscosity anomaly under pressure, thelatter being befairlywell reproduced from thesimulated
2 diffusionconstant. Inmodelmagmas(MORBliquid),wefindaplateaufollowedbyacontinuousincreaseof
theviscositywithpressure. Finally,havingcomputedbothdiffusionandviscosityindependently, wecandis-
n
cussthevalidityoftheEyringequationforviscositywhichrelatesdiffusionandviscosity. Itisshownthatit
a
canbeconsideredasvalidinmeltswithahighviscosity. Ontheoverall,theseresultshighlightthedifficultyof
J
establishingafirmrelationshipbetweendynamics,structureandthermodynamicsincomplexliquids.
1
3
PACSnumbers:

Keywords:
]
i
c
s
I. INTRODUCTION andBehrens,2009).Atconstanttemperature,theevolutionof
-
l viscosity with pressure is complexand dependson the com-
r
t positionofthesilicatemelt,inparticularonitsdegreeofde-
m Viscosityisoneofthekeypropertiesinfluencingtheoverall
polymerization. With highly polymerized melts (e.g. albite
. behaviorofmagmaticliquids. Itisafundamentalpropertyin and dacite) the viscosity decreases continuously with pres-
t
a EarthSciencesandhasthereforeledafternearlyfivedecades sure, though its behavior above the maximal pressure of in-
m ofintensiveresearchtoahugebodyofexperimentalandthe-
vestigation (≃7 GPa for albite and dacite) cannot be clearly
oreticalstudies and data base (Giordanoet al., 2008; Mysen
- anticipated. However, it is observedthat the decrease of the
d and Richet, 2005). Systematic investigations with compo-
viscositywithpressurebecomesmoregradualwhenthetem-
n sition and temperature in simple or complex silicates have
peratureisraisedandcanleadtoaplateauvalueathighpres-
o
been performed, and some generic trends have been iden-
c sure (e.g. jadeite). When the melt is less polymerized (e.g.
tified that have become quite popular. For instance, it has
[ albite-diopsidesystem),afteraninitialdecrease,andaccord-
been found that numerous systems were displaying an Ar-
ing to the temperature, the viscosity exhibits a plateau value
1
rhenius behavior (Micoulaut, 2010) with temperature of the
v formh =h exp[E /k T]whereE representstheactivation or a weakly pronounced minimum with the pressure. With
8 0 A B A depolymerizedmelts(e.g. diopsideandperidotite)adifferent
energyforviscousflow. Whenproperelyrescaledwitharef-
1 situationoccurs.Theviscosityfirstincreaseswithrisingpres-
5 erence temperature,Tg, at which the viscosity is 1012 Poise, sureandreachesamaximumvaluepriortodecreaseathigher
6 twocategoriesofliquidshavebeenidentified(Angell,1995) pressure, but higher the degree of depolymerization of the
. : those which remainArrhenius-likein the super-cooledliq-
1 meltmorepronouncedistheviscositymaximum. Moreover,
uid over the entire range of temperatures implying that E
0 A inamoregeneralway,thepressureevolutionoftheviscosity
2 does not depend on T, and those which display a curvature hasbeenfoundtobe(anti-)correlatedwiththatoftheoxygen
1 in the Arrhenius representation (a semi-log plot as a func-
andsilicondiffusioncoefficients(ShimizuandKushiro,1991;
: tion of inverse temperature) indicating that one has a tem-
v Poeetal.,1997;Reidetal.,2001;Tinkeretal.,2003),afind-
i peraturedependenceintheactivationenergyforviscousflow ingwhichhasbeeninterpretedasthesignatureofapressure-
X (Wang et al., 2006). Usually, such liquids are fitted quite
inducedstructuralrearrangement.
r successfullywithaTamman-Vogel-Fulcher(TVF)law(Path-
a Theunderstandingofchangesintheviscousflowwithap-
manathanandJohari,1990;GiordanoandDingwell,2003a,b)
plied pressure can be suitably investigated by molecular dy-
: h =h exp(A/(T−T ),AandT beingconstants.
0 0 0 namics(MD)simulations.Thereareseveralwaystocompute
There has been a continuing effort to measure viscosity viscosity. ReadersshouldrefertoAllenandTildesley(1987)
as function of pressure (Scarfe, 1973; Kushiro, 1976, 1978; foradetaileddiscussionandpresentation.Atequilibrium,the
1986; Kushiro et al., 1976; Scarfe et al., 1979; Brearley et computation of viscosity from MD simulations can be per-
al., 1986; Dunn and Scarfe, 1986; Moriet al., 2000; Suzuki formed by using the Green-Kubo (GK) formalism which is
et al., 2002, 2005, 2011; Reid et al., 2001; 2003; Tinker et based on the calculation of the stress auto-correlation func-
al.,2004;Liebskeetal.,2005;Ardiaetal.,2008;DelGaudio tion(BoonandYip,1980)givenby:
F(t)=hPab (t)Pab (0)i (1)
∗Correspondingauthor:[email protected] wherePab (t)istheba component(a ,b )=(x,y,z)ofthemolec-
2
ularstresstensordefinedby: a sodium disilicate (NS2) and a mid-ocean ridge basalt
(MORB). The molecular dynamics (MD) simulations were
Pab =(cid:229)N mivai vbi +(cid:229)N (cid:229)N Fiajribj a 6=b , (2) p1e9r9f6o)rmineNdPwTithEnthseemDbLl_e.POThLeYe2q.u0actioodnes(oSfmmitohtiaonndsfFoorrraetostmers,
i=1 i=1j>i weresolvedwitha time step of1-2fs(10−15 s) bytheleap-
whereFa isthecomponenta oftheforcebetweentheionsi frogalgorithm.
ij
LiquidNa O-2SiO hasbeensimulatedbyplacing666sil-
and j, rb and vb being the b componentof the distance be- 2 2
ij i icon, 666 sodium, and 1668 oxygen atoms in a cubic box
tween two atomsi and j, and the velocityof atom i, respec-
withperiodicboundaryconditions. Astheatomsarebearing
tively. Alternative forms for Pab exist (Allen and Tildesley, charges(ions),thelongrangecoulombicforcesareevaluated
1987)butthesearefoundtoleadtospuriouseffectswhenperi-
by a Ewald sum. The atoms interact via a two-body poten-
odicboundaryconditionsareconsidered(Lee,2007)orwhen
tial(Born-Majertype)parametrizedby Teter(2003, see also
the Ewald sum is notperformedin a neat way (Alejandreet
Cormacketal.,2003),whichwrites
al.,1995;Allenetal.,1994).Inthiscase,F(t)isfoundtonot
decayto zeroin thelong-timelimit, asexpected(Lee,2007; V (r)=A .exp(−r/r )−C /r6+zz /r, (4)
ij ij ij ij i j
Green, 1957; Kubo, 1957). One other possibility is to com-
pute the viscosity within the frameworkof Non-Equilibrium whereristhedistancebetweenatomsiand j,zi istheeffec-
Molecular Dynamics (Lees and Edwards, 1972). Here, the tive charge associated with the ion i, and where Aij, r ij and
ratioofshearstressto astrain rateiscomputedandextrapo- Cij are parametersdescribingrepulsiveanddispersiveforces
lated to the limit of zerodrivingforce (Borodinetal., 2009; betweentheionsiand j. Inthepresentreport,wedonotfocus
Fernandezetal. 2004;CherneandDeymier,1998). onstructureandthermodynamicpropertiesofliquidNS2sim-
Coming back to the GK formalism commonly used with ulatedwiththispotentialandforadetailedanalysisandcom-
MDsimulations,theviscosityiscalculatedfromthetimein- parisonwithexperimentaldatawereferthereadertoDuand
tegralofthestressauto-correlationfunctionF(t)(seeeqs. (1) Cormack(2004,2005).Noticethatatzeropressureand300K
and(2)), thedensityofthesimulatedNS2glassisequalto2.45g/cm3,
i.e. quiteclosetotheexperimentalvalue(Vaillsetal.,2001)
h = 1 ¥ F(t)dt. (3) 2.37g/cm3.
k TV Z The mid-ocean ridge basalt (MORB), a nine component
B 0
system (0.13 wt% K O, 2.94 wt% Na O, 11.87 wt% CaO,
2 2
Theviscosityofafewsilicatemeltshasbeenevaluatedby
7.77wt%MgO,8.39wt%FeO,1.15wt%Fe O ,15.11wt%
2 3
MD simulation in followingthis route. Ogawaet al. (1990)
Al O , 1.52wt% TiO and 50.59wt% SiO ), hasbeen sim-
2 3 2 2
were the first to evaluate the viscosity of a silicate melt (a
ulatedwithatwo-bodypotentialofthesameformaseq. (4).
sodiumdisilicate)byMDsimulationbuttheirresultsareun-
ThechemicalcompositionoftheMORBmeltandthepoten-
certain due to a poor statistics. For silica, Horbachand Kob
tial parametersare those given by Guillot and Sator (2007a)
(1999)determinedthetemperaturebehaviorandBarratetal.
in their simulationstudy of naturalsilicate melts (see Tables
(1997)studiedthe pressuredependenceathightemperatures
1and2therein). ThesimulatedsampleofMORBwas com-
butwithoutpointingoutexplicitlythattheirresultsshoweda posedof3,000atoms(555Si,12Ti,195Al,9Fe3+,60Fe2+,
viscosityminimumatabout20GPa. MorerecentlyLackset
126 Mg, 141 Ca, 63 Na, 3 K, and 1,836 O) contained in a
al.(2007)ininvestigatingbyMDthepropertiesofsilicateliq-
cubicboxwithperiodicboundaryconditions.Thecontrolpa-
uidsalongtheMgO-SiO joinfoundaviscosityminimumat
2 rameters(timestep,Ensemble,etc.) ofthesimulationarethe
about20GPainpuresilica,aviscosityminimumwhichdisap-
same as for liquid NS2. The density of MORB at liquidus
pearsprogressivelywithincreasingMgOcontent.Sincethen,
temperatureanditsequationofstate(i.e. densityversuspres-
these resultshave been confirmedby the MD studies of Ad-
sure) are well reproduced by the MD calculation (for a de-
jaoudetal. (2008,2011)onMg SiO melts,thoseofSperaet
2 4 tailed description of the properties of the MORB melt, see
al. (2011)onmoltenenstatiteandalsobyfirst-principlesMD
GuillotandSator2007a,b). Noticethatlongsimulationruns
calculationsof Karki and Stixrude (2010a,b)on liquid silica (≃10 ns or 107 time steps) were performed to reach a good
andonliquidenstatite. Wefinallymentiontherecentstudies
statisticswhenevaluatingtheviscositybytheGreen-Kuboin-
ofKarkietal. (2011)andGoshandKarki(2011)onanorthite
tegrand.
liquidandMg SiO liquid,respectively.
2 4 Beforepresentingresults,itisworthwhiletokeepinmind
somekeydatawhenoneisabouttocalculatetheviscosityofa
supercooledmelt. Inusingahighperformancecomputer,the
II. MODELINGSILICATES
longestMDrunthatcanbeperformed(inareasonabletime)
forasystemcomposedof≃3,000atomsisabout1µs(or109
A. Simulationdetails MDsteps).Toobserveonthistimescalethediffusionofatoms
inthesimulationbox,the meansquaredisplacementofeach
To complement the above MD studies and to bring some atomafter1µshastobe(onaverage)oftheorderof≃10Å2
new information, we have investigated in the present study (themeandistancebetweentwonearestneighborsisroughly
two silicate melts exhibiting comparable NBO/T ratio (≃ 3 Å in a melt, so the diffusive regime is reached when each
0.8 and 1) but with very different chemical compositions; atom has interchanged its position with a nearest neighbor).
3
Therefore,thesmallestdiffusioncoefficientthatcanbeeval- 100
uatedovera1µsMDrunisequalto≃2×10−14m2/s(from Na : Negodaev et al.
Einstein relation Dmin=hrm2ini/6tmax where hrm2ini=10 Å2 and 10 Na NNaa :: GJouhpntsao ann edt Kali.ng
t =10−6 s). InapplyingtheEyringtheorytoviscousflow O Si, O : Truhlarova et al.
max s 1 MD : Horbach et al.
,theviscosityh canbededucedfromthediffusioncoefficient 2/ Na MD : Present work
m
Dofnetworkformingionsthroughtheequation, c
5 0.1
-0 O
1
h = kBT (5) D. 0.01 Si
l D
0.001 Si
whereT isthetemperatureandl a jumpdistanceinthe vis-
A
cous flow. Experimentally it has been shown (Shimizu and 0.0001
0 2 4 6 8 10 12 14 16
Kushiro,1984)thatthisequationworkswellwithviscousliq-
4
10 /T
uids(e.g. silicateswithhighsilicacontents)providedthatthe
100
D value is assigned to that of O or Si atoms and l ≃2.8Å,
Si, O : Lesher et al.
the oxygen-oxygen mean distance in silicate melts. Using 10 Na : Henderson et al.
equ. (5), a value of D equal to 2×10−14 m2/s leads to a MD : Present work
O Na
corresponding viscosity af about 3000 Pa.s at 1473 K. We 2/s 1
m
can use the Maxwell relation for which the relaxation time
for viscous flow is related to the viscosity throughthe equa- -5 0c 0.1 Si
tion,t relax =h /G¥ ,whereG¥ ≃1010 Pa.sistheshearmod- D.1 0.01 O
ulusatinfinitefrequency(DingwellandWebb, 1989). So to
fullyrelaxaviscousmeltofviscosityh ≃104Pa.soneneeds 0.001
toperformaMDrunlongerthant =1µs,afindingcon- B
relax
0.0001
sistent with the evaluation based on the diffusion. Although 0 1 2 3 4 5 6 7 8 9 10 1112 13 14 15 16
qualitative, these estimates are usefulbecause they pointout 104/T
that only the low viscosity regime (i.e. h ≪ 104 Pa.s) can
beinvestigatedbyMDsimulation,thehighlyviscousregime Figure1: (Coloronline)A)ComputeddiffusionconstantsDNa,DSi
(104 -1012 Pa.s)beingunreachablewith usualnumericalre- andDO intheNS2liquidasafunctionofinversetemperature(blue
sources. In the present study our MD runs did not exceed curvesandsymbols),comparedtoexperimentaldataforDNa(Gupta
andKing,1967;Negodaevetal.1972;Johnsonetal.,1951)andD ,
≃10ns,whichmeansthatonlyaviscosityvaluesmallerthan Si
D (Truhlarova et al., 1970), and to the simulated values of D ,
≃10Pa.s(100Poise) couldbe investigated(noticethatmost O Na
D andD usinganalternativepotential(redcurvesandsymbols,
oftheMDstudiespublishedintheliteratureoftendonotex- Si O
Horbach et al., 2001). B) Computed diffusion constants D , D
Na Si
ceed1nsrunintervals).
and D in the MORB liquid as a function of inverse temperature
O
(bluecurvesandsymbols), comparedtoexperimentaldataforD
Na
(Hendersonetal.)andD ,D (Lesheretal.1996).
Si O
B. Diffusioncoefficients
Wefirstcomputethemean-squaredisplacementofatagged of Truhlarova et al. (1970) who determined the (Si,O) dif-
atomoftypea inthemelt,givenby fusionfroma SiO dissolvationreactionin a NS2 liquid. In
2
thisrespect,theTeterpotentialappearstobeveryaccuratefor
hr2(t)i= 1 (cid:229)Na h|r(t)−r(0)|2i , (6) reproducingthediffusionofelementsintheNS2melt. Were-
Na i i markthatthediffusioncoefficientsdisplayanArrhenius-like
i=1 dependencewhen 104/T>4with activationenergiesequalto
and extract from the dependence of hr2(t)i the long time 1.18eV,1.17eVand0.43eVforSi,OandNarespectively.
behavior where the dynamics becomes diffusive. Using the To be complete, one should stress that alternative poten-
Einsteinrelationlimt→¥ hr2(t)i/6t =D,onecanindeedhave tials (Kramer et al., 1991) lead to diffusion coefficients for
access to the diffusion constants D (i=Si, O, Na) from the (Na,Si,O)thatare found(Horbachetal., 2001)to beatleast
i
mean square displacement. These quantities are plotted for one order of magnitude lower than the experimental values
NS2 and MORB in Fig. 1 as a function of inverse temper- (for 104/T>5), and our own results on the MORB liquid as
ature. Because of the slowing down of the dynamicsat low well(seebelow), whereasthecorrespondingactivationener-
temperature, the computationof D is here restricted to T > gies also differ by a factor of ≃2 (Fig. 1). These findings
i
1500KwithNS2andT>1173KwithMORB.ForNS2,the illustrate how difficult the numerical reproduction of diffu-
diffusionofsodiumatomsisfoundtoberemarkablycloseto sion coefficients in silicate melts can be. This well-known
theexperimentaldataobtainedbyGuptaandKing(1967)and featurehasbeenhighlightedforthecaseofsilicabyHemmati
Negodaev et al. (1972), with a very good agreement some- andAngell(1997)whoshowedthatdiffusionishighlymodel
where around 104/T ≃6. For network forming ions (Si,O) dependentwith differencesbetween modelsincreasing up to
an excellent agreement is found with the experimental data at least three orders of magnitude when the liquid is cooled
4
on 1 reasonableagreementwithBockrisdata(≃1.65eV)measured
cti (Bockrisetal.,1954)inthetemperaturerange(1723-1373K
n
n fu0.8 or 5.8<104/T<7.3). Furthermore, in the temperature range
atio NS2 where calculated and experimental values can be compared
correl0.6 swmitahlleearcbhyoatfhaecrto(ir.eo.f≃5.28<th1an04t/hTe<e6x.p7e)ricmalecnutlaaltdedatavaolfuBesocakre-
o
ut ris et al. (1955). This deviation is surpring given the good
A0.4
ss MORB agreementobtainedfor the diffusioncoefficients(see Figure
Stre0.2 1),andgiventhefactthath andDcanbesimplyrelated(equ.
d 5).
e
z
ali In the case of MORB the temperature dependence of the
m 0
or calculated viscosity is found to be non Arrhenian and can
N 0,001 0.01 0.1 1 10 100 1000 be fitted accurately over a large temperature range by the
Time (ps) TVF equation with h ¥ =6.5× 10−3 Poise, A = 0.42 eV, and
T = 670 K. Experimentally, the viscosity of basaltic melts
0
Figure 2: (Color online) Stress auto-correlation function F(t) of a is known to be non Arrhenian, so the simulation reproduces
2000KliquidNS2asafunctionoftimeforr =2.5g/cm3,anda2273 this feature, but the simulated MORB is much less viscous
KMORBliquidatthesamedensity. than the real one (a factor of 50 at 1673 K, see Fig. 3), this
is why A and T values are slightly differentfrom those rec-
0
ommendedforbasaltsbytheregressionformulaofGiordano
down. So, the very good agreement between experimental andDingwell(2003a,b),and Giordanoet al. (2008). In fact
andsimulatedionicdiffusivitiesintheNS2 melt, allowsone the MORB melt simulated at 1273 K corresponds, as far as
tobereasonablyconfidentintheabilityoftheTeterpotential the viscosity is concerned, to the real melt at about 1673K,
to describe accurately other transport properties as the melt thistemperaturedeviationdiminishingathighertemperatures
viscosity. As for the MORB melt, the temperature depen- wheretheviscosityofanysilicatemelt(whateverisitscom-
dence of the O, Si and Na diffusion coefficientsis shown in position)tendstobeverylow(lessthan0.1Poiseabove3000
panelBofFig.1(amongallthenetworkmodifyingcationsin K).Moregenerally,theagreementbetweenexperimentaland
MORB,onlythediffusioncoefficientofNaisshownbecause calculatedviscositiesfor silicate melts is actuallydifficultto
it is the most mobile ion and because it is commonwith the obtainbecausetheviscosityseemstobeevenmoresensitive
NS2melt).NaisfoundtodiffusemoreslowlyinMORBthan to the details of the potential than the diffusion coefficients
inNS2whereasSiandOatomsarediffusingfaster. However, themselves (for a related discussion see Vuilleumier et al.,
the agreementbetweensimulated andexperimentalvaluesis 2009). Moreover, direct comparisons between experimental
poorerwithMORBthanwithNS2(seeFig. 1foracompar- andsimulatedvaluesofthe viscosityarescarce inthe litera-
ison). Actually, Na, Si and O atoms diffuse too rapidly in ture(e.g. Lacksetal.,2007;Adjaoudetal.,2008;Karkiand
thesimulatedMORB (byafactorof≃5forNa and≃50for Stixrude, 2010b, Spera et al., 2011) and are usually done at
Si and O ). So one may anticipate some deviation between veryhightemperature(>2500K)where the viscosity is very
calculatedandexperimentalviscositiesforthebasalticliquid, lowanditsexperimentalvaluebadlyknownorobtainedfrom
andabetteragreementforliquidNS2. uncontrolledhigh-temperatureextrapolations.
In the TVF equation, h ¥ represents the high temperature
limittosilicatemeltviscosity.Byanalyzingmeasuredviscos-
C. Viscosity ity curves for a great number of silicate liquids, it has been
shown that all silicate melts convergeto a commonvalue of
In Figure 2, we represent the stress auto-correlationfunc- viscosity at very high-T, about 10−2-10−3 Poise (Russell et
tion F(t) of NS2 and MORB at around zero pressure, com- al., 2003; Giordano et al., 2008, Zheng et al., 2011). But it
puted following equs. (1) and (2). As mentioned above, isdifficulttotestthisconjecturebecauseofa lackofviscos-
the stress auto-correlationfunction is indeed found to decay itydataforsilicateliquidsatsuper-liquidustemperatures(i.e.
to zero. Typical oscillations which are found at the subpi- for T ≥ 2000 K). In contrast, the simulation can shed some
cosecondtimescaleareassociatedwithfrequenciesintheTer- light on this topic because the high-T range is easily acces-
ahertz range and above [10-300 THz] corresponding to the sible in a numerical experiment. However it is worthwhile
main band found in the vibrational density of states. This toclarifywhatmeanstheconceptofahigh-temperaturelimit
band is attributed to stretching motions of Si-O and Na-O forasilicateliquid.Whenasilicatemeltisheatedisobarically
bonds(and other cation-oxygenbondsin MORB). Note that (e.g. underatmosphericpressure) it vaporizesincongruently
theseoscillationsaredampedasthepressureisincreased(not in forming a gaseous phase which is a complex mixture of
shown). Followingeq. (3), the time integrationof the stress atomsandmolecules(e.g. themajorvaporspeciesaboveliq-
auto-correlationfunctionleadsto theviscosity(Fig. 3). The uid silica are SiO, SiO , O and O, and about thirty species
2 2
viscosity of NS2 exhibits an Arrhenius-like behavior over a compose the vapor phase above a basaltic lava, for a review
broad temperature range (2<104/T<6). The calculated acti- see Schaefer and Fegley, 2004). Nevertheless this is not a
vationenergyforviscousflowisabout≃1.33eVwhichisina real problemas long as the resulting vaporpressure remains
5
low,i.e. thetemperatureisnottoohigh(e.g. thevaporpres- 1000
sure above basaltic lavas is 10−5-10−4 bar at 1900 K). This MORB (Villeneuve et al., exp.)
situationchangesdrasticallyathighertemperaturewherethe
100
vaporpressurebecomeshigh(e.g.≃10−2barat2500Kabove
NS2 (Bockris, exp.)
for a tholeiitic liquid). Thenthe vaporizedfraction becomes
rapidlynon-negligible(10%ormore)andthecompositionof ) 10
e
theresidualmeltdeviatessignificantlyfromtheoriginalone. s
oi
To avoid this problem, one only has to prevent vapor for- P
mation in applying on the liquid an isostatic pressure larger ( 1 NS2 (MD)
h
than the saturating vapor pressure. A convenientthermody-
namic path consists to follow (or to be as close as possible) MORB (MD)
0.1
theliquidbranchofthe liquid-vaporcoexistencecurveupto
the critical point. The critical point parameters for silicates
haveneverbeenmeasuredbecauseoftheveryhightempera-
0.01
turesinvolved(muchbeyond5000K,estimationsforSiO are 2 3 4 5 6 7 8
2
giveninGuissaniandGuillot(1994)andMelosh(2007)).By 104/T (K-1)
imposingthepressure,wehavebeenabletoevaluatethevis-
cosityofoursimulatedMORBalong(orcloseby)theliquid
Figure3:(Coloronline)SimulatedviscosityoftheNS2liquidatzero
branch of the liquid-vapor coexistence curve. However due pressure(filledredsquares),comparedtoexperimentaldata(redcir-
tothedensityfluctuationsgeneratedbytheMDcalculationit cles) fromBockris et al. (1955), together withsimulatedviscosity
becomesmoreandmoredifficulttopreservethecohesionof of the MORB liquid at zero pressure (filled black squares), com-
the liquid phase in rising the temperature because the liquid paredtoexperimentaldata(brokenblackcurve)fromVilleneuveet
thentendstotransformspontaneouslyintothecorresponding al.(2008).ThedottedlineisahightemperatureArrheniusfitforthe
vapor phase. So the highest temperaturereached in keeping MORBdatawhereasthesolidlineisaTVFfit(seetextforparame-
ters).
the cohesion of the liquid phase around 1 bar was 4673K, a
temperature at which the density is equal to 1.47 g/cm3 and
the viscosity is evaluated to 10−2.8 Poise. This latter value Pressure (GPa)
matches rather well the estimated range for h ¥ (10−2-10−3 0 5 10 20 30 35
Poise,seeabove). Itiscertainlypossibletopursuethecalcu-
lationathighertemperatureuptothecriticalpoint,butitisa 1 1.2e-08
tedioustasktolocateaccuratelythecriticalpointofourmodel -1a)
adnudrewweitlhealviqeutihdaNttSo2awfuethuarevewroerakc.hIendu4s5i0n0gKth,eastaemmepeprraotucere- 2/s) 0.1 8e-09 (PT
ahcto¥ swifthoyirc1hM0t−hO2e.R5lBiPqouaiinsdedd.NeWnSse2itayantitsitchaiespalcotrewittihacsaatl1tp.h6oe6incgto/cmcmomu3lodannbdveaitlasubevoiousf-t -5D (10 cm 0.01 04e-09 kpressibility
10−3 Poise or slightly less. Notice that in approaching the m
-4e-09 o
criticalpoint,thestructureofasilicateliquidisquitedifferent 0.001 C
ofthestructurefromthemeltneartheliquidus.Thenearcrit-
-8e-09
icalliquidisindeedfullydepolymerizedandiscomposedof
1.5 2 2.5 3 3.5 4 4.5
molecularclustersinvolvingvariousspecies(oxides). Soitis
r (gm/cm3)
expectedthat all silicate liquidswill exhibitsimilar flow be-
haviorinthecriticalregion,thatofaregularmolecularliquid.
Figure4: (Coloronline)Silicon(red),Oxygen(green),andsodium
We now follow both viscosity and diffusionwith pressure
(black)diffusionconstantsofaNS2liquidatT=2000Kasafunction
inourtwomelts(NS2andMORB).InFig. 4and5,weplot ofsystemdensity(bottomaxis)andpressure(topaxis). Rightaxis:
thebehaviorofthediffusionconstantswithdensityandpres- Isothermalcompressibilityoftheliquidcomputedfromtheequation
sure. Themostmobileatom,Na,isfoundtodisplaythesame ofstateat2000K.
behaviorin the two melts: its diffusioncoefficientis contin-
uously decreasing with the pressure. In the case of MORB
theothernetworkmodifyingcations(Ca,K,Mg,...) showthe 2002),germania(DreyfusandMicoulaut,2012),somesodo-
samebehaviorwiththepressureasNadoes(Fig. 5). aluminosilicates(Poeetal.,1997)andinwater(Erringtonand
The behavior of network forming ions (Si,O) is quite dif- Debenedetti,2001).
ferent as their diffusion coefficients go through a maximum Issuchananomalyalsovisibleintheevolutionofviscosity
value with the pressure. Furthermore this maximum is lo- withpressure? Actually,forNS2wedofindthath displaysa
cated at about 15 GPa in NS2 and about 50 GPa in MORB minimumofviscositybutinapressurerange4-10GPawhich
andit ismuchmorepronouncedin the NS2 meltthanin the does not coincide with the diffusivity maximum for Si and
basaltic composition. Notice that a diffusivitymaximumfor O (P =15GPa). ForMORB, the viscosityshowsnomini-
max
networkformingionsisalsoobservedwithsilica(Shelletal. mum,itisnearlyconstantupto5GPaandincreasesprogres-
6
10 of these data along the sequence of composition albite→
jadeite→rhyolite→andesite→basalts→albite -diopside
x 1−x
→jadeite -diopside →diopside→peridotite,whichcov-
x 1−x
Na ers a large NBO/T range (from 0 to 2.8), leads to a con-
)
s
2/ O trasted picture. In the initial pressure range (0-2 GPa), the
m slope (dh /dP) is negative for composition with NBO/T ≤
c T
5 1 1.2 and is positive for higher NBO/T value. Although the
-0 Si
1 NBO/T ratio is a good indicator of melt polymerization, its
(
D useforpredictingthepressureevolutionofviscosityrequires
somecaution(e.g. ToplisandDingwell,2004). Forinstance,
albiteandjadeitemeltsarebothconsideredasfullypolymer-
ized (NBO/T=0) but the pressure behavior of their viscosity
0.1
isdifferent(aviscositycurvepresentingaweakminimumlo-
0 5 10 15 20
catedabout2-3GPaforjadeiteascomparedwithacontinuous
Pressure (GPa)
decreaseoftheviscosityoverseveraldecadesupto7GPafor
albite,seeSuzukietal. 2002,2011). Inkeepinginmindthis
Figure5: (Coloronline)Silicon(red),Oxygen(green),andsodium
information,andaccordingtotheabovedataanalysis,thevis-
(black)diffusionconstantsofaMORBliquidatT=2273Kasafunc-
cosity of NS2 and MORB composition (NBO/T=1 and 0.8)
tionofpressure.
should exhibit a negative slope. For instance, Scarfe et al.
(1987) notice that the viscosity of the NS2 melt is decreas-
ing between 0 and 2 GPa. In the same way, the viscosity of
10 oceanic basalts (Kushiro et al., 1976; Kushiro, 1986; Ando,
2003)maintainedatfixedtemperaturedecreasesslightlywith
increasingpressureupto3GPawhereitreachesaminimum
valuepriorincreaseuponfurthercompression. Butasan in-
) crease of the temperature tends to decrease the viscosity of
e
s NS2 a meltmaintainedat a fixed pressure, the weak minimumof
oi
P viscosityinbasaltsisexpectedtodisappearinincreasingthe
( 1 temperature. Thisiswhya viscosityminimumisbarelyvis-
h ible with ourMORB melt simulated at 2273K,especially as
thelatteronedescribestheviscousflowofarealMORBata
highertemperature,maybeashighas2600-2700K(seeFig.
MORB 3).IncontrasttheviscosityofthesimulatedNS2meltat2000
Kpresentsaclearminimumasafunctionofpressure(notice
0.1 that the temperature dependenceof the viscosity of the NS2
0 10 20 30 meltismuchbetterreproducedbythesimulationthanin the
P (GPa) MORBcase,seeFig. 3).
Sincethefirstreportofananomalouspressuredependence
Figure 6: (Color online) Simulated viscosity of a NS2 liquid at oftheviscosity((dh /dP) <0)byKushiro(1976)andKushiro
T
T=2000 K (filled circles) and a MORB liquid at T=2273 K (open
et al. (1976), differentauthors have tried to explain the ori-
circles)asafunctionofpressure.
gin of this phenomenon. Using the Adam-Gibbstheory and
a pressure dependence for the degree of melt polymeriza-
tion,BottingaandRichet(1995)haveshownthatthedecrease
sivelyatahigherpressure.Thesesimplefindingssuggestthat
of the specific volume with alkali content was leading to an
thediffusivitiesofnetworkformingionsarenottheonlyones
increased sensitivity of the viscosity with pressure change.
responsible of the viscous flow, network modifying cations
Gupta (1987) has proposed that a system having a thermal
playalsoarole.
expansion coefficient in a glass larger than that of a liquid
coulddisplaysuchanomalies.However,thisdoesnotseemto
beverifiedexperimentally(Knocheetal., 1994). A relation-
III. DISCUSSION shipwithstructuralchangesunderpressurealsohasbeenpro-
posed,inparticularapressure-inducedcoordinationchangeof
As briefly recalled in the introduction, the pressure de- SiandAlorSi-OandAl-Obondweakeningduetopressure-
pendence of the viscosity has been reported experimentally induceddistortionoftheliquidstructure(Waff,1975;Mysen
for a number of silicates (Scarfe, 1973; Kushiro, 1976, etal.,1980).However,astudy(Sharmaetal.,1979)ondense
1978; 1986; Kushiro et al., 1976; Brearley et al., 1986; liquidGeO wherebothcoordinationchange(Micoulautetal.
2
Dunn and Scarfe, 1986; Mori et al., 2000; Suzuki et al., 2006)anda pressureinducedviscosity anomalyoccurcould
2002, 2005, 2011; Reid et al., 2003; Tinker et al., 2004; notestablish the correlation. This doesnot rule outthe pos-
Liebske et al., 2005; Ardia et al., 2008). The analysis sibility of a structural origin for the anomaly. For instance,
7
Suzukietal. (2002)dofindacorrelationbetweenthegradual
decreaseof the viscosity of albite with pressure upto 7 GPa T=2000K, change in P
20 P=0 GPa, change in T
andtheincreasingconcentrationof5-coordinatedAl.
Thevariationofh alonganisothermalsosuggeststhatthe
corresponding activation energy for viscous flow EA(P,T) is
minimum.Here,forselectedpressuresintheNS2system,we 15
have computed h with temperature and find that the activa- )
Å
tionforviscousflow shouldminimizeindeedin thepressure (

window, i.e. we obtain EA=1.33, 1.09 and 1.34 eV for P=0, l 10
6, and 22 GPa, respectively. The link between E minima
A
and the presence of a stress-free state in correspondinglow-
temparature glassy networks has been established recently
(Micoulaut, 2010), indicative of a flexible to rigid transition 5
(MicoulautandPhillips,2003). Sinceitisknownthatatam-
bientpressuretheNS2isflexible(Vaillsetal.,2005),onemay
0.1 1 10 100
obtain a pressure-inducedflexibleto rigid transitionwith in- h (Poise)
creasingP(Trachenkoetal. 2004),theminimuminE (and
A
h )beingoneoftheexpectedsalientfeatures. Figure 7: (Color online) Typical lengthscale l associated with the
Insummary,allthesestudiespointoutthesubtleinterplay Eyringrelation(5)asafunctionofthemeltviscosityofaNS2liq-
betweenstructure,dynamicsandthermodynamicsofthesili- uid, computed either from thepressure variation(at 2000 K,black
catemeltthatleadstoanomalousbehaviorintransportprop- symbols)orfromthetemperaturevariation(atP=0,bluesymbols).
erties.
T=2000Kfor NS2 (with varyingpressure). For convenience
A. Viscosityversusdiffusion this length is represented in Fig. 7 as function of the calcu-
latedmeltviscosity. Alongeachpath(isobarorisotherm)the
arrowindicatesthedirectionoftheincreasingthermodynamic
Theconnectionbetweendiffusionandviscositycanbere-
parameter(TorP).Noticethatintherepresentationl (h )(Fig.
alizedinasimplewayviatheEyringrelationshipgivenbyeq.
7), both curvesexhibita maximumvalue (22Å) for both the
(5)involvingl whichisatypicaljumpdistanceforthediffus-
isothermandtheisobarwhichexpressesanon-trivialbehavior
ingatom. AccordingtoEyring,therelationship(5)holds(i.e.
of h with diffusivity. Clearly, for the two meltcompositions
l is constant and equal to the average interatomic distance)
andwhateverthethermodynamicpath,l decreaseswhenthe
iftheactivatedprocessfordiffusioncanbeassumedidentical
meltviscosityincreases. Asimpleextrapolationsuggeststhat
with that of viscous flow. So for silicate melts the network
theEyringrelationshouldbeverifiedonlywhentheviscosity
formingions(SiandO)aregoodcandidatestofulfillthispre-
issufficientlyhigh(forh > 100Poise),i.e. whenl becomes
requisite.InfacttheEyringrelationshiphasbeentestedovera
of the orderof ≃ 5Å.. This conclusionis in agreementwith
largerangeofmeltcomposition,fromsimplebinaryoxidesto
theviscositydatadiscussedabove. On theotherhand,when
naturalmagmacompositions(Oishietal., 1975;Yinnonand
the viscosity is becoming very low and the diffusion coeffi-
Cooper,1980;Dunn,1982;ShimizuandKushiro,1984;Dunn
cientsofnetworkformingionshigh(>10−9m2/s),thelength
andScarfe,1986;Lesheretal., 1996;Reidetal., 2003;Tin-
l matchesbetter the valuespredicted by the Stokes-Einstein
keretal.,2004). Inpolymerizedmeltsofhighviscosity(e.g.
equation, l = 2p d, where d is the diameter of the diffusing
jadeite),itisobservedthatthediffusivityofoxygenionsand
molecule(≃ 3Åforoxygen).AstheStokes-Einsteinrelation
themeltviscosityarerelatedtoeachotherthroughtheEyring
workswellwithmolecularliquids(LiandChang,1955)itis
relation in introducing for the jump distance, l , the average
not too surprising that it also leads to reasonable values for
oxygen-oxygendistanceinthemelt(≃2.8Å).Indepolymer-
silicateliquidsatveryhightemperature(e.g. inapproaching
izedmeltswheretheviscosityisgenerallylower(e.g. basalts
thecriticalpoint)wherethemeltisessentiallydepolymerized
anddiopside),theEyringrelationisfulfilledwithl largerthan
andsharessomestructuralsimilaritieswithregularliquids.In
2-3timesthemeanoxygen-oxygendistance.Howeverevenin
summary,theEyringrelationisnotadaptedtothedescription
polymerizedmelts,theEyringrelationfailstofullyreproduce
of the transport properties of silicate melts at superliquidus
thecorrelatedevolutionofviscosityandoxygenself-diffusion
temperatures. Alternative modelshave been proposedin the
coefficient. Asamatteroffactithasbeenreportedthatoxy-
literature(seeMungall,2002)butaquantitativetheoryisstill
gen self-diffusion coefficient in albite (Poe et al., 1997) and
lacking.
indacitemelts(TinkerandLesher,2001)reachesamaximum
valueatabout5GPawhereastheviscosityofthecorrespond-
ing liquids is continuously decreasing with the pressure and
shows no viscosity minimum at 5 GPa (Suzuki et al., 2002; B. Structuralcorrelations
Tinker et al., 2004). In this context, we have evaluated the
length, l = k T/h D along two thermodynamic paths : Asmentionedbefore,therehavebeenmanyattemptstocor-
B O,Si
the isobar P=0 (with varying temperature) and the isotherm relate structuralfeatureswith the viscosity anomaliesin sili-
8
cateliquids(Waff,1975;WoodcockandAngell,1976;Mysen 100
a
et al., 1980, Suzukiet al., 2002). From the liquid structures IV
Si
obtainedforthetwomeltcompositions(NS2andMORB),we 80
%)
h6oa-nfveoelsdceoecmsootphruadtteindtahtteheedm.feraRltcetssitournultcsotfuarsreieliedcvoisonplvlaaetyosemfdrsoitmnhaFrtoiagur.egh48l-ay,.51-H0a0en%rde bility ( 60 SiV SiVI
4-coordinatedSiIV to a growing proportionof SiV and SiVI, oba 40
SiV increasingsubstantiallyforPlargerthan10GPaandgo- Pr
20
ing througha maximumfoundatabout25 GPa in NS2 (and
abovethispressureforMORB)whereasSiVI isbecomingsig-
0
nificantonlyabove15GPa. 0 10 20 30 40 50
From the probability distribution {p} of Sii species -1.1
i
b
(i=IV,V,VI), one can compute a configurational entropy ac-
-1.2
cordingto:
10
S =−(cid:229) p lnp, (7) e) -1.3
c i i i Pois (k)2B
which maximizeswhen the populationof SV is maximum ( -1.4s
h
(≃30GPa).ThisconfigurationalentropyisrepresentedinFig.
-1.5
1
8reblaatsioenxfpo[1r/vSics]c.oIsnityfaocft,thifeofnoermas:shum=ehstheaxtpt[hAe/TASda]mis-Gvaiblibds, exp[1/Sc]
0 c -1.6
one may expectfrom a simple inspection of the relationship 0 10 20 30 40 50
thata minimumin h iscorrelatedwith a maximumin S , as 22
c 150
also suggestedby Bottingaand Richet(1995). However, for
NS2 a quick look at figure 8b clearly shows that h (P) and 145 20
Sexcpd[o1e/Sscn]oatrceoninocticdoerrweiltahtetdheasptrheessluorceartiaonngoefwthheermeaaxviimscuomsitiyn deg) 140 18eg)
minimumisobserved. InthecaseofMORB,thegradualin- ( 135 (d
creaseoftheviscositywithpressurecouldbecorrelatedwith q>< 130 16sq
that of S . So it is hard to firmly conclude that a structural
c
reorganizationthroughcoordinationchangesofSiisdirectly 125 c 14
responsible of the pressure evolution of the viscosity. To be
120
complete,noticethatGoeletal. (2011)haverecentlyshown 0 10 20 30 40 50
P (GPa)
byMDsimulationsthatthediffusioncoefficientofSiatomsin
silicate melts alongthe MgO-SiO join, correlateswell with
2
the excess entropywhen the latter one is evaluated from the Figure8:(Coloronline)a)Fractionof4-,5-and6-foldsiliconatoms
asafunctionofpressureintheNS2(blackcurves)andMORBliquid
knowledgeofthepairdistributionfunctions. Butitisunclear
(redcurves). b)SimulatedviscosityoftheNS2liquidatT=2000K
ifthiscorrelationalsoholdswiththeviscosityastherelation-
(filledblackcircles, sameasFig. 3), compared toanAdam-Gibbs
shiplinkingdiffusivityandviscosityisfartobeobvious(see
viscositycomputedfromastructure-relatedconfigurationalentropy
theabovediscussion).
Sc. Rightaxis: Entropys2(brokenline)computedfromequ. (8). c)
We have therefore calculated an excess entropy s2 in us- Meananglearoundthebridgingoxygen(BO)atomandcorrespond-
ingtheexpressionofthetwo-bodyterm(BaranyaiandEvans, ingsecondmoment(standarddeviation)ofthebondangledistribu-
1989),whichwrites: tion(rightaxis).
kBr (cid:229) (cid:229)
s =− xx g lng −(g −1) dr (8)
2 2 i jZ (cid:20) ij ij ij (cid:21)
i j
afurtherscrutiny.Analternativewaytounderstandthetrends
r beingthedensity,andx theconcentrationofspeciesi. Re- of the viscosity with pressure may arise from the evolution
i
sults ofthis calculationare shownin Fig. 8b (rightaxis). s of the T-O-T intertetrahedral bond angle distribution (BAD)
2
is foundto decrease continuouslywith pressureand a corre- asitiswellknownthatbondanglesaremuchmorepressure
spondingAdam-Gibbsh doesnotshowaviscosityminimum. sensitive than bonds themselves (Micoulaut, 2004). Figure
ThisbehaviorisatvariancewiththatobtainedbyGoeletal. 8c shows the evolution of the average value of the BAD be-
(2011)intheirMDsimulationstudyofMgO-2SiO (asystem tweenabridgingoxygen(BO)anditstwoSiO tetrahedra,
2 4/2
related to NS2) where S is exhibitingan entropymaximum togetherwith the secondmoment(standarddeviation)of the
2
around 10 GPa, very well correlated with a diffusivity max- BAD.Thelatteronequantifiestheangularexcursionarounda
imum, as also found by Jabes et al. (2010) for liquid SiO , meanvalueanditisexpectedtobepressuresensitive,ashigh-
2
GeO ,BeFe andH O.Sooursearchforaneventualrelation- lighted in the case of densified GeSe (Bauchyet al., 2011).
2 2 2 2
shipbetweenviscosityandconfigurationalentropyatthesin- WhentheaveragevalueoftheBADdecreasessteadilyinNS2
gletemperature2000Kisinconclusiveandcertainlydeserves meltfrom150oto140obetweenP=0and20GPa,itvariesbe-
9
tween140oand128oinMORBoverthesamepressurerange. shouldholdonlyfordeeplysupercooledliquidsasitfailsfor
In the same way the angular excursion increases (with pres- liquidustemperatures.
sure) equally well in the NS2 and in MORB melt as shown Finally,wehaveattemptedtorelatetheobtainedanomalies
bythestandarddeviations q evolutionwithpressure.So,here inviscositywithstructuralchangesofthemelt. Withgrowing
again, it is not clear that viscosity and diffusion anomalies pressure, the population of higher coordinated silicon (SiV,
(e.g.inNS2)inducedbythepressurecanbeunivoquelyasso- SiVI) increases in a similar fashion in the NS2 and MORB
ciatedwithaparticularstructuralrearrangementevenifthere liquid which rules out the possibility of a direct relationship
isnodoubtthatstructuralchangesaretakingplaceinthesili- between e.g. coordination increase and viscosity minimum.
catemeltunderinvestigation. Moreover, in the NS2 liquid configurational entropies have
beencalculatedandtheirbehaviorwithpressuredoesnotfol-
low at all the trend obtained for the viscosity. More subtle
IV. SUMMARYANDCONCLUSION structuralandthermodynamicchangesareclearlyatplayand
onedoesnotrecoverforNS2thereportedcorrelationbetween
configurationalentropyanddiffusivitymaximafoundbyGoel
In summary, we have used Molecular Dynamics simula-
etal. (2011)intheMgO-SiO system.
tionsto studythe dynamicsof two systemsof importancein 2
Our attempts to find out a quantitative link between
materials science and geology : the sodium disilicate (NS2)
pressure-induced structural rearrangement and viscosity
and a basaltic (MORB) liquid. Results show that the diffu-
anomaliesisratherdeceptive.Althoughsomecorrelationsbe-
sion constant computed from a two-body Born-Majer type
tween diffusivity, viscosity and excess entropy certainly do
potential at ambient conditions exhibits a remarkable agree-
existandaresometimesemphasizedintheliterature,theiref-
mentwithexperimentaldataforallspecies(Na,Si,O)inthe
fectivenessdependsstronglyon the meltcomposition, a fact
NS2 liquid. This contrasts with the level of agreement ob-
which makes risky a clear conclusion. On the other hand, a
tained for the MORB system which shows a deviation by at
recenttheoreticalanalysisby Schmelzer et al. (2005)seems
least one order of magnitude when compared to experimen-
tobepromising. Inthisstudythepressuredependenceofthe
tal findings. The same situation is nearly reproduced in the
viscosityis analyzedfroma thermodynamicalpointofview.
comparisonofcalculatedviscosityusingtheGreen-Kubofor-
Theworkingequationexpressesthederivativeoftheviscosity
malism (based on the stress auto-correlation function) with
withrespecttothepressureunderthefollowingform,
experimental data from Bockris et al; (1955) for NS2, and
Villeneuveetal. (2008)forMORB. dh k h¶ h¶ dx
T
Wehavethenstudiedthebehaviorofboththediffusioncon- (cid:18)dP(cid:19) =−a (cid:18)¶ T(cid:19) +(cid:18)x¶ (cid:19) (cid:18)dP(cid:19) , (9)
stantandtheviscositywithchangingpressureinthetwoliq- T T P T T
uidswithtemperaturecloseto 2000K. IntheNS2 system,a wherex isanorderparameter(e.g. thedegreeofpolymer-
diffusivitymaximumfor the networkformingspecies (Si,O) izationgivenbyBO/(BO+NBO)whereBOandNBOarethe
isobtainedataround15GPa,aresultwhichparallelssimilar numbers of bridging and nonbridgingoxygen, respectively),
findingsforother densifiedtetrahedralliquidssuch aswater, k isthecompressibilityanda thethermalexpansioncoef-
T T
germania, or silica. The sodium cation is found to display ficient. If the order parameter is invariant with pressure (so
threedistinctiverégimesfordiffusion,theonecorresponding dx /dP=0)then dh /dPis governedby the first term in eq. (9)
to lower densities (negativepressure)being identified with a whichexpressesa freevolumeeffect(h onlyvarieswith the
stretchedmelt(BauchyandMicoulaut,2011). Thesefeatures total volume of the melt). But as k is negative in a liquid
T
arenotfoundinthemoredepolymerizedMORBliquidwhich at equilibrium,and ¶h¶ T is positive(the viscosity decreases
showsforSiandOdiffusionasteadybehaviorbetween0and whenthetemperatureincreases),thevariationofh withPwill
5 GPa, followed by a continuous decrease applied pressure. bepositiveornegativeaccordingtothesignofthethermalex-
Interestingly, a pressure window is found in the NS2 liquid pansioncoefficient, a . If the latter oneis negativethen the
T
withaminimumataround5GPathatdoesnotseemtobecor- viscositywilltendtodecreaseinincreasingpressureandifa
T
relatedwiththemaximumintheobtaineddiffusivity. Again, becomespositiveafterafurthercompressionthentheviscos-
this contrastswith the steadily increase of the MORB liquid itywillgotoaminimumvaluebeforetoincrease. Asshown
viscositywithpressure. bySchmelzeretal. (2005)thisbehaviorisencounteredwith
Having in hand both the computed diffusion and the vis- theviscosityofliquidwater. Howeverthefirsttermineq. (9)
cosity, we have investigated the validity of the Eyring equa- isnotsufficienttoreproduceallthemagnitudeoftheviscosity
tion which relates both quantities and involves also a jump minimumandthe contributionof structuralchangesinduced
distancel usuallyoftheorderoftheaverageoxygen-oxygen bythepressurehastobeaccountedfor(viathesecondtermin
bond distance (2.8 Å), as exemplified by numerous experi- eq.(9)).Inthecaseofsilicatemeltsitiswelldocumentedthat
mentalstudies. Thisinvestigationhas beenperformedalong theviscosityincreaseswhenthedegreeofpolymerizationof
twothermodynamicpaths:atfixedpressure(P=0)andchang- the melt increases(so dh /dx >0)wherasfor a givenmelt its
ing temperature, and at fixed temperature (T≃ 2000 K) and degree of polymerization decreases when pressure increases
changingpressure. Wehaveshownthatl calculatedfromthe (so dx /dP<0). Therefore the second term of eq. (9) is gen-
simulateddiffusionandviscositywasindeedoftheorderofa erally negative. For NS2 and MORB we have checked that
few Angstroems, butonly in the high viscosity régime (here a ispositiveatliquidtemperatureandoveralargepressure
T
100Poise).ItunderscoresthefactthattheEyringrelationship rangewhichimpliesthatthefirsttermineq.(9)ispositivefor
10
thesetwocompositions. AsourMDresultsforNS2leadtoa Bauchy M., MicoulautM. (2011). From pockets to chan-
decreaseoftheviscositybetween0and≃4GPa,thissuggests nels: density controlleddiffusionin silicates. Phys. Rev. B
that the second term in eq. (9) overbalances the first term. 83,184118.
ForMORBtheviscosityisfoundvirtuallyconstantinthe0- BauchyM.,MicoulautM.,CelinoM.,BoeroM.,LeRoux
4 GPa range so the two terms nearly compensate each other S., MassobrioC. (2011). Angularrigidityin tetrahedralnet-
in eq. (9). At high pressures(e.g. above10 or 20 GPa) the work glasses with changingcomposition. Phys. Rev. B 83,
two melts are more and more depolymerized and dx /dP be- 054201.
comessmallbecausethedegreeofdepolymerizationdoesnot Behrens H., Schulze F. (2003). Pressure dependence of
evolveverymuchwiththepressureduetomeltdensification. melt viscosity in the system NaAlSi O -CaMgSi O . Am.
3 8 2 6
Sothefreevolumeeffect(firstterm)isthedominantcontribu- Mineral.88,1351.
tiontoeq. (9)athighcompressionratesandtheviscosityof Bockris J.O’M., Mackenzie J.D., Kitchener J.A. (1955).
thetwomeltsincreasessteadilyinthehighpressurerange.In Viscousflowinsilicaandbinaryliquidsilicates. Trans. Fara-
summaryitappearsthateq. (9)isveryusefultodecipherthe daySoc. 51,1734-1748.
complexbehaviorexhibitedbytheviscosityofsilicatemelts BoonJ.P.,YipS.(1980. MolecularHydrodynamics. Mc-
asfunctionofpressure.Tobemorequantitativeandpredictive GrawHill.
itremainstomakethelinkbetweenthestructuralevolutionof Borodin O., Smith G. B., Kim H. (2009). Viscosity of a
themeltandthevariationoftheorderparameterwiththepres- RoomTemperatureIonicLiquid: PredictionsfromNonequi-
sure(x (P)).Thiswillbedoneinafuturepublication. libriumandEquilibriumMolecularDynamicsSimulations.J.
Phys.Chem. B113,4771.
BottingaY.,RichetP.(1995). Silicatemelts: The"anoma-
Acknowledgments lous" pressure dependence of the viscosity. Geochim. Cos-
mochim.Acta59,2725-2731.
TheauthorsthankG. Henderson,D.R. Neuville, P. Richet BrearleyM.,DickinsonJr. J.E.,ScarfeC.M.(1986). Pres-
forusefuldiscussions. sure dependence of melt viscosities on the join diopside-
albite. Geochim.Cosmochim.Acta50,2563-2570.
CherneF.J.,DeymierP.A.(1998).CalculationofViscosity
of Liquid Nickel byMolecular Dynamics Methods. Scripta
Mater39,1613.
AdjaoudO.,Steinle-Neumann,JahnS.(2008).Mg2SiO4liq- Cormack A. N., Du J., Zeitler T. R. (2003). Sodium ion
uid under high pressure from molecular dynamics. Chem. migrationmechanismsinsilicateglassesprobedbymolecular
Geol. 256,185-192. dynamicssimulations.J.Non-Cryst.Solids323,147.
Adjaoud O., Steinle-Neumann, Jahn S. (2011). Transport DelGaudioP., BehrensH.(2009). Anexperimentalstudy
propertiesofMg2SiO4liquidathighpressure: Physicalstate on the pressure dependence of viscosity in silicate melts. J.
ofamagmaocean.EarthandPlanet. Sci. Lett. 312,463-470. Chem. Phys.131,044504-1-14.
AbramoM.C., CaccamoC., PizzimentiG.(1992). Struc- DingwellD.B.,WebbS.L.(1989). Structuralrelaxationin
tural properties and medium-range order in metasilicate silicate melts and non-newtonian melt rheology in geologic
(CaSiO )glass: Amoleculardynamicsstudy.J.Chem. Phys. processes.Phys. Chem. Minerals16,508-516.
3
96,9083. DreyfusC.,MicoulautM.(2012),unpublished
Alejandre J., Tildesley D., Chapela G. A. (1995). Molec- DuJ.,CormackA.N.(2004.) Themediumrangestructure
ularDynamicssimulationoftheorthobaricdensitiesandsur- ofsodiumsilicateglasses: amoleculardynamicssimulation.
facetensionofwater. J.Chem. Phys.102,4574. J.Non-Cryst.Solids349,66.
AllenM.P.,TildesleyD.J.(1987). ComputerSimulations DuJ.,CormackA.N.(2005). MolecularDynamicsSimu-
ofLiquids,OxfordUniv.Press,p. 64. lationoftheStructureandHydroxylationofSilicaGlassSur-
Allen M. P., Brown D., Masters A. J. (1994). Use of the faces. J.Am. Ceram. Soc. 88,2532.
McQuarrieequationforthecomputationofshearviscosityvia Du J., Corrales L. R. (2006). Compositional dependence
equilibriummoleculardynamics-comment.Phys.Rev.E49, ofthefirstsharpdiffractionpeaksinalkalisilicateglasses: A
2488. moleculardynamicsstudy.J.Non-Cryst.Solids,352,3255.
AndoR., OhtaniE., SuzukiA., RubieD.C., FunakoshiK. Dunn T. (1982). Oxygen diffusion in three silicate melts
(2003). Pressuredependencyof viscosities of MORB melts. along the join diopside-anorthite. Geochim. Cosmochim.
AGUFallMeeting2003,abstractV31D-0958. Acta46,2293-2299.
Angell C.A. (1995). Formation of Glasses from Liquids Dunn T., Scarfe C.M. (1986). Variation of the chemical
andBiopolymers.Science267,1924. diffusivity of oxygen and viscosity of an andesite melt with
Ardia P., Giordano D., Schmidt M.W. (2008). A model pressureatconstanttemperature.Chem. Geol. 54,203-215.
fortheviscosityofrhyoliteasafunctionofH2O-contentand ErringtonJ.R.,DebenedettiP.G.(2001).RelationshipBe-
pressure: A calibration based on centrifuge piston cylinder tween Structural Order and the Anomalies of Liquid Water.
experiments.Geochim.Cosmochim.Acta72,6103-6123. Nature410,259.
Barrat J.-L., Badro J., Gillet P. (1997). Strong to Fragile FernandezG.A.,VrabecJ.,HasseH.(2004).AMolecular
transitioninamodelofliquidsilica. Mol. Simul. 20,17-25. Simulation Study of Shear and Bulk Viscosity and Thermal