TWOLIQ.FOR: a FORTRAN77 program for simulating immiscibility in silicate liquids

PERGAMON

Computers & Geosciences 25 (1999) 151±159

TWOLIQ.FOR: a FORTRAN77 program for simulating immiscibility in silicate liquids H. Ma *, Y. Hu, T. Fang School of Materials Science and Technology, China University of Geosciences, Beijing 100083, People's Republic of China Received 9 September 1997; revised 7 July 1998

Abstract The program TWOLIQ.FOR is designed for predicting immiscibility in silicate liquids, by the thermodynamic Ho Ho criterion: a(Ai/T + CiP/T)XHo and for calculating compositions and amounts of the i rÿ a(aDiXi ÿBi)Xi conjugate liquids from oxide partition coecients between the coexisting Si- and Fe-rich melts, expressed as: ln(XSi i / Ho XFe i ) = ai/T + bi+ciP/T + adiXi . Where T and P denote temperature (in Kelvin) and pressure (in GPa), respectively, Xi mole fraction of oxide i, Ho, Si and Fe refer to homogeneous, Si- and Fe-rich melt phases, respectively and Ai to Di, ai to di are constants. Uncertainties of calculated oxide compositions in the liquids are 3.0±4.0 mol% for SiO2, Al2O3 and FeO, less than 1.0 mol% for the other oxides and predicted amounts around 1.0 mol% for the coexisting two liquids. Oreforming processes of magnetite±apatite deposits, therefore, can be numerically simulated by the program. # 1999 Elsevier Science Ltd. All rights reserved. Code available at http://www.iamg.org/CGEditor/index.htm Keywords: Silicate liquid; Immiscibility; Thermodynamics; Magnetite±apatite deposit; Numerical simulation

1. Introduction Liquid immiscibility has been demonstrated by both ®eld and experimental studies. Greig (1927) investigated the system FeO±SiO2 and found a second liquid in equilibrium with an almost pure SiO2 melt. Roedder (1951) discovered low-temperature liquid immiscibility in the system K2O±FeO±Al2O3±SiO2. Since then, experimental studies of immiscibility in a number of chemically distinct systems have been carried out to explore the e€ects of temperature, pressure, oxygen fugacity and chemical compositions on immiscibility (Visser et al., 1979a,b,c; Philpotts, 1982; Philpotts and Doyle, 1983; Hess and Wood, 1982; Naslund, 1983; Freestone and Powell, 1983).

* Corresponding author. E-mail: [email protected].

A major goal of theoretical petrology is the development of quantitative models for the magmatic di€erentiation processes that produce the observed chemical diversity of igneous rock suites. Currie (1972) presented a criterion for predicting immiscibility in silicate melts, but the lack of laboratory evidence and of a thermodynamic basis led to the abandonment by most later petrologists. Ghiorso et al. (1983) developed an expanded regular solution model for calculating the Gibbs free energy of silicate liquids, and used that to predict liquid immiscibility in a tholeiite, but this approach predicted coexisting immiscible liquid compositions, that deviated signi®cantly from the experimental results of Philpotts (1979). Upon examining the experimentally determined coef®cients of REE partitioning between immiscible silicate liquids of Ellison and Hess (1989) Nielsen and Gallahan (1990) argued that application of the Nielsen (1985) two-lattice melt model removes 90±96% of the

0098-3004/99/$ - see front matter # 1999 Elsevier Science Ltd. All rights reserved. PII: S 0 0 9 8 - 3 0 0 4 ( 9 8 ) 0 0 1 1 5 - 0

152

H. Ma et al. / Computers & Geosciences 25 (1999) 151±159

dependence of partitioning on melt compositions. Ma (1993) further found that the partition coecients of REE can be predicted within uncertainties around 1.5% of the experimentally determined values by the two-lattice melt model and an added term of temperature dependence. But the two-lattice melt model can not be applied to correctly predict liquid immiscibility in silicate liquids. TWOLIQ.FOR is a FORTRAN77 program incorporating the thermodynamic model and algorithms of Ma et al. (1998) for simulating liquid immiscibility in silicate liquids. The program is divided conceptually into two main functions: predicting the total apparent Gibbs free energy of silicate liquids and calculating chemical compositions and amounts of immiscible conjugate liquids. The functions for deriving oxygen fugacity from bu€er reactions (Ballhaus et al., 1991) and calibrating Fe2O3 and FeO contents in silicate melts (Kress and Carmichael, 1991) are also included in the program.

The equilibrium constant for Eq. (4) is de®ned as a Fe b Ho Ki ˆ …aSi i † …ai † =ai :

…5†

The free energy change of ith component for Eq. (4) must be zero, while the immiscible reaction is at equilibrium Dgi ˆ 0 ˆ DH 0i ÿ TDS 0i ‡ PDV 0i ‡ RT ln Ki :

…6†

Homogeneous, Si- and Fe-rich phases are all silicate liquids with various amounts of oxides. Having postulated a regular solution model, the ith component activity is de®ned as ln ai ˆ ln Xi ‡ fi =T ‡ PDVi =RT

…7†

ln gi ˆ ‰fi ‡ PDVi =RŠ=T,

…8†

The change of Gibbs free energy of an immiscible equilibrium can be written

where fi is the compositional dependence of ith component activity coecient in the melt, DVi denotes the di€erence between the partial molar volume of the ith oxide in silicate liquids and that in the pure ith component liquid. It is determined by ®tting calculations that the extent of fi Si, fi Fe and fi Ho variations with di€ering melt compositions have little di€erence, being incorporated into fi and similarly DVi Si, DVi Fe, DVi Ho and DVi 0 in Eq. (6) into DVi, therefore

DG ˆ DH 0 ÿ TDS 0 ‡ PDV 0 ˆ 0

ln Kdi ‡ DH 0i =RT ÿ DS 0i =R ‡ PDVi =RT

2. Thermodynamic model

…1†

0

assuming that DCp=0, DV is constant and the entropy of mixing is ideal. DH 0, DS 0 and DV 0 are the changes of the enthalpy, entropy and volume for the immiscible equilibrium reaction in the silicate melts. The total Gibbs free energy of a mixture of homogeneous and Si- and Fe-rich melts can be de®ned as X X Gˆ ni gi ˆ ni mi …i ˆ 1, n† …2† mi ˆ m0i ‡ RT ln ai ˆ m0i ‡ RT ln Xi ‡ RT ln gi

…3†

where mi 0 refers to the standard state chemical potential of the ith component, Xi, gi and ai are mole fraction, activity coecient and activity of the ith component, respectively and R the universal gas constant. The oxides are chosen as components to represent the compositions of melts. In a closed system, mass balance must be maintained for an immiscible reaction at equilibrium, i.e. X

Ho i

ˆ aX

Si i

‡ bX

Fe i

…a ‡ b ˆ 1:0†

…4†

where a and b are the mole fractions of Si- and Fe-rich liquid phases, respectively, Xi Ho, Xi Si and Xi Fe denote mole fractions of the ith component in homogeneous, Si- and Fe-rich liquids, respectively.

a Fe b Ho ÿ ln‰…gSi i † …gi † =…gi †Š ˆ 0:

…9†

Assuming fi is a linear function of the liquid compositions, we obtain Fe Ho a lnX Si i ‡ b ln X i ÿ ln X i

ˆ Ai =T ‡ Bi ‡ Ci P=T ‡

X Di X Ho i ,

…10†

where Ai, Bi and Ci are constants corresponding to DHi 0/R, DSi 0/R and DVi/R, respectively and aDiXi is the compositional dependence of the activity coecients. When an immiscible reaction is at equilibrium, the following conditions are satis®ed. 1. The chemical potentials of the ith component in the conjugate liquids are equal, i.e. Fe mSi i ˆ mi

…11†

and from Eq. (3), we have Si RT ln X Si i ‡ RT ln gi Fe ˆ RT ln X Fe i ‡ RT ln gi :

…12†

ÿ8.7815

11.0671 7.10842

ÿ1.5058 ÿ55.991

1.41108 3.18122 4.02085 11.2950

ÿ10.574 41.7399 4.95122

0.87420

ÿ0.2186 0.78158 ÿ1.4874

1.42321 ÿ4.4945 0.29619 5.12725

20.2554

0.97004 3.74236 3.49033 ÿ1.8264 3.21503

3.81819 9.99417 23.3753 4.05015

2.14780 15.4775 27.3440 3.06022 ÿ0.8892 3.60831 4.27522

2.42860

ÿ1.2553 ÿ9.8385 ÿ7.6604 ÿ2.2940

1.42321 3.45345

ÿ903.77 ÿ0.2369 0.02783 ÿ1257.9 1.34896

0.89122

0.40325 ÿ3.7149 ÿ1.9625 ÿ0.9129 ÿ1.6670

ÿ1.20131 2.39021 ÿ1.71627 ÿ0.87946 0 0.715328 ÿ3.66322

Al2O3 TiO2 SiO2

Table 1 Fitted parameters in Eq. (10)

The experimental data on immiscibility in silicate systems used in ®tting calculations were taken from Visser et al. (1979a,b,c), Dixon and Rutherford (1979), Ryerson and Hess (1980), Hess and Wood (1982), Philpotts (1982), Philpotts and Doyle (1983), Naslund (1983) and Freestone and Powell (1983). The compositional ranges of the investigated systems are as follows (wt%): SiO2 44.76±73.25, TiO2 0±12.70, Al2O3 2.06±18.99, FeO
0.90±40.90, MgO 0±10.30, MnO 0± 3.00, CaO 0±12.29, Na2O 0±6.12, K2O 0±11.17, P2O5 0±13.70. The temperatures range from 960 to 15508C, pressures from 1 atm to 1.5 GPa and the oxygen fugacities from in air to the IW bu€er. These experimental studies of immiscible equilibria as a function of temperature, pressure and composition have gone a long way towards establishing the broad outlines of the controls on the partitioning of major oxides between coexisting immiscible liquids. The data also set limits on the intensive and extensive variables, within which the program can model liquid immiscibility in natural magmas.

Fe2O3

3. Model limits and ®tting parameters

0.89122 ÿ2.8114 0.35970 6.79089 2.06001

FeO

MnO

where ai, bi and ci are constants corresponding to DHi 0/R, DSi 0/R and DVi/R, and adiXi, the compositional dependence of the ith oxide activity coecient. Constants Ai to Di and ai to di can be resolved by ®tting experimentally determined partition coecients of oxides in the immiscible conjugate melts to Eqs. (10) and (16).

ÿ0.3584 ÿ6.7457 ÿ1.2803 ÿ2.134 ÿ4.3889

…16†

ÿ3.06266

X di X Ho i

4.09419

Fe ln…X Si i =X i † ˆ ai =T ‡ bi ‡ ci P=T ‡

MgO

By analogy to Eq. (10), the coecient of the ith oxide partitioning between the two conjugate liquids can be ®tted to the following expression:

ÿ5419.6 ÿ2.9583

CaO

where DHi 0, DSi 0 and DVi 0 are the di€erences of molar enthalpy, entropy and volume of the ith oxide between the two conjugate liquids.

ÿ2480.2 2.24954 0.01263 ÿ0.3584 ÿ2.9667

…15†

ÿ3470.5 1.83403 0.00610 ÿ2.3333 ÿ7.8693 ÿ1.1697 ÿ9.4396 ÿ7.9183 2.46494 ÿ12.616 ÿ1.7105 ÿ0.6906 ÿ4.9179 4.41969

Fe ‡ PDV 0i =RT ÿ ln…gSi i =gi †,

Na2O

Fe 0 0 ln…X Si i =X i † ˆ DH i =RT ÿ DS i =R

ÿ514.50 0.31697

…14†

K2O

0 0 0 ÿRT ln gFe i ˆ DH i ÿ TDS i ‡ PDV i

ÿ3068.8 4.98002

Si Fe RT ln X Si i ‡ RT ln gi ÿ RT ln X i

ÿ255.04 0.13765

…13†

Ai Bi Ci DSiO2 DTiO2 DAl2 O3 DFe2 O3 DFeO DMnO DMgO DCaO DNa2 O DK2 O DP2 O5

G Si ˆ G Fe

P2O5

2. The Gibbs free energies of the coexisting liquids are equal, i.e.

153

ÿ2223.4 3.86944 0.01526 ÿ1.9539

H. Ma et al. / Computers & Geosciences 25 (1999) 151±159

ÿ33.646 21.0437

11.9163 ÿ10.171 3.81792

21.3443 ÿ33.374 ÿ14.977 7.15110 ÿ22.479 6.77789 12.3630

1.33260 38.8446 43.3430 19.0308 1.94907

19.7490

52.0016

ÿ1.2684

ÿ4.6760 ÿ10.685 ÿ2.3960 ÿ5.8983

14.1927 ÿ5.3031 7.66464 5.30037 90.9966

ÿ7.9283

ÿ599.13 5.48521 0.05396 14.5310 ÿ12.074 7.45673 ÿ17.409 ÿ26.043 4936.04 ÿ3.6747 ÿ0.0294 ÿ1588.5 0.90265 ÿ6048.0 0.34747 ÿ9920.8 5.20305 ÿ4928.2 34.3771

ÿ14.749 18.6494

P2O5 K2O Na2O CaO MgO

H. Ma et al. / Computers & Geosciences 25 (1999) 151±159

4. Prediction of immiscibility When immiscible silicate liquids are at equilibrium, we can de®ne X DG ˆ Dgi Xi ˆ 0 X Di X Ho ÿ a ln X Si i i

…17†

23.8501 ÿ326.94

17.1656 ÿ14.402 ÿ5.7269 ÿ58.978 6.90140

ÿ5203.7 5.24249 0.05393

ÿ3940.0 3.70936 0.02514 11.6677 ÿ11.771 7.87989 ÿ15.934 ÿ27.186 136.977 ÿ30.590 3.21513 7.39564 37.1314 3985.38 ÿ2.3565 ÿ0.0438 ÿ2.4042 34.0143 ÿ12.924 9.20668 14.3543 23.3747

13.7736

27.7561 3.59214

2.18274 ÿ5.4068

ÿ2.5735 ÿ4.1481

ai bi ci dSiO2 dTiO2 dAl2 O3 dFe2 O3 dFeO dMnO dMgO dCaO dNa2 O dK2 O dP2 O5

1333.88 ÿ0.9899 ÿ0.0175 0.73414 ÿ3.4815

At the moment that an immiscible reaction begins in a homogeneous silicate liquid, the amount of one of the two conjugate liquids, say the Fe-rich liquid, must be close to zero, i.e. b 1 0, a 1 1. The composition of the Si-rich liquid is then nearly the same to the homogeneous liquid, i.e. ln Xi Si1ln Xi Ho, and therefore we obtain X Dgi ˆ 0 ˆ Ai =T ‡ Bi ‡ Ci P=T ‡ Di X Ho …18† i and give the thermodynamic criterion of immiscibility in silicate liquids X DG ˆ Dgi Xi R0: …19†

4.62721 ÿ19.213 4.91313 ÿ7.2803 ÿ14.586 ÿ46.596 ÿ13.735 ÿ11.233 6.42954

ÿ0.4873 20.2199

MnO FeO

ÿ613.34

Ho ÿ b ln X Fe i ‡ ln X i :

Fe2O3

ÿ9622.2 7.5262

On the basis of the data set, the e€ects of composition, temperature and pressure on the partitioning behavior of major oxides were evaluated, and the parameters in Eqs. (10) and (16) were calculated by the least squares analysis (Tables 1 and 2).

Dgi ˆ Ai =T ‡ Bi ‡ Ci P=T ‡

Al2O3 TiO2 SiO2

Table 2 Fitted parameters in Eq. (16)

154

For all components in the system, the criterion becomes  X XX Ho …Ai =T ‡ Ci P=T †X Ho r ÿ D X ÿ B X Ho i i i i i : …20† Immiscibility in natural silicate magmas can therefore be predicted. For a data set of 239 immiscible experiments used in ®tting the calculation, 232 equilibrated homogeneous silicate melts are correctly predicted to be separating into Si- and Fe-rich liquids. The other seven are all samples deleted in the ®tting calculation owing to abnormally large deviations. Using Eq. (16) and constants in Table 2, the compositions of the Si- and Ferich liquids equilibrated with the homogeneous silicate melts are calculated. The results show standard deviations 3.0±4.0 mol% for SiO2, Al2O3 and FeO, less than 1.0 mol% for TiO2, MgO, CaO, Na2O, K2O and P2O5 and around 1.0 mol% for amounts of both Si-

H. Ma et al. / Computers & Geosciences 25 (1999) 151±159

155

Fig. 1. Correlation of experimentally determined oxide contents with those calculated from Eq. (16), showing program's ability for duplicating oxide contents in equilibrated immiscible liquids (after Ma et al., 1998). Diagonal line represents 1:1 correlation. Data source is given in text. Abscissa: experimental mol%; Ordinate: calculated mol%.

156

H. Ma et al. / Computers & Geosciences 25 (1999) 151±159

and Fe-rich liquids, compared with the experimentally equilibrated compositions. Fig. 1 shows ability of the program to predict compositions consistent with immiscible experiments. 5. Simulation of magnetite±apatite mineralization It was experimentally determined that liquid immiscibilities in lunar basalts, in the Skaergaard intrusion, in `Kiruna type' ore bodies, and in alkaline magmas are likely to be best developed at local liquidus±temperature minima associated with the iron±fayalite±silica, fayalite±magnetite±silica, magnetite±hematite±silica and fayalite±magnetite±feldspathoid or magnetite±hematite±feldspathoid equilibrium curves, respectively (Naslund, 1983). Immiscibilities between alkaline±gabbroic to ferro±pyroxenitic magmas and monzonitic to potassium-rich granitic magmas are ®eld analogues of the immiscibility in the system KAlSi3O8±NaAlSi3O8± FeO±Fe2O3±SiO2, and are responsible for formation of iron±titanium oxide and apatite rocks (Philpotts, 1967). The process for forming magnetite±apatite ore deposits can be numerically simulated by the program TWOLIQ.FOR. As an example, we took the magnetite±apatite ore deposit in Yangyuan pluton, Hebei province, north China. The ore bodies were inferred to be formed by liquid immiscibility (Hou, 1990). Occurrence of biotite±pyroxene±syenite with spherulitic structures was considered to be evidence of liquid immiscibility. Yangyuan pluton is a ring complex with an exposed area about 1.3 km2 and an isotopic age of 188.2 Ma by biotite K/Ar method. The pluton is composed of an alkaline pyroxenite suite and a syenite suite, rich in P2O5, TiO2 and K2O. The sequence of magmatic emplacement was pyroxenite±biotite pyroxenite, biotite±orthoclase±pyroxenite, biotite±pyroxene±syenite and ®nally syenite. The spherulitic biotite±pyroxene± syenite is ring shaped, intruding into biotite±pyroxenite and intruded by syenite. Both globule and mesostasis phases in the spherulitic biotite±pyroxene±syenites are SiO2 undersaturated. The globules are mainly composed of orthoclase megacrystals including pyroxene, biotite and magnetite, the mesostasis consists of pyroxene, biotite and apatite, with interstitial orthoclase. Chemical compositions and optical properties of pyroxene and biotite in both globule and mesostasis phases are similar. The content of apatite in the rocks ranges from 6 to 10%. It was assumed that the pressure at which immiscibility occurred was 0.5 GPa and the redox state of the magmas was bu€ered by FMQ bu€er. Calculations for three samples of the spherulitic biotite±pyroxene±syenites, resembling the hypothetic homogeneous melts in chemical composition, show that the total apparent

free energy changes DG associated with liquid immiscibility are between ÿ0.021 and ÿ0.103 at temperatures ranging from 950 to 12508C. Immiscibility was, therefore, likely to take place in the magmas. The simulated correlation of compositions and amounts of the conjugate melts with temperature by the program TWOLIQ.FOR are shown in Fig. 2. It can be seen from Fig. 2 that the concentrations of SiO2, Al2O3 and K2O are much higher in Si-rich melts, compared with Fe-rich melts and tend to increase, while FeO, MgO and CaO are lower and to decrease, with decreasing temperature. Therefore the simulation satisfactorily duplicates experimentally determined compositional trends of immiscible silicate melts (e.g. Philpotts, 1982; Naslund, 1983). At a temperature around 12508C, the amounts of the calculated Si- and Fe-rich melts are less than 6 mol% and more than 94 mol%, respectively (Fig. 2H) and concentrations of MgO, CaO and K2O in the two melts tend to converge (Fig. 2D±F), indicating that immiscibility is likely to take place at about 12508C. This is a slightly higher temperature than the 12008C inferred by Hou (1990) from homogeneous temperature measurements of glassy inclusions. When the temperature decreases to 11008C, the amounts of the Siand Fe-rich melts go up and down to about 25 and 75 mol%, respectively (Fig. 2H). The amounts of the two coexisting melts are consistent with the occurrence of the Si-rich phase as globules and the Fe-rich phase as mesostasis in the biotite±pyroxene±syenite. The simulated coexisting melts resemble the alkaline pyroxenite and syenite suites, respectively. It is easily seen from Fig. 2 that the magnetite±apatite deposits in the Yangyuan pluton may be formed mainly at temperature ranging from 1250 to 11508C (Fig. 2G). Below 11008C, solubility of P2O5 in the Ferich melts is dramatically decreased, certainly unfavorable to magnetite±apatite ore formation, while above 11508C, P2O5 is highly enriched in the Fe-rich liquids. This result is consistent with the occurrences of the Yangyuan magnetite±apatite bodies in the pyroxenite facies (Hou, 1990).

6. Program structure and operation The program TWOLIQ.FOR has been designed to be intuitive thereby making the program easy to use. The program structure and calculation algorithm are shown in Fig. 3. The program is divided into three sections. 1. Data processing: the program is designed to read oxide contents of a silicate melt by a subroutine READIN.FOR and intensive variables from keyboard. After the input pressure, iterative tempera-

H. Ma et al. / Computers & Geosciences 25 (1999) 151±159

157

Fig. 2. Correlation of predicted oxide compositions (A±G) and amounts (H) of immiscible conjugate liquids with temperature for biotite±pyroxene±syenite magmas from Yangyuan pluton, Hebei province (after Ma et al., 1998), showing program's ability for simulating ore-forming process of magmatic magnetite±apatite deposit by liquid immiscibility.

158

H. Ma et al. / Computers & Geosciences 25 (1999) 151±159

tine Gunmix.FOR and compositions and amounts of the conjugate liquids by subroutine LIQMIX.FOR. The results are output into a user assigned ®le. The data format can be rearranged by rewriting the write statements in the program. For TWOLIQ to run eciently, an IBM PC-compatible 486, DOS 6.0 or later is required. Those who request TWOLIQ should send an unformatted 3.5 0 diskette to the ®rst author, or an e-mail request, or download the code from FTP.IAMP.ORG by anonymous FTP.

Acknowledgements The generous support of the China National Education Committee grant (to MHW) is gratefully acknowledged. We would like to thank Professor Zhou Xunruo and Professor Du Yangsong for their critical readings and comments of an early Chinese version of this manuscript.

Appendix A

Fig. 3. Flowchart showing program structure and calculation algorithm.

ture and bu€er reaction are entered, the oxygen fugacity is derived from the expressions of Ballhaus et al. (1991) by a function Bu€o2.FOR. The contents of Fe2O3 and FeO in the melt are calibrated using the model of Kress and Carmichael (1991) by a subroutine Fe3clb.FOR. These subroutines are parts of a professional library PETROL.LIB, coded by the ®rst author (Ma and Tao, 1998). 2. Free energy prediction: using the subroutine Gunmix.FOR, the total apparent free energy of unmixing equilibrium in silicate liquids is derived from Eq. (10) and the constants in Table 1. Whether liquid immiscibility can be developed in the melt is then predicted using Eq. (20). 3. Calculation of compositions and amounts of the two conjugate liquids: using subroutine LIQMIX.FOR, compositions and amounts of the immiscible conjugate liquids are calculated according to Eq. (16) and constants in Table 2. The calculation is iteratively executed by selecting a new iterative temperature, bu€ered oxygen fugacity, calibrating Fe2O3 and FeO contents of the silicate melt, calculating total apparent free energy by subrou-

The TWOLIQ.FOR program source code is written in FORTRAN77 language and compiled by Microsoft FORTRAN Optimizing compiler Version 5.00 and linked by Microsoft Segmented-Excutable Linker Version 5.03. This highly structured TWOLIQ.FOR source code has 343 lines. Each subroutine or function has been tested independently. The program listing on the IAM G server is followed by an example of calculated result for the mesostasis melt of tholeiite sample (67670) from Kohala volcano, Hawaii (Philpotts, 1982).

References Ballhaus, C., Berry, R.F., Green, D.H., 1991. High pressure experimental calibration of the olivine±orthopyroxene±spinel geobarometer: implications for the oxidation state of the upper mantle. Contributions to Mineralogy and Petrology 107 (1), 27±40. Currie, K.L., 1972. A criterion for predicting liquid immiscibility in silicate melts. Nature Physical Science 240 (1), 66± 68. Dixon, S., Rutherford, M.J., 1979. Plagiogranites as late-stage immiscible liquids in ophiolite and mid-ocean ridge suites: an experimental study. Earth and Planetary Science Letters 45 (1), 45±60. Ellison, A.J.G., Hess, P.C., 1989. Solution properties of rare earth elements in silicate melts: inferences from immiscible liquids. Geochimica et Cosmochimica Acta 53 (8), 1965± 1974.

H. Ma et al. / Computers & Geosciences 25 (1999) 151±159 Freestone, I.C., Powell, R., 1983. The low temperature ®eld of liquid immiscibility in the system K2O±Al2O3±FeO±SiO2 with special reference to the join fayalite±leucite±silica. Contributions to Mineralogy and Petrology 82 (3), 291± 299. Greig, J.W., 1927. Immiscibility in silicate melts. American Journal of Science 5th Series 13 (1), 1±44. Ghiorso, M.S., Carmichael, I.S.E., Rivers, M.L., Sack, R.O., 1983. The Gibbs free energy of mixing of natural silicate liquids: an expanded regular solution approximation for the calculation of magmatic intensive variables. Contributions to Mineralogy and Petrology 84 (1), 107±145. Hess, P.C., Wood, M.I., 1982. Aluminum coordination in metaaluminous and peralkaline silicate melts. Contributions to Mineralogy and Petrology 81 (1), 103±112. Hou, Z., 1990. Origin of immiscibility in the magmas of Yangyuan and Fanshan ring complexes, Hebei province, with special reference to formation of the magnetite±apatite deposits. Deposit Geology 9 (2), 119±128 (in Chinese). Kress, V.C., Carmichael, I.S.E., 1991. The compressibility of silicate liquids containing Fe2O3 and the e€ect of composition, temperature, oxygen fugacity and pressure on their states. Contributions to Mineralogy and Petrology 108 (1/ 2), 82±92. Ma, H., 1993. Introduction to Thermodynamics in Crystalline Petrology, (in Chinese). Geological Publishing House, Beijing, 188 pp. Ma, H., Hu, Y., Yuan, J., Fang, T., 1998. Simulation of immiscibility in silicate liquids: thermodynamic model and numerical method. Earth Sciences 23 (2), 41±48 (in Chinese). Ma, H., Tao, C., 1998. Software of Thermodynamics in Crystalline Petrology, (in Chinese). Geological Publishing House, Beijing, ca. 350 pp, in press. Naslund, H.R., 1983. The e€ect of oxygen fugacity on liquid immiscibility in iron-bearing silicate melts. American Journal of Science 283 (10), 1034±1059.

159

Nielsen, R.L., 1985. A method for elimination of compositional dependence of trace element distribution coecients. Geochimica et Cosmochimica Acta 49 (8), 1775± 1779. Nielsen, R.L., Gallahan, W.E., 1990. In defense of the twolattice melt model: a comment on Ellison and Hess. Geochimica et Cosmochimica Acta 54 (8), 2293±2295. Philpotts, A.R., 1967. Origin of certain iron±titanium oxide and apatite rocks. Economic Geology 62 (3), 303±315. Philpotts, A.R., 1979. Silicate liquid immiscibility in tholeiitic basalts. Journal of Petrology 20 (1), 9±118. Philpotts, A.R., 1982. Compositions of immiscible liquids in volcanic rocks. Contributions to Mineralogy and Petrology 80 (2), 201±218. Philpotts, A.R., Doyle, C.D., 1983. E€ect of magma oxidation state on the extent of silicate liquid immiscibility in a tholeiitic basalt. American Journal of Science 283 (9), 967±986. Roedder, E., 1951. Low temperature liquid immiscibility in the system K2O±FeO±Al2O3±SiO2. American Mineralogist 36 (3/4), 282±286. Ryerson, F.J., Hess, P.C., 1980. The role of P2O5 in silicate melts. Geochimica et Cosmochimica Acta 44 (4), 611±624. Visser, W., Koster, V., Groos, A.F., 1979a. Phase relations in the system K2O±FeO±Al2O3±SiO2 at 1 atmosphere with special emphasis on low temperature liquid immiscibility. American Journal of Science 279 (1), 70±91. Visser, W., Koster, V., Groos, A.F., 1979b. E€ects of P2O5 and TiO2 on liquid±liquid equilibria in the system K2O± FeO±Al2O3±SiO2. American Journal of Science 279 (8), 970±988. Visser, W., Koster, V., Groos, A.F., 1979c. E€ect of pressure on liquid immiscibility in the system K2O±FeO±Al2O3± SiO2±P2O5. American Journal of Science 279 (10), 1160± 1175.