A thermodynamical model for the surface tension of silicate melts in contact with H2O gas

Available online at www.sciencedirect.com

ScienceDirect Geochimica et Cosmochimica Acta 175 (2016) 113–127 www.elsevier.com/locate/gca

A thermodynamical model for the surface tension of silicate melts in contact with H2O gas Simone Colucci a,⇑, Maurizio Battaglia b,c, Raffaello Trigila b a

Istituto Nazionale di Geofisica e Vulcanologia INGV-sezione di Pisa, via della Faggiola 32, 56126 Pisa, Italy b Sapienza – Universita` di Roma, Piazzale A. Moro 5, 00185 Roma, Italy c USGS – Volcano Science Center, 345 Middlefield Rd, Menlo Park, CA 94025, USA Received 28 April 2015; accepted in revised form 27 October 2015; Available online 11 December 2015

Abstract Surface tension plays an important role in the nucleation of H2O gas bubbles in magmatic melts and in the time-dependent rheology of bubble-bearing magmas. Despite several experimental studies, a physics based model of the surface tension of magmatic melts in contact with H2O is lacking. This paper employs gradient theory to develop a thermodynamical model of equilibrium surface tension of silicate melts in contact with H2O gas at low to moderate pressures. In the last decades, this approach has been successfully applied in studies of industrial mixtures but never to magmatic systems. We calibrate and verify the model against literature experimental data, obtained by the pendant drop method, and by inverting bubble nucleation experiments using the Classical Nucleation Theory (CNT). Our model reproduces the systematic decrease in surface tension with increased H2O pressure observed in the experiments. On the other hand, the effect of temperature is confirmed by the experiments only at high pressure. At atmospheric pressure, the model shows a decrease of surface tension with temperature. This is in contrast with a number of experimental observations and could be related to microstructural effects that cannot be reproduced by our model. Finally, our analysis indicates that the surface tension measured inverting the CNT may be lower than the value measured by the pendant drop method, most likely because of changes in surface tension controlled by the supersaturation. Ó 2015 Elsevier Ltd. All rights reserved.

1. INTRODUCTION Magmatic vesiculation is a complex process involving bubble nucleation (Toramaru, 1989; Navon and Lyakhovsky, 1998; Gonnermann and Gardner, 2013), growth (Proussevitch and Sahagian, 1998; Gardner et al., 2000; Lensky et al., 2002; Masotta et al., 2014), ripening (Lautze et al., 2011) and coalescence (Sahagian et al., 1989; Castro et al., 2012). The supersaturation pressure triggering the nucleation of gas bubbles in a magmatic melt hinges on the melt–gas surface tension (Navon and

⇑ Corresponding author.

E-mail address: [email protected] (S. Colucci). http://dx.doi.org/10.1016/j.gca.2015.10.037 0016-7037/Ó 2015 Elsevier Ltd. All rights reserved.

Lyakhovsky, 1998). For many magmas the dominant volatile is H2O and we focus on H2O to develop a model of melt–gas surface tension. Furthermore, we know from rheology theories for soft inclusion suspensions that the Capillary number (i.e., the ratio of viscous stress and surface tension) characterizes the ability of bubbles to deform in a shear flow. So, knowledge of the surface tension provides input for the estimation of the time-dependent rheology of bubble-bearing magmas as well (Mader et al., 2013). The experimental method of the pendant drop offers an attractive means of determining surface tensions for polymer–gas and polymer–polymer combinations. A full discussion of the method is available in Couper (1993). This technique, widely used in industrial applications (Harrison et al., 1996; Jaeger et al., 2002; Park et al., 2007), is not a

114

S. Colucci et al. / Geochimica et Cosmochimica Acta 175 (2016) 113–127

standard method in the study of surface tension in magmatic melts. As far as we know, only three studies (Epelbaum et al., 1973; Khitarov et al., 1979; Bagdassarov et al., 2000) measured the high-pressure/ high-temperature surface tension of hydrous magmas employing the pendant drop method, with somewhat different results. For example, although a systematic decrease in the surface tension with increased H2O pressure is always observed, Epelbaum et al. (1973) found smaller values compared to Bagdassarov et al. (2000) for similar compositions. Recent experimental studies (Mangan and Sisson, 2000, 2005; Gardner and Ketcham, 2011; Gardner, 2012; Gardner et al., 2013) evaluated the surface tension by determining first, through decompression experiments, the critical supersaturation pressure for bubble nucleation as a function of temperature and dissolved water content, and then inverting for the surface tension using Classical Nucleation Theory (CNT). Gonnermann and Gardner (2013) found, by integrating homogeneous bubble nucleation experiments with numerical modelling based on a nonClassical Nucleation Theory, that the surface tension between critical bubble nuclei and the surrounding melt depends on the degree of supersaturation. They proposed a thermodynamical formulation where the non-classical surface tension, rnc , is defined as a function of r1 , the surface tension value at equilibrium condition, and a parameter, n, determined by fitting homogeneous bubble nucleation experiments rnc ¼ r1 ð1  n2 Þ

1=3

:

ð1Þ

This model was verified with a hydrous rhyolitic melt, using as input for r1 the values measured by Bagdassarov et al. (2000) by the pendant drop method. Results of the work by Gonnermann and Gardner (2013) imply that the surface tension value inferred by inverting the CNT may be lower than the value measured by the pendant drop method. This is consistent with the observation that far from equilibrium the interface between a nucleus and a surrounding metastable bulk phase is not sharp but diffuse (Gonnermann and Gardner, 2013). Despite several experimental studies, a physics based model of the surface tension of magmatic melts is still lacking. By definition, the surface tension is the difference (per unit area of interface) between the actual free energy of the non-uniform system and the free energy the system would have if the properties of the phases were continuous throughout the interface. The gradient theory (Cahn and Hilliard, 1958), extended by Poser and Sanchez (1981) to binary mixtures, can be used for modelling the interfacial properties of a planar interface between a silicate melt and a gas phase (here H2O gas, but the model might be extended also to other volatile species, for example CO2). This approach has been successfully applied in several studies of industrial liquid mixtures in contact with supercritical gases (Poser and Sanchez, 1981; Harrison et al., 1996; Kahl and Enders, 2000; Enders et al., 2005) but has never been used before to study magmatic compositions. In the present work, we apply the gradient theory to study the equilibrium surface tension (r1 in Eq. (1)) of hydrous silicate melts in contact with supercritical H2O at

high temperature and low to moderate pressures. Even if magmatic melts are not exclusively silicate, as for example carbonate melts, and other volatile components, as for example CO2, can be present in magmatic systems, in this first application of the model we only consider silicate melts in contact with H2O gas. We calibrate and verify the model against experimental data on different melt compositions obtained by the pendant drop method and by inverting bubble nucleation experiments using the CNT. Our thermodynamical model successfully reproduces the decrease of surface tension with increase of H2O pressure observed in experiments. On the other hand, the effect of temperature on the surface tension is confirmed by the experiments only at high pressure, maybe because of microstructural effects that cannot be reproduced by our model. Furthermore, our analysis indicates that the experimental method for the calculation of the surface tension by inverting bubble nucleation experiments using the CNT may differ by pendant drop method experiments, because of changes in the surface tension controlled by the supersaturation as suggested by Gonnermann and Gardner (2013). 2. METHODS 2.1. Thermodynamical modelling Using a procedure similar to that proposed by Cahn and Hilliard (1958) and Poser and Sanchez (1981) developed a thermodynamical model of the surface tension of binary mixtures by considering the change in the partial densities of the two mixture components within a planar interface. The surface tension equation between the hydrous melt and H2O gas can be written as a function of temperature and pressure (Poser and Sanchez, 1981; Enders et al., 2005) " # Z mol
melt;l q pffiffiffiffiffiffiffiffiffi pffiffiffi pffiffiffiffiffiffiffiffiffiffi D
qmol H2 O rðP ; T Þ ¼ 2 M melt k melt þ M H2 O k H2 O mol
melt;l q 0 pffiffiffiffiffiffi mol  Dad
qmelt ; ð2Þ where M melt and k melt are the magmatic melt molecular weight and influence parameter, M H2 O and k H2 O are the H2O molecular weight and influence parameter, D
qmol H2 O is the difference between the H2O partial volume molar densi
mol ties of the two phases (gas and liquid) at equilibrium, q melt;l is the partial volume molar density of the melt at equilibrium, Da is the Grand Thermodynamic Potential of the binary mixture, that takes into account the excess of Free
mol Energy across the liquid–gas interface, and q melt is the partial volume molar density of the melt across the interface (further details in Step 5). The application of (2) to calculate the surface tension of a binary mixture of magmatic melt and H2O gas requires five steps: (1) determine the parameters of the equation of state (EOS) for the two pure components (i.e., anhydrous silicate melt or H2O); (2) determine the influence parameter of pure H2O; (3) determine the influence parameter of the pure melt; (4) determine the parameters of the EOS of the binary mixture (i.e., anhydrous magmatic melt and dis-

S. Colucci et al. / Geochimica et Cosmochimica Acta 175 (2016) 113–127

solved H2O); and (5) calculate the surface tension at the interface between the binary mixture and the H2O gas. 2.1.1. Step 1: EOS parameters for pure components We estimate the parameters of the EOS of the pure components (i.e., anhydrous silicate melt or H2O) by fitting the appropriate EOS to experimental data of density as a function of pressure and temperature. One EOS suitable for this purpose is the Sanchez–Lacombe equation of state (SLEOS; Sanchez and Lacombe, 1976; Poser and Sanchez, 1981; Harrison et al., 1996; Kahl and Enders, 2000; Enders et al., 2005). The SL-EOS is a lattice fluid equation (Sanchez and Lacombe, 1976) calibrated on experimental data    1 ~ ~ 2 2 ~ 2 ~Þ þ 1  q þ P þ T lnð1  q ¼ 0; ð3Þ r ~; Pe ; Te are the reduced density, pressure and temwhere q perature, and r is the number of lattice sites occupied by single H2O or melt molecules. The reduced density, pressure and temperature are defined by the ratio between the actual _ P_ ; T_ density, pressure, temperature and the parameters q; ~¼ q

q q_

P Pe ¼ P_

T Te ¼ : T_

v_ ¼

RT_ P_

r¼

M P_ ; RT_ q_

ð5Þ

where R is the gas constant and M is the molecular weight. _ represents the interaction energy term, v_ the characteristic volume and r the number of lattice sites occupied by single H2O or melt molecules (see also Tables 1 and 2 for notation). Machida et al. (2010) proposed, for polar liquids (e.g., water), an alternative formulation that takes into account the temperature dependence of hydrogen bonding and ionic interactions in the interaction energy term _ _ ðT Þ ¼ 0

aT ; 1 þ aT

ð6Þ

where the parameter a controls the temperature dependence and 0 is the asymptotic value of the interaction energy. The SL-EOS can be calibrated by either fitting the stan_ P_ ; T_ – or, dard equation (Eqs. (4) and (5)) – parameters q; Table 1 Notation. k k melt;H2 O M P R r T v_ w

Influence parameter Interaction parameter molecular weight Pressure Gas constant Number of lattice sites occupied by a molecule Temperature Characteristic volume Weight fraction

Table 2 Notation. Greeks a Da _ 0 l q qmol r /

Temperature dependence term Grand Thermodynamic Potential Interaction energy term Asymptotic value of the interaction energy Chemical potential Density Molar density Surface tension Volume fraction

K1 J m3 J mol1 J mol1 J mol1 kg m3 mol m3 J m2 dimless

Subscripts l g eq h mix nc

Liquid phase Gas phase Equilibrium Homogeneous Mixture Non-classical

Accent marks  . –

Reduced Parameter Partial

ð4Þ

_ P_ and T_ are related to each other The three parameters q; by the following relations _ ¼ RT_

115

J m5 kg2 dimless kg/mol Pa J mol1 K1 dimless K m3 mol1 dimless

alternatively, the modified formulation (Eqs. (4)–(6)) – parameters v_ ; r 0 ; a – to experimental data of density as a function of temperature and pressure. _ P_ and T_ of the We determine, first, the parameters q; SL-EOS of the pure, anhydrous silicate melt by fitting the standard equation (Eqs. (4) and (5)) to the model density curve by Papale et al. (2006) (online calculator from http://ctserver.ofm-research.org/Papale/Papale.php). This model curve is based on a very large database of silicate liquid systems. To improve the accuracy of the fit between the SL-EOS and the model curve, we assume that q_ and T_ are linear functions of the temperature (see formulas in Table 3). Fig. 1 shows a comparison between the densities determined by the SL-EOS and those from the model curve by Papale et al. (2006) for the melt composition discussed in this work. The results (Fig. 1) show a maximum deviation from the model curve by Papale et al. (2006) between 0.3% and 0.4%. Finally, we estimate the parameters 0 ; a; v_ and r of the SL-EOS of pure H2O by fitting the modified formulation (Eqs. (4)–(6)) to experimental data from Vukalovich et al. (1961, 1962), Hilbert et al. (1981), NISTIR (1998) and IAPWS (2007) (see Table 3, last two lines). The results shown in Fig. 2 indicate a maximum deviation of 2:5% between the SL-EOS and the experimental data for H2O. 2.1.2. Step 2: Influence parameter of H2O (k H2 O ) Although influence parameters have a molecular theoretical definition, a more practical approach is to estimate them by fitting the surface tension model for pure components to experimental data (Cahn and Hilliard, 1958; Poser and Sanchez, 1979; Kahl and Enders, 2000). The equation for the surface tension r of a liquid water–H2O gas system (Cahn and Hilliard, 1958), is given, in dimensional form, by

116

S. Colucci et al. / Geochimica et Cosmochimica Acta 175 (2016) 113–127

Table 3 Pure component parameters of melts (P_ ; q_ ¼ aq_ T þ bq_ ; T_ ¼ aT_ T þ bT_ ) and H2O (a; 0 ; v; r). Reference chemical compositions: (1) Bagdassarov et al. (2000); (2) Khitarov et al. (1979); (3) Gardner and Ketcham (2011); (4) Mangan and Sisson (2005); (5) Gardner (2012); (6) Iacono Marziano et al. (2007); (7) Gardner and Ketcham (2011); (8) Colucci (2012). Melt composition 1

HPG8 Basalt2 Dacite3 Dacite4 Na-Phonolite5 K-Phonolite6 Rhyolite7 Trachyte8 H2O

r ¼ 2M H2 O

pffiffiffiffiffiffiffiffiffiffi k H2 O

Z

qmol l

qmol g

P_ [MPa]

aq_ [kg/(m3T)]

bq_ [kg/(m3T)]

4900 4900 4900 4900 4900 4900 4900 4900

0:11859 0:22456 0:15372 0:15 0:20161 0:22163 0:12564 0:19744

2757:9 3225:6 2903:8 2908 2960:5 3050:3 2790 2987:7

2:202 2:5592 2:2915 2:2762 2:2045 2:1921 2:2088 2:2011

a [1/K]

0 [J mol1]

v_ [c3/mol]

r [dimless]

0:01

7026:12484

3:36314

4:77865

pffiffiffiffiffiffi mol Dadq ;

ð7Þ

where M H2 O is the H2O molecular weight, k H2 O the H2O influence parameter, qmol and qmol are the molar densities g l at equilibrium in the coexisting liquid and gas phases and Da is the Grand Thermodynamical potential. The molar and qmol are given by densities qmol g l ¼ qmol l

q~l q_ ; M H2 O

qmol ¼ g

q~g q_ ; M H2 O

ð8Þ

where q_ is the parameter for pure H2O calculated in Step 1 and q~l and q~g are the reduced densities of the coexisting H2O liquid and H2O gas phases at the equilibrium pressure P eq and given temperature T from the modified SL-EOS (3– 6). The equilibrium pressure P eq must satisfy the following criterion ll ðP eq Þ ¼ lg ðP eq Þ;

ð9Þ

where l represents the chemical potential. The chemical potential for a pure liquid or gas component (liquid and gas H2O in this case) is obtained by the Gibbs free energy of the SL-EOS and is given by Poser and Sanchez (1979) " # ~Þlogð1  q ~Þ 1 ~ Pe ð1  q q ~Þ ; ð10Þ l ¼ rRT  þ þ þ logðq ~ r q ~ Te Te q ~ ; Pe and Te are the reduced density, pressure and where q temperature from the SL-EOS (3) and r is given by Eq. (5). The Grand Thermodynamical potential Da takes into account the excess of Free Energy across the liquid–gas interface. According to the Free Energy theory of nonuniform systems (Cahn and Hilliard, 1958), the local free energy per molecule in a region of non-uniform composition (i.e., the interface between liquid and gas) will depend on both the local composition and the composition of the immediate environment, expressed by the local composition derivatives. The Grand Thermodynamical potential is given by the difference between the local composition term of the

aT_ [dimless]

bT_ [dimless] 516:76 687:43 559:57 555:41 533:4 534:36 525:66 533:56

free energy qlh (where the subscript h stands for homogeneous) and the free energy qleq , calculated assuming equilibrium conditions between the two phases Da ¼ qmol ½lh  leq :

ð11Þ

The equilibrium chemical potential, leq , is calculated imposing the equilibrium condition (9), while the homogeneous chemical potential lh is given by (10) at given T and P eq , and qmol is the integration variable in (7). To estimate the influence parameter k H2 O , we fit the model surface tension of H2O (7) to the experimental data from Vargaftik et al. (1983) (Table 4, Fig. 3). 2.1.3. Step 3: Influence parameter of the magmatic melt (k melt ) Because of the high vapour tension of silicate melts, the surface tension of pure melt is hardly measurable. For this reason, we cannot employ the same method to estimate k melt as used to find k H2 O . For low molecular weight polymers, the influence parameter k is a function of the molecular weight and approaches a constant value for increasing values of the molecular weight (Enders et al., 2005). Using this observation, we can estimate k melt by fitting the equation for the surface tension between the hydrous melt and H2O gas (Eq. (2); see next Step 5) to the surface tension measured by the pendant drop experiments of Khitarov et al. (1979) and Bagdassarov et al. (2000) (Fig. 4) in order to find a relation with the molecular weight (see Table 4). This approach is consistent with our initial assumption of equilibrium conditions, since experiments with the pendant drop method provide values of the surface tension at equilibrium (Gonnermann and Gardner, 2013). 2.1.4. Step 4: EOS parameters of the binary mixture The EOS parameters of the two pure components can be employed to calculate the parameters of the EOS of the binary mixture using the following mixing rules (Sato et al., 1996; Enders et al., 2005)

S. Colucci et al. / Geochimica et Cosmochimica Acta 175 (2016) 113–127

117

Fig. 1. Melt density as a function of pressure at different temperatures from the SL-EOS proposed in this work (solid lines) compared against the calibration curves (dashed lines) from the model by Papale et al. (2006). The legend in the upper, left panel is valid for all the other panels as well. Chemical compositions used as input for the models are from (1) Bagdassarov et al. (2000), (2) Khitarov et al. (1979), (3A) Gardner (2012), (3B) Iacono Marziano et al. (2007), (4A, 5) Gardner and Ketcham (2011), (4B) Mangan and Sisson (2005), (6) Colucci (2012).

118

S. Colucci et al. / Geochimica et Cosmochimica Acta 175 (2016) 113–127

Fig. 2. H2O gas density as a function of pressure at different temperatures reproduced by the SL-EOS (solid lines). Symbols refer to the experimental data used for calibration (see Step 1 for references).

P_ ¼

XX /i /j P_ij ; i

ð12Þ

j

qffiffiffiffiffiffiffiffiffi P_ij ¼ ð1  k ij Þ P_ i P_j ; X /i;0 T_ i ; T_ ¼ P_ P_ i i 1 X /i;0 ¼ ; r ri i / P_ i =T_ i /i;0 ¼ P i ; _ _ j /j P j =T j wi =qi ; /i ¼ P j wj =qj X v_ ¼ /i v_i ;

ð13Þ ð14Þ ð15Þ ð16Þ

H2O solubility from the model of Papale et al. (2006). We calculate the H2O solubility for each magmatic melt composition used to verify our model (Fig. 5). The model solubility, at a given P and T, is given by the /H2 O that satisfies the equilibrium criterion (9) applied to a binary system lH2 O;l ¼ lH2 O;g ;

ð19Þ

where the chemical potential of H2O is calculated in the liquid binary mixture (lH2 O;l ) and in the pure gas phase (lH2 O;g ). The equation for the chemical potential of H2O in the binary mixture is given by Sato et al. (1996).

ð17Þ ð18Þ

i

where the indexes i; j loop over the two components (magmatic melt and H2O), / is the volume fraction, w the weight fraction and k ij the binary interaction parameter. The only k ij that is different from 0 is k melt;H2 O ¼ k H2 O;melt . In the next equations, we indicate the EOS parameters of the pure components with the indexes melt and H2O (e.g., P_ H2 O and P_ melt ; T_ H2 O and T_ melt ; q_ H2 O and q_ melt ), and the parameters of the binary mixture without any index _ The volume fraction / (or w, the weight (e.g., P_ ; T_ ; q). fraction), pressure P and temperature T are independent variables of the SL-EOS (3) of the binary mixture. The binary interaction parameter k melt;H2 O is calculated by fitting the model solubility of H2O to the curves of

Table 4 Influence parameters of melt and H2 O : k melt parameters obtained by calibrating the model for the surface tension of a binary mixture with the pendant drop experiments of (1) Bagdassarov et al. (2000) and (2) Khitarov et al. (1979) in the temperature interval 1173– 1873 K (results of calibration in Fig. 4). The resulting linear relation used to find k melt for other molecular weight compositions is given by: k melt ¼ 8:3333  1015 –5:025  1016 ; k H2 O parameter obtained by calibrating the model for the surface tension of pure H2O to the experimental data of Vargaftik et al. (1983) (results of calibration in Fig. 3). Melt composition

k melt [J m5 kg2]

HPG81 basalt2

3:5  1017 2:5  1017 k H2 O [J m5 kg2]

H2O

2:5531  1017

S. Colucci et al. / Geochimica et Cosmochimica Acta 175 (2016) 113–127

119

Fig. 3. H2O surface tension calculated by the model using the Gradient theory combined with the standard SL-EOS formulation (dashed line) and the modified SL-EOS formulation (solid line). Symbols refer to the experimental data of Vargaftik et al. (1983) used for calibration. The modified SL-EOS formulation, employed in this work, shows a lower absolute average deviation (AAD) compared to the standard formulation. ( )   P_ melt þ P_ H2 O  2 Pe melt;H2 O 2 rH O ~ lH2 O;l ¼ RT ln/H2 O þ 1  2 /melt þ rH2 O q /melt rmelt P_ H2 O Te H2 O (  ) e ~ ~ P H2 O 1 q q ~Þlnð1  q ~Þ þ þ rH2 O RT  þ þ ð1  q ln~ q ; ~ r H2 O ~ Te H2 O q Te H2 O q ð20Þ

~ is the reduced density calculated from the SL-EOS where q (Eq. (3)) for a binary mixture with the parameters defined by the mixing rules (12)–(18). The above Eq. (20) reduces to (10) – i.e., the chemical potential of H2O in the pure gas phase (lH2 O;g ) – for /melt ¼ 0. The solubility values used for fitting are calculated from the model curves by Papale et al. (2006) – see Fig. 5. To improve the model best fit to the experimental data, we assumed that the binary interaction parameter k melt;H2 O is a polynomial function of pressure and temperature (see formulas in Table 5). In Fig. 5, the density of the binary mixture obtained by the SL-EOS is compared with the density from the model curves by Papale et al. (2006). Our model predicts lower solubility at low pressure (i.e., 6100 MPa, Fig. 5) – or an increase of the solubility with increasing pressure. The largest deviation from the model curves by Papale et al. (2006) is around 5% (Fig. 5). 2.1.5. Step 5: Surface tension of the binary mixture We can now employ (2) to calculate the surface tension of the binary mixture, hydrous melt and H2O gas. The partial volume molar density of the melt at equilibrium is given by
melt;l ¼ /melt qmol ¼ ð1  /H2 O Þqmol ; q

ð21Þ

where /H2 O (i.e., solubility) satisfies the equilibrium criterion (19) and qmol is calculated by the EOS of the binary

mixture for the same P, T and /H2 O . D
qmol H2 O is the difference between the H2O partial molar densities at equilibrium in the two phases mol
mol
mol D
qmol  qmol H2 O ¼ q H2 O;l  q H2 O;g ¼ /H2 O q H2 O ;

ð22Þ

where qmol H2 O is given by the EOS of pure H2O for a given P and T. The Grand Thermodynamic Potential of the binary mixture is
mol
mol Da ¼ q H2 O ðlH2 O;h  lH2 O;eq Þ þ q melt ðlmelt;h  lmelt;eq Þ;

ð23Þ


mol
mol where q melt and q H2 O are the partial molar densities of the melt and H2O across the interface, leq is the equilibrium chemical potential (see Eqs. (19) and (20)) and lh the actual homogeneous chemical potential across the interface calculated by (20) for a given T and P and the partial densities calculated across the interface. Calculations across the interface assume a linear partial density profile across the interface (Harrison et al., 1996; Enders et al., 2005). The
mol partial density of the melt q melt is a linear function of the
mol partial density of H2O q H2 O
melt ð

melt;l  q qH2 O Þ ¼ q


melt;l q
H2 O;l Þ: ð
qH2 O  q D
qH2 O

ð24Þ


melt , we can In this way when we solve (2), integrating for q
H2 O from (24) and we can calculate q ~ and /H2 O obtain q from
melt ¼ ð1  /H2 O Þ q

~q_ q ; M mix

~q_ q ; M mix 1 wmelt wH2 O ¼ þ : M mix M melt M H2 O


H2 O ¼ /H2 O q

ð25Þ ð26Þ ð27Þ

120

S. Colucci et al. / Geochimica et Cosmochimica Acta 175 (2016) 113–127

Fig. 4. Surface tension of the binary mixture melt–H2O obtained through the calibration of the model with the pendant drop experiments of Bagdassarov et al. (2000) (upper panel-haplogranite HPG8) and Khitarov et al. (1979) (lower panel-basalt) in the temperature interval 1173– 1873 K. The errors of the experimental data are about the size of the symbols.

~ and /H2 O to calculate the homoFinally, we can employ q geneous chemical potential across the interface (20) and the Grand Thermodynamical Potential (23). 3. RESULTS Model results predict a surface tension decrease with increasing H2O pressure (see Figs. 4 and 6) as has been observed in several experiments (Epelbaum et al., 1973; Khitarov et al., 1979; Bagdassarov et al., 2000; Mangan and Sisson, 2000, 2005; Mourtada-Bonnefoi and Laporte,

2004; Iacono Marziano et al., 2007; Gardner and Ketcham, 2011; Colucci, 2012; Gardner, 2012). The model has been calibrated with the pendant drop experiments of Bagdassarov et al. (2000) (haplogranite) and Khitarov et al. (1979) (basalt) in the temperature interval 1173–1873 K. In Fig. 4 we show the results of the calibration of the model. For the haplogranite at 1273 K (see Fig. 4, upper panel), the model is able to reproduce the changes in surface tension with pressure with a maximum error of 13%. The increase in surface tension with temperature is well reproduced by our model at high pressure (300 MPa) but the model shows a

S. Colucci et al. / Geochimica et Cosmochimica Acta 175 (2016) 113–127

121

Fig. 5. Model solubility (upper panel) and mixture density (lower panel) obtained by the SL-EOS for the binary mixture (solid lines) compared to the results obtained by the model of Papale et al. (2006) used for calibration (dashed lines). Colours are the same as in Fig. 1. Chemical compositions used as input for the models are from (1) Bagdassarov et al. (2000), (2) Khitarov et al. (1979), (3-thin lines) Gardner (2012), (3-thick lines) Iacono Marziano et al. (2007), (4-thin lines, 5) Gardner and Ketcham (2011), (4-thick lines) Mangan and Sisson (2005), (6) Colucci (2012).

decrease of surface tension with temperature at atmospheric pressure, a behaviour opposite to that recorded in the experiments. The basalt experiments of Khitarov et al. (1979) at 1473 K (see Fig. 4, lower panel) are reproduced by the model within the experimental error. Again, the model shows a negative temperature dependence of the surface tension at atmospheric pressure and a positive dependence at high pressure.

The two k melt parameters calculated by calibrating the model with the drop method experiments (Table 4), have been linearly interpolated to find the k melt parameters for other compositions, on the base of the molecular weight (see equation in caption of Table 4). In Fig. 6, we compare our model results for surface tension in the temperature range 1173–1873 K against experimental data on rhyolite (Mourtada-Bonnefoi and Laporte, 1999, 2002, 2004;

122

S. Colucci et al. / Geochimica et Cosmochimica Acta 175 (2016) 113–127

Fig 5. (continued)

Table 5 Interaction parameters for the binary mixture melt–H2O. Parameters depend on pressure and temperature through the relation k melt;H2 O ¼ aP 2 þ bP þ c , a ¼ a1 T þ a2 ; b ¼ b1 T þ b2 ; c ¼ c1 T þ c2 . Reference chemical compositions: (1) Bagdassarov et al. (2000); (2) Khitarov et al. (1979); (3) Gardner and Ketcham (2011); (4) Mangan and Sisson (2005); (5) Gardner (2012); (6) Iacono Marziano et al. (2007); (7) Gardner and Ketcham (2011); (8) Colucci (2012). Melt composition

a1 [K1 MPa2]

a2 [MPa2]

b1 [K1 MPa1]

b2 [MPa1]

HPG81 Basalt2 Dacite3 Dacite4 Na–Phonolite5 K–Phonolite6 Rhyolite7 Trachyte8

4:0857  1010 4:5857  1010 4:0057  1010 4:0  1010 6:65  1010 6:65  1010 4:1857  1010 4:2857  1010

9:6292  107 9:6292  107 9:6292  107 9:6292  107 9:6292  107 9:6292  107 9:6292  107 9:5292  107

1:9071  107 2:0  107 1:9271  107 1:9271  107 0:6571  107 0:6571  107 1:8571  107 1:8571  107

0:00017681 0:00017681 0:00017681 0:00017681 0:00017681 0:00017681 0:00017681 0:00017681

c1 [dimless]

c2 [dimless]

9:5714  105 8:5714  105 9:0014  105 9:6014  105 8:60  105 8:45  105 9:3714  105 9:9714  105

0:052177 0:052177 0:052177 0:052177 0:052177 0:052177 0:052177 0:052177

Mangan and Sisson, 2000; Gardner and Ketcham, 2011), phonolite (Iacono Marziano et al., 2007; Gardner, 2012), dacite (Mangan and Sisson, 2005; Gardner and Ketcham, 2011) and trachyte (Gardner et al., 2013; Colucci, 2012) obtained by inverting bubble nucleation experiments using the CNT. The equilibrium surface tension predicted by our model is higher than the values measured by this method (Fig. 6). Only two experiments with rhyolites (Fig. 6b,

upper panel) show a surface tension larger than surface tension predicted by our thermodynamical model. Furthermore, the equilibrium surface tension predicted by our model for Na- and K-phonolite (Fig. 6a, upper panel) and trachyte (Fig. 6b, lower panel) is slightly higher (ca. 20%) than the surface tension predicted for basalt (Fig. 4, lower panel), rhyolite (Fig. 6b, upper panel) and dacite (Fig. 6a, lower panel).

S. Colucci et al. / Geochimica et Cosmochimica Acta 175 (2016) 113–127

123

Fig. 6. Model results for the surface tension of the binary mixture melt–H2O in the temperature range 1173–1873 K against literature experimental data obtained by inverting bubble nucleation experiments using the Classical Nucleation Theory. The errors of the experimental data are about the size of the symbols. The arrow indicates that the value of surface tension is a minimum. The equilibrium surface tension predicted by our model is higher than the experimental values (see the text for discussions).

4. DISCUSSIONS The dependence of surface tension on pressure and temperature can be explained by our model in terms of thermodynamics. Surface tension (Eq. (2)) is the product of a term positively dependent on the density contrast between the hydrous melt and the H2O gas phase (in square brackets before the integral in Eq. (2) – function C in the upper panel of Fig. 7) that multiplies an integral term taking into account the excess of Free Energy across the liquid–gas interface – function F in the upper panel of Fig. 7. In

Fig. (7), we show the contribution from these terms on the surface tension. When the H2O pressure increases both these terms decrease (see Fig. 7, upper panel), hence the surface tension decreases (Fig. 7, lower panel). In fact, when the H2O pressure increases, the density of the H2O increases (Fig. 2) and the density of the hydrous melt decreases (Fig. 5, lower panels) because of the dominant effect of the increased solubility of H2O (Fig. 5, upper panels), reducing the density contrast between the hydrous melt and the H2O gas phase. An increase in temperature has the opposite effect (Fig. 7 upper panel): the density contrast

124

S. Colucci et al. / Geochimica et Cosmochimica Acta 175 (2016) 113–127

Fig 6. (continued)

between the hydrous melt and the H2O gas phase (function C) increases, while the excess of Free Energy (function F) across the interface decreases. At low pressure (Fig. 7, upper panel) the lower effect of temperature on the H2O gas density (Fig. 2) makes the decrease of the excess of Free Energy (function F) larger than the increase of the term accounting for the density contrast, (function C), forcing a decrease of surface tension with temperature (Fig. 7, lower panel). At high pressure (Fig. 7, upper panel), when the H2O gas density is more sensitive to temperature (Fig. 2), the increase in C dominates and the surface tension increases (Fig. 7, lower panel).

Comparing model results for haplogranite and basalt with the pendant drop experiments used for the calibration (Fig. 4), we observe that, while the trend of variation of surface tension with pressure is well reproduced by our model, the effect of temperature is validated by the experiments only at high pressure. In fact, the model shows a decrease of surface tension with temperature at atmospheric pressure, a behaviour opposed to that measured in the experiments by Bagdassarov et al. (2000) – hydrous haplogranitic melts – and by Walker and Mullins (1981) – anhydrous basaltic melts. The increase of surface tension with increasing temperature at low pressure, registered in

S. Colucci et al. / Geochimica et Cosmochimica Acta 175 (2016) 113–127

125

Fig. 7. Surface tension (upper panel) is the product of a term positively dependent on density contrast between the hydrous melt and the H2O gas phase, C, and a term that takes in account the excess of Free Energy across the liquid–gas interface, F (lower panel). By increasing pressure both of these terms decrease decreasing the surface tension. The increase of temperature has an opposite effect, it increases the density contrast between the hydrous melt and the H2O gas phase, C, and decreases the excess of Free Energy, F, across the interface. At low pressure F is predominant on C, determining a decrease of surface tension with temperature. At high pressure the increase of C dominates and surface tension increases.

the experiments on silicate melts, is uncommon, since the surface tension of pure liquids and binary mixtures of polymer melts and gas typically decreases with increasing temperatures (e.g., see H2O in Fig. 3, Park et al. (2007)). For example, in the experiments reported in Park et al. (2007) the dependence of surface tension on temperature, for polystyrene in carbon dioxide, is negative and becomes slightly positive with increasing pressure (from 0.1 to ca. 17 MPa). Regarding the unusual positive temperature effect produced by our model, Bagdassarov et al. (2000) suggest that the positive temperature dependence might be explained either by the dissolution of gas molecules in the surface layer or by the differentiation of the free drop surface into sub-layers with different chemical compositions. Walker and Mullins (1981) conclude that the increase of surface tension with

temperature, at atmospheric pressure, is determined by structural dissociation that reduces the molar volume of the average structural unit by at most 5% for a 373 K increase in temperature. Finally, the lack of a comprehensive experimental dataset of the surface tension at equilibrium (r1 in Gonnermann and Gardner (2013)) for hydrous magmas represents the major drawback for a rigorous calibration of the model and minimization of the mismatch between model end experiments. In this work, the k melt parameter of the model has been calibrated with the few available experiments on equilibrium surface tension for basalt and haplogranite. Nevertheless, since the influence parameter k is a smooth function of the molecular weight across several order of magnitude for a number of industrial polymers (Enders et al., 2005) and

126

S. Colucci et al. / Geochimica et Cosmochimica Acta 175 (2016) 113–127

the molecular weights of the studied silicate melts fall in a very limited interval (63.3–67.7 g/mol), we believe that the linear interpolation of the k melt adopted here to extrapolate the value of the parameter to other compositions is adequate, at least as a first order approximation. When we compare the equilibrium surface tension calculated with our model, equivalent to r1 in Gonnermann and Gardner (2013), against the values obtained by inverting bubble nucleation experiments using the CNT, we can see that the values of the equilibrium surface tension predicted by our model are higher than those estimated using the CNT (Fig. 6). This result is consistent with the observation by Gonnermann and Gardner (2013) that the surface tension between critical bubble nuclei and the surrounding melt depends on the degree of supersaturation because – far from equilibrium – the interface between a nucleus and a surrounding metastable bulk phase is not sharp but diffuse. The slightly higher equilibrium surface tension predicted by our model for Na- and K-phonolite (Fig. 6a, upper panel) and trachyte (Fig. 6b, lower panel) with respect to the surface tension predicted for basalt (Fig. 4, lower panel), rhyolite (Fig. 6b, lower panel) and dacite (Fig. 6a, lower panel) is explained by the higher molecular weight of phonolites and trachyte (respectively, 67.2 g/mol and 67.7 g/mol compared to basalt – 63.3 g/mol – rhyolite – 64.8 g/mol – and dacite – 65.2 g/mol). This is in agreement with the conclusion by Gardner et al. (2013) that the composition of silicate melts should have a limited impact on the surface tension (ca. 20% in the non-equilibrium surface tension calculated by Gardner et al. (2013)). Two experiments with rhyolites (Fig. 6b, upper panel) show a surface tension larger than surface tension predicted by our thermodynamical model. Following Gardner and Ketcham (2011), these relatively cold rhyolites could have been viscous enough (viscosity > 104 Pa s) to partially delay the bubble nucleation, giving a biased higher value of the surface tension. Finally, the equilibrium surface tension (r1 ) provided by the model could be used to infer the non-classical surface tension (rnc ) once determined the experimental parameter n as suggested by Gonnermann and Gardner (2013) (Eq. (1)). Since disequilibrium degassing may be an integral component of explosive volcanic eruptions (Mangan and Sisson, 2000), in some cases, as for example for high-viscous magmas rising up into the volcanic conduit, non-classical surface tension is required to take in account the effect of supersaturation. 5. SUMMARY AND CONCLUSIONS We employed the gradient theory to develop a thermodynamical model of the equilibrium surface tension (r1 in Gonnermann and Gardner (2013)) of H2O gas in contact with hydrous silicate melts. The model has been calibrated against the pendant drop experiments by Bagdassarov et al. (2000) (haplogranite) and Khitarov et al. (1979) (basalt) in the temperature interval 1173–1873 K. The changes in surface tension with pressure are reproduced by our model with a maximum error of 13%. The effect of temperature

is confirmed by the experiments only at high pressure. At atmospheric pressure, the model shows a decrease of surface tension with temperature. This is in contrast with the experimental observations on hydrous haplogranitic melts by Bagdassarov et al. (2000) and those on anhydrous basaltic melts by Walker and Mullins (1981). The increase in surface tension with increasing temperature at low pressure, observed in the experiments on silicate melts, is decidedly uncommon if compared to pure liquids and polymer–gas mixture and may be related to microstructural effects that cannot be reproduced by our model. Since for different industrial polymers the influence parameter k melt is a smooth function of the molecular weight across several orders of magnitude and the molecular weights of the studied silicate melts fall in a limited interval (63.3–67.7 g/mol), k melt can be linearly interpolated to estimate the parameters for other compositions. Model results in the temperature range 1173–1873 K have been compared against literature experimental data on rhyolite, phonolite, dacite and trachyte obtained by inverting bubble nucleation experiments using the Classical Nucleation Theory (CNT). The equilibrium surface tension predicted by our model is higher than the values measured by the method of inverting the CNT. This is consistent with the observation that the surface tension between critical bubble nuclei and the surrounding melt depends on the degree of supersaturation because far from equilibrium the interface between a nucleus and a surrounding metastable bulk phase is diffuse instead of sharp (Gonnermann and Gardner, 2013). ACKNOWLEDGMENTS We would like to thank T. Sisson for the helpful suggestions and the constructively critical review of this paper. The authors are grateful to J. Gardner and two anonymous reviewers for their careful review and useful and constructive comments.

REFERENCES Bagdassarov N. S., Dorfman N. A. and Dingwell D. B. (2000) Effect of alkalis, phosphorus, and water on the surface tension of haplogranite melt. Am. Mineral. 85, 33–40. Cahn J. W. and Hilliard J. E. (1958) Free energy of a nonuniform system. I. Interfacial free energy. J. Chem. Phys. 28, 258–267. Castro J. M., Burgisser A., Schipper C. I. and Mancini S. (2012) Mechanisms of bubble coalescence in silicic magmas. Bull. Volcanol. 74, 2339–2352. Colucci, S. (2012) Modelling of bubble nucleation in trachyphonolitic magmas: implications for the dynamics of ash-rich eruptions. Ph. D thesis. Sapienza-Universita` di Roma (Italy). pp. 35–37, available at: http://hdl.handle.net/10805/1653. Couper A. (1993). In Investigations of Surfaces and Interfaces, vol. IXA (eds. B. Rossiter and R. Baetzold). John Wiley, New York. Enders S., Kahl H. and Winkelmann J. (2005) Interfacial properties of polystyrene in contact with carbon dioxide. Fluid Phase Equil. 228–229, 511–522. Epelbaum M. B., Babashov I. V. and Salova T. P. (1973) Temperatures and pressures. Geochim. Int. J. 10(12), 343–345. Gardner J. E. (2012) Surface tension and bubble nucleation in phonolite magmas. Geochim. Cosmochim. Acta 76, 93–102.

S. Colucci et al. / Geochimica et Cosmochimica Acta 175 (2016) 113–127 Gardner J. E. and Ketcham R. A. (2011) Bubble nucleation in rhyolite and dacite melts: temperature dependence of surface tension. Contr. Mineral. Petrol. 162, 929–943. Gardner J. E., Ketcham R. A. and Moore G. (2013) Surface tension of hydrous silicate melts: constraints on the impact of melt composition. J. Volcanol. Geotherm. Res. 267, 68–74. Gardner J. E., Hilton M. and Carroll M. R. (2000) Bubble growth in highly viscous silicate melts during continuous decompression from high pressure. Geochim. Cosmochim. Acta. 64(8), 1473–1483. Gonnermann H. M. and Gardner J. E. (2013) Homogeneous bubble nucleation in rhyolitic melt: experimental and nonclassical theory. Geochem. Geophys. Geosyst. 14, 4758–4772. Harrison K. L., Johnston K. P. and Sanchez I. C. (1996) Effect of surfactants on the interfacial tension between supercritical carbon dioxide and polyethylene glycol. Languimir 12, 2637– 2644. Hilbert R., Tdheide K. and Frank E. U. (1981) PVT Data for water in the ranges 20 to 600 °C and 100 to 4000 bar. Ber. Bunsenges. Phys. Chem., 85–636. Iacono Marziano G., Schmidt B. C. and Dolfi D. (2007) Equilibrium and disequilibrium degassing of a phonolitic melt (Vesuvius AD 79 white pumice) simulated by decompression experiments. J. Volcanol. Geotherm. Res. 161, 151–164. IAPWS (International Association for the Properties of Water and Steam) (2007) Revised Release on the IAPWS Industrial Formulation 1997 for the Thermodynamic Properties of Water and Steam, Lucerne, Switzerland, August 2007. Jaeger Ph. T., Eggers R. and Baumgartl H. (2002) Interfacial properties of high viscous liquids in a supercritical carbon dioxide atmosphere. J. Supercrit. Fluids 24, 203–217. Kahl H. and Enders S. (2000) Calculation of surface properties of pure fluids using density gradient theory and SAFT-EOS. Fluid Phase Equil. 172, 27–42. Khitarov N. I., Lebedev Y. B., Dorfman A. M. and Bagdasarov N. S. (1979) Effect of temperature, pressure and volatiles on the surface tension of molten basalt. Geochem. Int. J. 16, 78–86. Lautze N. C., Sisson T. W., Mangan M. T. and Grove T. L. (2011) Segregating gas from melt: an experimental study of the Ostwald ripening of vapor bubbles in magmas. Contrib. Mineral. Petrol. 161, 331–347. Lensky N., Navon O. and Lyakhovsky V. (2002) Bubble growth during decompression of magma: experimental and theoretical investigation. J. Volcanol. Geotherm. Res. 129(1–3), 7–22. Machida H., Sato Y. and Smith, Jr., R. L. (2010) Simple modification of the temperature dependence of the SanchezLacombe equation of state. Fluid Phase Equil. 297, 205–209. Mader H. M., Llewellin E. W. and Mueller S. P. (2013) The rheology of two-phase magmas: a review and analysis. J. Volcanol. Geotherm. Res. 257, 135–158. Mangan M. T. and Sisson T. (2000) Delayed, disequilibrium degassing in rhyolite magma: Decompression experiments and implications for explosive volcanism. Earth Planet Sci. Lett. 183, 441–455. Mangan M. T. and Sisson T. (2005) Evolution of melt–vapor surface tension in silicic volcanic systems: Experiments with hydrous melts. J. Geophys. Res. 110, 1–9. Masotta M., Ni H. and Keppler H. (2014) In situ observations of bubble growth in basaltic, andesitic and rhyodacitic melts. Contrib. Mineral. Petrol, 167–976.

127

Mourtada-Bonnefoi C. C. and Laporte D. (1999) Experimental study of homogeneous bubble nucleation in rhyolitic magmas. Geophys. Res. Lett. 26, 3505–3508. Mourtada-Bonnefoi C. C. and Laporte D. (2002) Homogeneous bubble nucleation in rhyolitic magmas: an experimental study of the effect of H2O and CO2. J. Geophys. Res. 107(B4). Mourtada-Bonnefoi C. C. and Laporte D. (2004) Kinetics of bubble nucleation in a rhyolitic melt: an experimental study of the effect of ascent rate. Earth Planet Sci. Lett. 218, 521–537. Navon O. and Lyakhovsky V. (1998) Vesiculation processes in silicic magmas. Geol. Soc. Lond. Spec. Publ. 145, 27–50. NISTIR (National Insitute of Standards and Technology) 5078 (1998) Thermodynamic properties of water: tabulation from the IAPWS formulation 1995 for the thermodynamic properties of ordinary water substance for general and scientific use. Papale P., Moretti R. and Barbato D. (2006) The compositional dependence of the saturation surface of H2 O þ CO2 fluids in silicate melts. Chem. Geol. 229, 78–95. Park H., Thompson R. B., Lanson N., Tzoganakis C., Park C. B. and Chen P. (2007) Effect of Temperature and Pressure on Surface Tension of Polystyrene in Supercritical Carbon Dioxide. J. Phys. Chem. B 111, 3859–3868. Poser C. I. and Sanchez I. C. (1979) Surface tension theory of pure liquids and polymer melts. J. Coll. Interf. Sci. 69–3, 539–548. Poser C. I. and Sanchez I. C. (1981) Interfacial tension theory of low and high molecular weight liquid mixtures. Macromology 14, 361–370. Proussevitch A. A. and Sahagian D. L. (1998) Dynamics and energetics of bubble growth in magmas: analytical formulation and numerical modelling. J. Geophys. Res. 103, 18223–18251. Sahagian D. L., Anderson A. T. and Ward B. (1989) Bubble coalescence in basalt flows: comparison of a numerical model with natural examples. Bull. Volcanol. 52, 49–56. Sanchez I. C. and Lacombe R. H. (1976) An elementary molecular theory of classical fluids: pure fluids. J. Chem. Phys. 80, 2352. Sato Y., Yurugi M., Fujiwara K., Takishima S. and Masuoka H. (1996) Solubilities of carbon dioxide and nitrogen in polystyrene under high temperature and pressure. Fluid Phase Equil. 125, 129–138. Toramaru A. (1989) Vesciculation process and bubble size distributions in ascending magmas with constant velocities. J. Geophys. Res. 94, 17523–17542. Vargaftik N. B., Volkov B. N. and Voljak L. D. (1983) International tables of the surface tension of water. J. Phys. Chem. Ref. Data 12(3), 817–820. Vukalovich M. P., Zubarev W. N. and Alexandrov A. A. (1961) Experimental determination of specific volume of water to pressures at 1200 kg/cm2 and temperature 400–650 °C; 700– 900 °C. Teploenergetika 8(10), 79. Vukalovich M. P., Zubarev B. N. and Alexandrov A. A. (1962) Experimentelle Bestimmung des spezifischen Volumens von Wasserdampf bei 700 bis 900 °C und bis zu Drcken von 1200 kp/cm2 (russisch). Teploenergetika 9, 49. Walker D. and Mullins, Jr., O. (1981) Surface tension of natural silicate melts from 1,200–1,500 °C and implications for melt structure. Contrib. Mineral. Petrol 76, 455–462. Associate editor: Lawrence M. Anovitz