The COMAGMAT-5: Modeling the Effect of Fe–Ni Sulfide Immiscibility in Crystallizing Magmas and Cumulates

The COMAGMAT-5: Modeling the Effect of Fe–Ni Sulfide Immiscibility in Crystallizing Magmas and... ABSTRACT Details of the calibration and testing of the first ‘sulfide version’ of the COMAGMAT magma crystallization model (version 5·2, 2012–2014) are presented. The model’s updated empirical basis includes new mineral–melt geothermometers for olivine, plagioclase, high-Ca pyroxene, pigeonite, and orthopyroxene (calibrated at 1 atm. pressure), which are combined with a recently proposed Fe–Ni sulfide solubility model. This allows COMAGMAT-5 to be used for calculations of equilibrium and fractional crystallization of S-saturated and S-undersaturated magmas, including changes in the Fe/Ni ratio in silicate melts, femic minerals, and coexisting sulfides, as well as sulfide-silicate (±Fe–Ti oxides) proportions for multiply-saturated mineral assemblages. Based on our experience in tests of experimental data and modeling crystallization of mafic magmas, the possible range of application of COMAGMAT-5 may be extended up to 1–2 kbar pressure. The new model suggests a strong dependence of sulfide liquid immiscibility on the Ni content of the melt, as the increase of Ni is shown to decrease sulfide solubility, thus stabilizing the sulfide liquid in crystallizing mineral assemblages. This effect was tested for a Ni-rich sulfide-saturated glass dredged from the the southern Mid-Atlantic Ridge. The sulfide COMAGMAT was found to accurately predict the low S content observed in the glass (∼600 ppm); other existing ‘FeS’ solubility models yield higher values, mostly > 1000 ppm. In addition, using an experimentally studied high magnesia andesite composition, the results of calculations with two versions of COMAGMAT (3·72 & 5·2) were compared with those produced by MELTS family models. Application examples of COMAGMAT-5 include the modeling Fe–Ni sulfide saturation during equilibrium and fractional crystallization of ultramafic systems, approximating the most primitive magmas and cumulates from the Bushveld Complex in South Africa. INTRODUCTION The thermodynamic aspect of the genetic reconstructions of sulfide saturation history in mafic to ultramafic systems requires evaluating sulfide immiscibility timing and sulfide liquid composition at different stages of magmatic evolution. This involves constraints on the phase and chemical compositions of the crystallizing magma and, or, cumulates, particularly focusing on the equilibrium between the residual silicate melt and immiscible sulfide liquid. Some of these issues have been addressed in the ‘sulfide version’ of the COMAGMAT-5 magma crystallization model (Ariskin et al., 2009, 2013). This version combines an updated sulfide solubility model with newly-calibrated silicate mineral-melt geothermometers in a single algorithm designed to assess the effect of the Fe/Ni sulfide immiscibility on the compositional evolution of silicate melts, Fe–Mg minerals, and coexisting sulfides. Another feature of the new COMAGMAT is its ability to calculate sulfide-silicate (±Fe–Ti oxides) proportions for multiple-saturated mineral assemblages, and their effect on rock/mineral geochemistry during post-cumulus solidification of cumulate piles. Using these unique capabilities, the sulfide model has been used recently to quantify the geochemical evolution of an immiscible sulfide liquid which originated in olivine cumulates from the Yoko-Dovyren layered intrusion (Ariskin et al., 2016), as well as for quantifying temperature-compositional parameters of a parental magma calculated for the Chiney Complex in Southern Siberia, Russia (Gongalsky et al., 2016). In this study, we provide details on the calibration and other capabilities of the most recent COMAGMAT-5·2·2 model (http://geo.web.ru/∼ariskin/soft.html-id=comagmat.htm) and propose a methodology for its application to sulfide-saturated and sulfide-undersaturated systems. MAIN FEATURES OF THE EARLIER COMAGMAT MODELS COMAGMAT (Ariskin et al., 1993; Ariskin, 1999; Ariskin & Barmina, 2004) and MELTS (Ghiorso & Sack, 1995; Ghiorso et al., 2002; Gualda et al., 2012) are two families of popular thermodynamic programs, developed to model crystallization of mafic to ultramafic and silicic magmas. Both series of models are internally consistent and have some limitations when applied to volcanic systems and layered intrusions (Ghiorso, 1997; Ariskin & Barmina, 2004; Gualda & Ghiorso, 2015). Initially, COMAGMAT-3 was restricted to modeling low-pressure tholeiitic systems, such as mid-ocean ridge basalt (MORB) and large igneous province (LIP) magmas that follow the following order of crystallization: Ol → Pl → Aug → Pig (or Opx) → Fe–Ti oxides (Ariskin et al., 1993). Later it was updated to model crystallization of hydrous calc-alkaline magmas, starting from Ol + Aug cotectics at elevated pressures; see applications of COMAGMAT-3·5 to high-Mg magmas of the Klyuchevskoi volcano (Ariskin, 1999; Ariskin & Barmina, 2004) and to recent calculations for the parental magmas of the Bezymianny volcano in Kamchatka (Almeev et al., 2013). An important feature of COMAGMAT is that it allows for modeling crystallization of pigeonite and orthopyroxene, which is important for modeling the late crystallization stages of silica-oversaturated magmas. All previous versions of COMAGMAT (3·0–3·72, http://geo.web.ru/∼ariskin/soft.html-id=comagmat.htm) allowed for an alternative choice of ‘Pig’ or ‘Opx’ crystallization to model tholeiitic and calc-alkaline magmas. Thermodynamic fundamentals of COMAGMAT Incremental simulation of silicate melt equilibrium crystallization involves calculating mineral-melt equilibria in an evolving heterogeneous system, which consists of a residual melt (l) and crystallizing minerals (1 ≤ j ≤ M): nil=nibulk−∑j=1M∑r=1R(j)νi(r)jnrj, (1) where i is an oxide component (1 ≤ i ≤ N), r refers to end-members of a given j-mineral solid solution (1 ≤ r ≤ R(j)), nibulk,nil,nrj are the amounts of i-component in the system, melt, and crystallizing end-members of j-mineral, respectively, and νi(r)j represents the number of moles of i-component in one formula unit of an end-member. At a constant pressure, the mass-balance constraint (1) underlies calculations of the Gibbs free energy of a partially crystallized system at a set of independent (mainly P–T) parameters: G=∑i=1Nnilμil+∑j=1M∑r=1R(j)nrjμrj, (2) where μil and μrj are the chemical potentials of i-component in the melt and r-end-member of j-mineral. There are two main approaches to the search for the minimum Gibbs free energy at each step of crystallization (Frenkel & Ariskin, 1984). The first approach allows the isobaric potential at a given temperature to be minimized. This can be achieved through a number of algorithms based on nonlinear mathematical programming (e.g. Shvarov, 1981; Ghiorso, 1994; Karpov et al., 2001) utilizing an internally-consistent set of thermodynamic properties for the chosen components (e.g. Berman, 1988), thus ensuring that the model is applicable to a range of compositions and P–T conditions. An alternative approach is based on minimizing G at a given extent of crystallization (F), i.e. at fixed proportions of melt and crystals, searching for the equilibrium temperature (Frenkel & Ariskin, 1984; Ariskin et al., 1993). An advantage of this method is that it can account for the strong nonlinear dependence of T vs. F for multiply-saturated systems. Construction of crystallization trajectories as a function of the crystallization degree Fk=Fk−1+ΔF (where ΔF is an increment at each k-step of modeling) via small increments (e.g. ΔF ≤ 0·01) results in fast and accurate identification of inflection points which correspond to changes in crystallizing mineral assemblages. This capability is important for modeling the evolution of multi-mineral cotectics (e.g. Ol + Pl + Cpx → Ol + Pl + Cpx + Opx  → Ol + Pl + Cpx + Opx + Mt) characterized by a small dependence of their crystallization temperature on F. Another advantage of this approach is that it allows a finite-difference computational algorithm to be implemented when constructing magma differentiation models, which combine thermal and dynamical constraints of magma chambers with calculated crystallization sequences (Frenkel et al., 1989; Ariskin et al., 1993; Ariskin & Yaroshevsky, 2006). Basic principles of the COMAGMAT algorithm In a general form, the magma crystallization process can be described as a sequence of interdependent crystallization reactions producing end-members of rock-forming minerals EMrj (e.g. Fo, An, En) from the melt components Li (e.g. SiO2, MgO, NaAlO2 or Na2SiO3): ∑i=1Nνi(r)jLi=EMrj, (3) where νi(r)j are the stoichiometric coefficients of reaction (3), which correspond to the end-member components assumed in Eq. (1). Values of the equilibrium constants for the crystallization reactions (3) are given by the law of mass action: Krj=arj/∏i=1Naiνi(r)j, (4) where arj are the activities of end-members in the crystallizing minerals and ai are the activities of the melt components. In fact, the choice of thermodynamic components in the melt is arbitrary, as it depends on the accepted melt activity model. The lack of a comprehensive thermodynamic theory for multicomponent silicate melts leads to utilization of simplified descriptions of silicate melts, which involve a range of components unlikely to be physically present (e.g. NaAlO2; Nielsen & Dungan, 1983; or Fe2O3 and Na2SiO3; Ghiorso & Sack, 1995). The semi-empirical character of thermodynamic models for silicate melts is a feature of all current models of mineral–melt equilibria. The initial COMAGMAT model employed an ideal melt mixing model modified after Nielsen’s ‘two-lattice model’ (Ariskin et al., 1993). Assuming also ideal mixing for end-member components in mineral solutions, the link between the equilibrium constant (4) and the equilibrium temperature may be simplified to: ln Krj=ln xrj/∏i=1Nxiνi(r)j=−ΔG(3)/RgT=Arj/Trj+Brj, (5) where xrj are the concentrations of end-members in the crystallizing minerals, xi are the concentrations of the postulated melt components, and Arj>0 is the Arrhenius term representing the normalized enthalpy of r-phase melting, which is equivalent to the negative change in enthalpy of reaction (3), Arj=−ΔH(3)rj/Rg ( Rgis the gas constant). Modeling a k-step of equilibrium crystallization, Fk=Fk−1+ΔF initially involves the use of equation (1) to calculate a new melt composition ni(k)l by forming new amounts of the crystallizing mineral end-members nr(k)j (the total amount of new minerals is ΔF), using the same proportions as existed at the previous step Fk−1. This procedure can be described as a transition from the initial equilibrium state Fk−1 to a new non-equilibrated state due to a change in the melt composition following equation (1). Based on equations (1–5), a search for the minimum of the free energy function (2) for the new metastable system can be performed at the given crystallization degree Fk using a relationship described by Frenkel & Ariskin (1984): ∂G/∂nrj=−RgT ln Krj(o)+RgT ln Krj, (6) where Krj(o) is the reference state for pure endmembers r (i.e. forsterite or anorthite at their melting points) and Krj=Kr(k)j is calculated from (5) at each k-step of crystallization. Equation (6) is given to demonstrate that the decrease in the free energy of a metastable melt–mineral system at given conditions (P and Fk) can occur via iterative cycles of crystallization ( ∂nrj > 0) and/or dissolution ( ∂nrj < 0) of a number of end-members, leading to the equilibrium values of Kr(k)j following equation (4). Equation (6) forms the basis of the original COMAGMAT algorithm for free-energy minimization (Ariskin et al., 1993). The role of pseudo-liquidus mineral–melt temperatures Substituting K from Eq. (5), equation (6) can be re-written as ∂G/∂nrj=RgArj(T/Tj−1), (7) where Arj>0 represent parameters of corresponding mineral–melt geothermometers, T is the target (initially unknown) equilibrium temperature at step k, and Tj is the temperature of the current, notionally metastable state of each j-mineral ( Trj=Tj) with the current melt composition. Eq. (7) is used iteratively to determine the order and proportions of crystallization, including peritectic reactions at each step. For the completed derivation of expressions (6) and (7), see equations (25–45) in Ariskin & Barmina (2004). The COMAGMAT algorithm includes a special subroutine designed for calculating and comparing Tj (the so-called ‘pseudo-liquidus temperatures’, see Ariskin et al. (1993) and Danyushevsky & Plechov (2011)). Metastable mineral–melt temperatures Tj during iterations are calculated for each end-member component of every mineral present on the liquidus (i.e. nrj  > 0) following stoichiometry constraints: ∑r=1R(j)∏i=1Nxiνi(r)j exp ⁡(Arj/Tj+Brj)=1, (8) which require that the concentrations of end-members in each mineral sum to 1 ( ∑r=1R(j)xrj=1), where xrj=∏i=1Nxiνi(r)j exp ⁡(Arj/Tj+Brj) (9) The target temperature T must be within the range between the maximum and minimum of all pseudo-liquidus mineral temperatures Tj for the minerals currently present on the liquidus ( nj> 0 →nrj> 0). In this case, minimization in G is equivalent to minimization of the difference Tjmax⁡−Tj(nj>0)min⁡<εT, (10) where εT is the chosen computational accuracy during modelling. Detailed flowcharts of the main COMAGMAT algorithm describing its thermodynamic background are given by Ariskin & Barmina (2004). DEVELOPMENT OF COMAGMAT-5 COMAGMAT-5 (v. 5·2) is a new generation of the COMAGMAT programs designed to model crystallization of five silicate minerals (Ol–Pl–Aug–Pig–Opx), two Fe–Ti oxides (ilmenite and Fe–Ti spinel), and immiscible Fe–Ni sulfide liquids. This code was essentially rewritten, using more sophisticated compiling programs, as compared to previous versions of COMAGMAT (Ariskin et al., 1993; Ariskin, 1999); however, it preserved the main features of the basic algorithm described above. It includes a set of updated 1-atm mineral–melt geothermometers for olivine, plagioclase, augite, pigeonite, and orthopyroxene. Several new subroutines have been added, including SULSAT, which estimates the solubility of Fe–Ni sulfides based on the model of Ariskin et al. (2013). COMAGMAT-5 is also capable of modeling simultaneous crystallization of three pyroxenes including augite, pigeonite, and orthopyroxene. New silicate mineral–melt geothermometers The term ‘geothermometer’ is applied to a system of thermodynamic equations (5) which describe the partitioning of major and some minor elements between rock-forming minerals and silicate melts over a range of P–T–fO2–H2O conditions (Ariskin et al., 1993). The geothermometers are calibrated on a large dataset of melting and crystallization experiments in natural and synthetic systems ranging from ultramafic and high-Mg basalts to ferrobasalts and dacites (i.e. the INFOREX experimental database). Now the INFOREX melting-experiment database provides information on 15,418 sub-liquidus runs from 442 experimental studies published in the literature, see the INFOREX Reference List in Supplementary Data Electronic Appendix 1. The available information includes pressure, temperature, fO2, run duration, containers, 17415 compositions of 40 minerals and 8404 melt (glass) compositions, as well as signatures of nominally anhydrous and hydrous experiments in a wide range of natural and synthetic systems (Ariskin & Barmina, 2004; Nikolaev et al., 2016). Ariskin & Barmina (2004) have demonstrated that this approach results in geothermometers with a general precision of 10–15°C for each silicate mineral; these geothermometers can reproduce liquidus mineral compositions with a precision of 1–3 mol% for olivine and pyroxenes, and 2–4 mol% for plagioclase. Two major constraints underlie this approach. First, all mineral–melt geothermometers must be calibrated with respect to the equilibrium constants (5): ln Krj=a/T+b1 ln cl1+b2 ln cl2+⋯+d, (11) where Krj characterizes the reaction of r-end-member formation from melt components (3); cl1,cl2,… are empirically selected melt structural-chemical parameters; and a, b1,b2,… and d are the regression coefficients. Second, the Krj values need to be calculated using the same melt component activity model for each end-member and mineral. The above constraints are derived from the fundamental thermodynamic relationships (Ghiorso, 1987) to ensure that the crystallization model represents an internally-consistent formulation. Ignoring these constraints by selecting an arbitrary system of mineral–melt partitioning equations (Danyushevsky & Plechov, 2011) may result in incorrect proportions of crystallized minerals and large errors in the calculated liquid lines of descent (Ghiorso, 1987; Ariskin & Barmina, 2004). The Krj values in Eq. (11) are calculated using a two-lattice melt components model (Ariskin & Barmina, 1990) modified after Nielsen & Dungan (1983) and Nielsen (1990). Use of melt structural-chemical parameters cl (e.g. Si/O, Al/O, Al/Si) significantly improves the fit to the experimental data (Ariskin, 1999). The calibration data set To update geothermometers for olivine, plagioclase, augite, and pigeonite, we used an expanded set of coexisting mineral–melt compositions. The set includes ∼270 runs conducted at 1 atm (thus nominally anhydrous) for > 48 hours at temperatures from 1050–1250°C over a wide range of oxygen fugacity (IW - NNO + 1). These experimental glasses contain 7 ≤ FeO ≤ 18 wt %, 45 ≤ SiO2 ≤ 60 wt %, and 2 ≤ Na2O + K2O ≤ 5 wt %. This experimental array includes ∼150 mineral–melt pairs for olivine, 187 pairs for plagioclase, 125 pairs for high-Ca Cpx, and 43 pairs for pigeonite. To calibrate the new orthopyroxene–melt geothermometer an extra set of 96 experiments was added, which covers higher temperatures (to 1300°C) and more Si-rich compositions (48 ≤ SiO2 ≤ 68 wt %, see Supplementary Data Electronic Appendix 1; supplementary data are available for downloading at http://www.petrology.oxfordjournals.org). The calibration dataset was limited to 1065°C ≤T ≤1300°C as it is rare for volcanic suite and layered intrusion parental magmas to have higher temperatures (e.g. Ariskin & Barmina, 2004). The lower temperature range is included in the calibration dataset because protracted magma fractionation observed in many differentiated intrusions leads to highly evolved compositions (e.g. in the Skaergaard intrusion, the Layered Series rocks are composed of mineral assemblages corresponding to temperatures ∼1115–1080°C; McBirney, 1996; Ariskin, 2003). Since COMAGMAT-3 was calibrated at higher temperatures (∼1125–1350°C), calculations at and below 1100°C often resulted in over-expansion of the Ca-clinopyroxene field at the expense of plagioclase and pigeonite. Moreover, plagioclase compositions calculated for the most evolved melts were too sodic and potassic. Methods of calibration and results Calibrations of olivine and plagioclase geothermometers were carried out in the form of Eq. (11) by a least-squares method, using a two-stage procedure of multiple regression analysis (Ariskin, 1999; Ariskin & Barmina, 2004). For high-Ca clinopyroxene, pigeonite, and orthopyroxene a modified method was used, which is based on multiple regression analysis for normalized values: ln (Krj)′=ln Krj−aPxj/T=b1 ln cl1+b2 ln cl2+⋯+d, (12) where the Arrhenius terms aPxj for enstatite (En), ferrosilite (Fs), wollastonite (Wo), and the Al-component are postulated to be the same for all pyroxenes. Detailed comments about both approaches are given in Supplementary Data Electronic Appendix 2 (Part 1). Note that in addition to Mg, Fe, Ca, and Mn end-members, geothermometers for olivine and pyroxenes include the Ni end-members NiSi0·5O2 and NiSiO3, leading to a total system of 22 end-member geothermometers (Table 2·1 and Part 2 in Supplementary Data Electronic Appendix 2). The equations for each end-member were then combined into five mineral–melt equilibria models for olivine, plagioclase, augite, pigeonite, and orthopyroxene (Ariskin & Barmina, 2004). The proposed mineral–melt models have uncertainties for T of 10–15°C and for mineral compositions of ∼0·7–2 mol% of the end-member concentrations in olivine, orthopyroxene and augite, and of 3–4 mol% for pigeonite and plagioclase. Table 1: Results of sulfide solubility calculations for the BTJ glass at P = 1 atm SCSS model Melt parameters QFM QFM + 0·5 QFM + 1 Ni, ppm H2O, ppm T, °C SCSS, ppm T, °C SCSS, ppm T, °C SCSS COMAGMAT-5 (Ariskin et al., 2013) 310 0 1258·2 514* 1256·9 566* 1255·3 710* 0 0 1256·4 1241* 1255·0 1368* 1253·3 1710* Li & Ripley, 2009 0 515 1257 1054 1257 1042 1257 1028 Liu et al., 2007 0 515 1257 1027 1257 1002 1257 973 Fortin et al., 2015 (Model B) 0 515 1257 1121 – – – – SCSS model Melt parameters QFM QFM + 0·5 QFM + 1 Ni, ppm H2O, ppm T, °C SCSS, ppm T, °C SCSS, ppm T, °C SCSS COMAGMAT-5 (Ariskin et al., 2013) 310 0 1258·2 514* 1256·9 566* 1255·3 710* 0 0 1256·4 1241* 1255·0 1368* 1253·3 1710* Li & Ripley, 2009 0 515 1257 1054 1257 1042 1257 1028 Liu et al., 2007 0 515 1257 1027 1257 1002 1257 973 Fortin et al., 2015 (Model B) 0 515 1257 1121 – – – – * COMAGMAT-5 calculates the total sulfur content as SCSS (Sulfur Content at Sulfide Saturation), by adding S6+ using the model of Jugo (2009). Other models calculate only S2- as SCSS. Note that Model B by Fortin et al. (2015) does not take into account the effect of fO2. The major element composition of BTJ is given in Kamenetsky et al. (2013). Table 1: Results of sulfide solubility calculations for the BTJ glass at P = 1 atm SCSS model Melt parameters QFM QFM + 0·5 QFM + 1 Ni, ppm H2O, ppm T, °C SCSS, ppm T, °C SCSS, ppm T, °C SCSS COMAGMAT-5 (Ariskin et al., 2013) 310 0 1258·2 514* 1256·9 566* 1255·3 710* 0 0 1256·4 1241* 1255·0 1368* 1253·3 1710* Li & Ripley, 2009 0 515 1257 1054 1257 1042 1257 1028 Liu et al., 2007 0 515 1257 1027 1257 1002 1257 973 Fortin et al., 2015 (Model B) 0 515 1257 1121 – – – – SCSS model Melt parameters QFM QFM + 0·5 QFM + 1 Ni, ppm H2O, ppm T, °C SCSS, ppm T, °C SCSS, ppm T, °C SCSS COMAGMAT-5 (Ariskin et al., 2013) 310 0 1258·2 514* 1256·9 566* 1255·3 710* 0 0 1256·4 1241* 1255·0 1368* 1253·3 1710* Li & Ripley, 2009 0 515 1257 1054 1257 1042 1257 1028 Liu et al., 2007 0 515 1257 1027 1257 1002 1257 973 Fortin et al., 2015 (Model B) 0 515 1257 1121 – – – – * COMAGMAT-5 calculates the total sulfur content as SCSS (Sulfur Content at Sulfide Saturation), by adding S6+ using the model of Jugo (2009). Other models calculate only S2- as SCSS. Note that Model B by Fortin et al. (2015) does not take into account the effect of fO2. The major element composition of BTJ is given in Kamenetsky et al. (2013). Table 2: Compositions of magmatic melts and cumulate rocks used for testing COMAGMAT and MELTS on experimental data and example calculations with COMAGMAT-5 Melt components, wt % HMA-85-41c Proposed Bushveld magmas Cumulate rocks B1* UM** UM + 56 wt % Ol-91·6 Ol pyr*** SiO2 57·79 56·09 53·38 46·29 46·03 TiO2 0·60 0·28 0·39 0·17 0·10 Al2O3 14·46 11·31 9·73 4·28 4·86 FeO* 5·74 9·17 10·19 9·18 9·21 MnO 0·11 0·17 0·17 0·16 0·16 MgO 9·14 13·58 19·00 36·39 35·13 CaO 8·17 6·34 5·36 2·36 3·04 Na2O 3·11 1·43 1·34 0·59 0·35 K2O 0·71 1·05 0·54 0·24 0·14 P2O5 0·15 0·07 0·08 0·04 0·02 NiO – 0·04 0·09 0·32 0·29 Cr2O3 – 0·25 – – 0·90 Total 99·98 99·78 100·27 100 100·23 S, ppm – 438**** 368**** 162**** – Melt components, wt % HMA-85-41c Proposed Bushveld magmas Cumulate rocks B1* UM** UM + 56 wt % Ol-91·6 Ol pyr*** SiO2 57·79 56·09 53·38 46·29 46·03 TiO2 0·60 0·28 0·39 0·17 0·10 Al2O3 14·46 11·31 9·73 4·28 4·86 FeO* 5·74 9·17 10·19 9·18 9·21 MnO 0·11 0·17 0·17 0·16 0·16 MgO 9·14 13·58 19·00 36·39 35·13 CaO 8·17 6·34 5·36 2·36 3·04 Na2O 3·11 1·43 1·34 0·59 0·35 K2O 0·71 1·05 0·54 0·24 0·14 P2O5 0·15 0·07 0·08 0·04 0·02 NiO – 0·04 0·09 0·32 0·29 Cr2O3 – 0·25 – – 0·90 Total 99·98 99·78 100·27 100 100·23 S, ppm – 438**** 368**** 162**** – * HMA-85–41c is a high-Mg andesite used in melting experiments by Grove et al. (2003). Proposed magmas and cumulate rocks represent mafic to ultramafic (UM) compositions assumed to be parental for the Lower Zone and the Basal Ultramafic Sequence of the Bushveld Complex (Wilson, 2012, 2015). See column 4 in table 3 of Wilson (2012), **UM liquid composition from Schwerin fold chill listed as CH12/7&14UM in column 3 of table 4 in Wilson (2015), ***Olivine pyroxenite listed as CH7/91 in table 3 of Wilson (2015). ****S concentrations used in the COMAGMAT-5 modeling: 438 ppm S in B1 is taken from Barnes et al. (2010), 368 ppm initial S concentration estimated based on allowing UM to crystallize 16% olivine, which produces S concentration of 438 ppm, similar to that seen in B1; S in ‘UM + 56% Ol-91·6’ corresponds to a mixture of 44 wt.% of the UM composition with 56 wt.% of olivine mg#91.6. Table 2: Compositions of magmatic melts and cumulate rocks used for testing COMAGMAT and MELTS on experimental data and example calculations with COMAGMAT-5 Melt components, wt % HMA-85-41c Proposed Bushveld magmas Cumulate rocks B1* UM** UM + 56 wt % Ol-91·6 Ol pyr*** SiO2 57·79 56·09 53·38 46·29 46·03 TiO2 0·60 0·28 0·39 0·17 0·10 Al2O3 14·46 11·31 9·73 4·28 4·86 FeO* 5·74 9·17 10·19 9·18 9·21 MnO 0·11 0·17 0·17 0·16 0·16 MgO 9·14 13·58 19·00 36·39 35·13 CaO 8·17 6·34 5·36 2·36 3·04 Na2O 3·11 1·43 1·34 0·59 0·35 K2O 0·71 1·05 0·54 0·24 0·14 P2O5 0·15 0·07 0·08 0·04 0·02 NiO – 0·04 0·09 0·32 0·29 Cr2O3 – 0·25 – – 0·90 Total 99·98 99·78 100·27 100 100·23 S, ppm – 438**** 368**** 162**** – Melt components, wt % HMA-85-41c Proposed Bushveld magmas Cumulate rocks B1* UM** UM + 56 wt % Ol-91·6 Ol pyr*** SiO2 57·79 56·09 53·38 46·29 46·03 TiO2 0·60 0·28 0·39 0·17 0·10 Al2O3 14·46 11·31 9·73 4·28 4·86 FeO* 5·74 9·17 10·19 9·18 9·21 MnO 0·11 0·17 0·17 0·16 0·16 MgO 9·14 13·58 19·00 36·39 35·13 CaO 8·17 6·34 5·36 2·36 3·04 Na2O 3·11 1·43 1·34 0·59 0·35 K2O 0·71 1·05 0·54 0·24 0·14 P2O5 0·15 0·07 0·08 0·04 0·02 NiO – 0·04 0·09 0·32 0·29 Cr2O3 – 0·25 – – 0·90 Total 99·98 99·78 100·27 100 100·23 S, ppm – 438**** 368**** 162**** – * HMA-85–41c is a high-Mg andesite used in melting experiments by Grove et al. (2003). Proposed magmas and cumulate rocks represent mafic to ultramafic (UM) compositions assumed to be parental for the Lower Zone and the Basal Ultramafic Sequence of the Bushveld Complex (Wilson, 2012, 2015). See column 4 in table 3 of Wilson (2012), **UM liquid composition from Schwerin fold chill listed as CH12/7&14UM in column 3 of table 4 in Wilson (2015), ***Olivine pyroxenite listed as CH7/91 in table 3 of Wilson (2015). ****S concentrations used in the COMAGMAT-5 modeling: 438 ppm S in B1 is taken from Barnes et al. (2010), 368 ppm initial S concentration estimated based on allowing UM to crystallize 16% olivine, which produces S concentration of 438 ppm, similar to that seen in B1; S in ‘UM + 56% Ol-91·6’ corresponds to a mixture of 44 wt.% of the UM composition with 56 wt.% of olivine mg#91.6. Figure 1 shows the fit of the new models to the experimental data used for calibration. The new olivine model was also tested at higher temperatures on a set of 66 experiments at 1220°C  < T <1350°C (9–16 wt % MgO in the melt, 83–89% Fo), which were not included in the calibration dataset (Fig. 1a), demonstrating that the new model for olivine is applicable to at least 1350°C. The complete set of calculated temperatures and compositions and experimental data is given in Supplementary Data Electronic Appendix 1. The new calibrations for plagioclase and clinopyroxene resulted in a better fit of the calculated mineral compositions to the experimental data than achieved by the previous version of COMAGMAT-3.5 (Fig. 1c, d). Fig. 1. View largeDownload slide Comparison of the accuracy of the COMAGMAT-5 and COMAGMAT-3.5 calculations for equilibrium temperatures and mineral compositions. The dashed lines represent equal values. The temperature–compositional calculations by COMAGMAT were conducted on the same experimental arrays as those extracted from the INFOREX database (Ariskin & Barmina, 2004) and used in the calibration of mineral–melt geothermometers for the COMAGMAT-5 model, see details in the text (‘The calibration data set’ section). To test the validity of the model at higher temperatures, an additional set of 66 experimental Ol–melt pairs were used (1220°C < T <1350°C, 9–16% MgO in the melt, 83–89% Fo, see crosses in Fig. 1a). Note that the previous version of COMAGMAT-3.5 (Ariskin, 1999) worked well for low-alkaline tholeiitic systems; however, it essentially over-estimated the composition of plagioclase modeled in evolved crystallization products of mildly-alkaline magmas. Fig. 1. View largeDownload slide Comparison of the accuracy of the COMAGMAT-5 and COMAGMAT-3.5 calculations for equilibrium temperatures and mineral compositions. The dashed lines represent equal values. The temperature–compositional calculations by COMAGMAT were conducted on the same experimental arrays as those extracted from the INFOREX database (Ariskin & Barmina, 2004) and used in the calibration of mineral–melt geothermometers for the COMAGMAT-5 model, see details in the text (‘The calibration data set’ section). To test the validity of the model at higher temperatures, an additional set of 66 experimental Ol–melt pairs were used (1220°C < T <1350°C, 9–16% MgO in the melt, 83–89% Fo, see crosses in Fig. 1a). Note that the previous version of COMAGMAT-3.5 (Ariskin, 1999) worked well for low-alkaline tholeiitic systems; however, it essentially over-estimated the composition of plagioclase modeled in evolved crystallization products of mildly-alkaline magmas. Incorporation of Fe–Ni sulfide solubility model Recently, an additional subroutine (SULSAT) has been developed for calculating Sulfur Contents at Sulfide Saturation (SCSS, Campbell & Naldrett, 1979) in mafic magmas, as a function of pressure, temperature, fO2, major components, and Ni content in the melt (Ariskin et al., 2013). In contrast to other SCSS models, which approximate the immiscible sulfides by a pure FeS pyrrhotite liquid, see comments and references in Supplementary Data Electronic Appendix 2 (Part 3), we argued for a strong effect of minor components (such as Ni and Cu) which may be present in the natural sulfide liquids in amounts up to 10–20 wt %. Our SCSS model proposes the existence of positively charged Fe–Ni sulfide complexes in the silicate melt, which are formed as a result of complexation reactions between the sulfide-forming ions (Fe2+, Ni2+, S2-) and (Fe, Ni)S species. The proposed mechanism of sulfide solubility was calibrated and tested against 290 nominally anhydrous experimental and natural glass compositions (both Ni-free and Ni-bearing), producing a good fit between the observed and calculated variations of S content vs. temperature and Fe content in the melt (Ariskin et al., 2013). Integration of the SULSAT model into COMAGMAT-5 Incorporation of the SULSAT subroutine into the COMAGMAT-5 model required the addition of Ni and S as the principal thermodynamic components. This resulted in several changes to the algorithm, which allowed for combining the Fe–Ni sulfide model with equilibrium temperature and composition modeling for silicate minerals and Fe–Ti oxides. For initially S-undersaturated compositions, the algorithm uses SCSS calculations to determine the onset of sulfide immiscibility. This is achieved by comparing the modeled S content in the melt (which behaves as an incompatible component until sulfide saturation is reached) with the calculated SCSS at each stage of crystallization (e.g. Barnes & Lightfoot, 2005; Li & Ripley, 2005; Barnes, 2007; Barnes et al., 2010). When the calculated SCSS is the same as the S content in the melt, it is taken to indicate that the system has become sulfide saturated (see figs 9 and 10 in Ariskin et al. (2013)). Starting from this point, mass-balance constraints are used to calculate both the proportions and the compositions of all phases, including Fe–Ni sulfide. The algorithm also takes into account possible S over-saturation at the onset of crystallization, thus allowing for modeling systems containing an excess of sulfide phases. To account for the effect of fO2 on S speciation in the melt, we followed the approach of Jugo (2009), which allows for calculating the proportion of S6+ in a silicate melt. The value of bulk SCSS becomes a function of fO2 and increases with fO2 for a given major element composition. This is due to the addition of some sulfate component (as S6+) to the sulfide species ((Fe, Ni)S-complexes) at more oxidized conditions. The main options and capabilities of COMAGMAT-5 There are three main options in COMAGMAT-5: (i) calculating sulfur solubility for a given temperature and melt composition; (ii) calculating equilibrium temperatures for silicate and oxide minerals and their compositions for a given melt composition; and (iii) modeling crystallization processes for a given magma composition. A detailed overview of the main options and the updated structure of COMAGMAT-5 is given as Supplementary Data Electronic Appendix 2 (Part 4); examples of output files are available as spreadsheets in Supplementary Data Electronic Appendix 1. The first option allows SCSS to be calculated for a number of melt compositions, if their temperatures and fO2 are known. The SCSS calculations can be performed at a range of pressures up to 1 GPa. The second option allows both the temperature of crystallization and the composition of the mineral on the liquidus at 1 atm pressure to be calculated for a given melt composition. This routine is useful for testing the accuracy of the thermometric and compositional calculations in COMAGMAT-5 using experimental data. The third option allows crystallization to be modeled for multiply-saturated magmatic systems, including five silicate minerals (olivine, plagioclase, augite, pigeonite, orthopyroxene), two Fe–Ti oxides (ilmenite, magnetite), and a Fe–Ni sulfide liquid. Calculations can be performed for the cases of: (1) equilibrium crystallization in a closed system; (2) perfect fractional crystallization (assuming complete separation of crystallizing solids from the melt); and (3) an intermediate case where a certain portion of crystallizing minerals are postulated to re-equilibrate with the residual melt, whereas the remainder of the minerals are removed from the system (Nielsen, 1990; Ariskin & Barmina, 2004; Danyushevsky & Plechov, 2011). Changes in the user-interface and reorganization of the input/output files made COMAGMAT-5 more flexible and easy to use (see examples of its interface in Supplementary Data Electronic Appendix 2, Part 4). Currently, COMAGMAT-5 is the only algorithm capable of modeling co-crystallization of Aug±Pig±Opx, including calculations of peritectic relationships between orthopyroxene, pigeonite, and olivine, as well as between magnetite and ilmenite in plagioclase-pyroxene saturated systems. As compared to previous versions of COMAGMAT-3.n (Ariskin et al., 1993; Ariskin, 1999), the updated model for orthopyroxene provides more accurate calculations of its melting/crystallization proportions relative to olivine. Two examples of testing the COMAGMAT-5 on natural and experimental data are given in the section below. VERIFICATION OF THE SULFIDE COMAGMAT Test using a high-Ni natural glass composition The SCSS model used in COMAGMAT-5 (Ariskin et al., 2013) predicts that Ni has a pronounced effect on decreasing SCSS, despite its relatively low concentrations in silicate melts. Here we present a demonstration of the effect of Ni on SCSS using a sulfide-saturated glass composition dredged from the southern Mid-Atlantic Ridge near the Bouvet Triple Junction, sample S18–60/1 (the BTJ glass, Kamenetsky et al., 2001, 2013). This is a high-Mg (8·2 wt % MgO, mg# 67·1) and high-Ni (310 ppm) andesite (57·3 wt % SiO2) which is depleted in Cu (33 ppm) and S (on average 600 ppm), compared to typical MORB glasses (Mathez, 1976; Ariskin et al., 2013). A distinctive feature of this sample is the presence of numerous spherical Ni-rich sulfide globules, which have been interpreted to indicate sulfide-saturation of the melt (Kamenetsky et al., 2013). A detailed examination of the largest globule (∼180 μm in diameter), combining X-ray mapping and microprobe data, allowed the authors to estimate the average Ni and Fe contents in the sulfide melt; this melt is characterized by a high Ni/(Ni + Fe) molar ratio of ∼0·40. Rare micro-phenocrysts of olivine (Fo86·9) suggest the BTJ melt has olivine on its liquidus, thus allowing its liquidus temperature to be estimated by using olivine-melt geothermometers. Calculations of SCSS for the BTJ glass Falloon et al. (2007) have demonstrated that the olivine-melt geothermometer of Ford et al. (1983) accurately estimates crystallization temperatures of MORB melts at low pressure. Calculations for the BTJ glass, using the model of Ford et al. (1983) and COMAGMAT-5 at P = 1 atm and QFM buffer, resulted in 1255·2°C (Fo86·3) and 1258·2°C (Fo85·9), respectively. The ‘Sulfur solubility’ option of COMAGMAT-5 (see Supplementary Data Electronic Appendix 2) was used to calculate the SCSS for the BTJ glass composition to compare the calculated values with the measured S content in the glass (Table 1). The calculations were performed at the calculated liquidus temperature at 1 atm pressure and fO2 ranging from QFM to QFM + 1. Relatively oxidized conditions were selected to account for the oxygen contents observed in the sulfide globules from BTJ: 0·33–0·55 wt % O (Kamenetsky et al., 2013). Assuming 310 ppm Ni in the melt, COMAGMAT-5 calculated bulk S solubility ranging from 514 ppm S at QFM (4·5 % S presented as sulfate species S6+), to 566 ppm S at QFM + 0·5 (12·5% S6+), and 710 ppm S at QFM + 1 (29·3% S6+), see Table 1. The value at QFM + 0·5 is close to the 600 ppm S measured in the BTJ glass. Similar calculations based on the Li & Ripley (2009) equation resulted in SCSS values around 1040 ppm. We also tested the Liu et al. (2007) equation, which postulates a strong ‘dropping effect’ of H2O on SCSS at H2O < 0·1 wt%, see figure 3 in Ariskin et al. (2013). Using the H2O content determined in the BTJ glass (515 ppm H2O), these calculations resulted in the SCSS ranging from 1027 to 973 ppm (Table 1). The use of Model-B from Fortin et al. (2015) at the same parameters resulted in ∼1121 ppm S. As a result, none of the SCSS models, except COMAGMAT-5, was found to replicate the low sulfide solubility observed in the high-Ni BTJ glass. If we assume that the BTJ glass is Ni-free, COMAGMAT-5 calculates an SCSS of >1240 ppm (Table 1), which is similar to other models that do not take Ni content in the silicate melt into account. Thus, the 1·5–2 times overestimation of the SCSS based on pure FeS models is clearly due to the high Ni content observed in the glass. Results of COMAGMAT-5 modeling at QFM + 0·5 indicate that this melt should be in equilibrium with a high Ni sulfide (Ni/(Ni + Fe)=0·370), which is consistent with the Ni/(Ni + Fe)∼0·408 estimated for the BTJ sulfide globule (Kamenetsky et al., 2013). Test using experimental data The second test was carried out using the ‘Modeling crystallization’ option of COMAGMAT-5 (Supplementary Data Electronic Appendix 2, Part 4). To test the validity of the phase equilibria calculations, we selected results of eight 1 atm. synthesis melting experiments carried out on the 85–41c primitive magnesian andesite from Mt. Shasta (Grove et al., 2003), see Table 2. Two simulations of the equilibrium crystallization for this melt were carried out. The first was done using the basic version of the COMAGMAT-5, whereas the second modeling was performed using a slight correction of the calculated temperatures (see captions to Fig. 2a). In addition, we used a previous version of COMAGMAT-3.n to calculate the crystallization trajectory for the same melt. All simulations were carried out at 1 atm. and QFM, similar to conditions of the experimental study (Grove et al., 2003). Fig. 2. View largeDownload slide Experimental and modeled sequences for equilibrium crystallization of the high-Mg andesite 85–41c (Grove et al., 2003) at P = 1 atm and QFM. The dashed lines represent experimental liquidus temperatures for olivine, orthopyroxene, and clinopyroxene & plagioclase. (A) Calculations using the COMAGMAT-5·2 model were conducted without corrections for the modeled mineral temperatures and using small temperature shifts with respect to the original calculated values (−10°C for olivine, -12°C for plagioclase, -3°C for high-Ca pyroxene, and +5°C for orthopyroxene), see option 4·5 in Supplementary Data Electronic Appendix 2 (Part 4). For comparison, results of simulations by COMAGMAT-3·72 are shown. Both models were downloaded from http://geo.web.ru/∼ariskin/soft.html-id=comagmat.htm. All crystallization trajectories were constructed with crystallization increments from 0·5% up to 80% crystallized. The composition of the HMA is listed in Table 2. (B) MELTS calculations (Ghiorso & Sack, 1995) were carried out using the Adiabat_1ph version 1·8 run by an EXCEL-based front-end using a shell command (developed by J.-I. Kimura, JAMSTEC, Japan), whereas pMELTS (Ghiorso et al., 2002) and RMELTS (Gualda et al., 2012) codes were downloaded from http://melts.ofm-research.org/unix.html. The MELTS crystallization trajectory was modeled with a 10°C temperature increment in the 1300–1000°C range; calculations using pMELTS and RMELTS were carried out in the 1265–1000°C range with a 5°C temperature increment. The calculations using pMELTS and RMELTS resulted in a wide field of spinel (with onset of crystallization temperatures of 1160°C and 1220°C, correspondingly), not shown in the Figure 2b plots. A complete series of EXCEL-files including results of the modeling is available upon request from the authors. Fig. 2. View largeDownload slide Experimental and modeled sequences for equilibrium crystallization of the high-Mg andesite 85–41c (Grove et al., 2003) at P = 1 atm and QFM. The dashed lines represent experimental liquidus temperatures for olivine, orthopyroxene, and clinopyroxene & plagioclase. (A) Calculations using the COMAGMAT-5·2 model were conducted without corrections for the modeled mineral temperatures and using small temperature shifts with respect to the original calculated values (−10°C for olivine, -12°C for plagioclase, -3°C for high-Ca pyroxene, and +5°C for orthopyroxene), see option 4·5 in Supplementary Data Electronic Appendix 2 (Part 4). For comparison, results of simulations by COMAGMAT-3·72 are shown. Both models were downloaded from http://geo.web.ru/∼ariskin/soft.html-id=comagmat.htm. All crystallization trajectories were constructed with crystallization increments from 0·5% up to 80% crystallized. The composition of the HMA is listed in Table 2. (B) MELTS calculations (Ghiorso & Sack, 1995) were carried out using the Adiabat_1ph version 1·8 run by an EXCEL-based front-end using a shell command (developed by J.-I. Kimura, JAMSTEC, Japan), whereas pMELTS (Ghiorso et al., 2002) and RMELTS (Gualda et al., 2012) codes were downloaded from http://melts.ofm-research.org/unix.html. The MELTS crystallization trajectory was modeled with a 10°C temperature increment in the 1300–1000°C range; calculations using pMELTS and RMELTS were carried out in the 1265–1000°C range with a 5°C temperature increment. The calculations using pMELTS and RMELTS resulted in a wide field of spinel (with onset of crystallization temperatures of 1160°C and 1220°C, correspondingly), not shown in the Figure 2b plots. A complete series of EXCEL-files including results of the modeling is available upon request from the authors. Comparisons of the COMAGMAT calculation results with those obtained in experiments are shown in Figure 2a. These data demonstrate that the basic version of COMAGMAT-5·2 is able to predict the major minerals crystallizing in the experimental range; however, it overestimates by ∼20°C the liquidus temperature for olivine and underestimates by ∼20°C that for orthopyroxene. The modeled temperatures for plagioclase and augite are accurate, if the range of 15°C between two experimental points where these phases are absent and present (Grove et al., 2003) is taken into account. At temperatures below 1120°C, COMAGMAT-5 predicts a peritectic reaction between orthopyroxene and pigeonite and the occurrence of ilmenite. COMAGMAT-3·72 demonstrates a similar sequence of crystallization; however, the modeled olivine liquidus temperature is slightly higher as compared to COMAGMAT-5·2, whereas the modeled onset of orthopyroxene crystallization is delayed by 15°C (Fig. 2a). Examples of calculations with the corrected mineral temperatures are given in Figure 2a (see plot labeled ‘COMAGMAT 5·2 corrected’) to demonstrate the COMAGMAT-5 capability for ‘fine tuning’ the model, in order to make the modeling more consistent with available experimental data or natural observations. It is possible to see that even small shifts of the modeled temperatures around ±10°C or less allow construction of the crystallization order exactly as observed in the experiments (Grove et al., 2003): Ol → Ol + Opx → Opx only → Opx + Pl + Aug. Note that it would be difficult (or perhaps impossible) to do this using crystallization programs, which calculate phase equilibria at predetermined temperatures. Figure 2b represents similar comparisons of the thermodynamic calculations on the same HMA (high magnesium andesite) composition, using three models of the MELTS family, including the original MELTS (Ghiorso & Sack, 1995), pMELTS (Ghiorso et al., 2002), and RMELTS (Gualda et al., 2012). Similar to COMAGMAT-3·72 and COMAGMAT-5·2, both MELTS and RMELTS are accurate enough to predict the experimental crystallization sequence, whereas pMELTS strongly underestimates crystallization temperatures for clinopyroxene and, particularly, orthopyroxene. Note that the pMELTS calculations were carried out at pressures far from the recommended 1–4 GPa range of pressures (Ghiorso et al., 2002). Generally, both series of plots in Figure 2a and b demonstrate that modern magma crystallization models should have a user-friendly capability to slightly change or to correct input thermodynamic parameters, allowing petrologists to make the results of their calculations more consistent with results obtained using experimental data. Data for the corrected COMAGMAT-5 (see captions to Fig. 2a) show that even insignificant temperature shifts using option 4·5 in COMAGMAT-5·2 (Supplementary Data Electronic Appendix 2, Part 4) may produce essential changes in the relative volumes of crystallized pyroxenes and olivine. The reasons for the strong effect of the modeled temperatures are discussed in Supplementary Data Electronic Appendix 2 (Part 5). APPLICATION OF COMAGMAT-5 TO THE BUSHVELD COMPLEX Below two examples of sulfide COMAGMAT use are given which demonstrate the genetic potential of the model, when it is applied to the sulfide saturation history of large layered intrusions and their parental magmas. We focus on modeling sulfide immiscibility in initially sulfide-undersaturated systems; these systems approximate the most primitive magmas and cumulates from the Basal Ultramafic Sequence of the Bushveld Complex in South Africa (Wilson, 2015). In addition, essential differences in the sulfide saturation history during equilibrium and fractional crystallization of the same ultramafic (UM) magma are considered. The Bushveld Complex and its parental magmas The Bushveld Complex was formed within a stable cratonic shield ∼2·06 Ga ago, as a result of three distinctive events; it covers an area up to 120 000 km2, including intrusions of mafic/ultramafic to granitic rocks into the Transvaal sedimentary sequences of the Kaapvaal craton (Naldrett et al., 2012). The first suite of rocks consists of early mafic sills, followed by felsic volcanic rocks forming the roof of a 7–8 km thick mafic to ultramafic cumulate succession, which represents the second intrusive suite, broadly referred to as the ‘Bushveld Complex’ or, more specifically, the so-called ‘Rustenburg Layered Suite’ (RLS; South African Committee for Stratigraphy, 1980). The radially inward dipping layered rocks cover an area of 40 000 km2 or even more (Cawthorn, 2015), thus comprising the world’s largest mafic/ultramafic intrusion with giant Ni–Cu–PGE, Cr, and Ti–V–Fe deposits (Eales & Costin, 2012). The third component of the Bushveld intrusive event is the Lebowa Granite Suite, which is possibly 2 km thick (Cawthorn, 2015). The detailed structure of the RLS has been examinated in five peripheral areas of the Bushveld Complex, including three (initially probably interconnected) basin-like lobes, the Eastern Bushveld, the Southeastern Bushveld and the Western Bushveld, as well as the Far Western Bushveld and the Northern Limb. A classical stratigraphic succession through the RLS succession is divided into five major zones, including the Marginal Zone composed of norite and minor pyroxenite (up to 800 m thick), the 800–1300 m thick ultramafic Lower Zone (mostly orthopyroxenite, harzburgite, and minor dunite), the Critical Zone marked by incoming chromitites in addition to orthopyroxenite (1·3–1·8 km thick), which is overlain by the ∼3·2 km thick Main Zone (norite, gabbros and minor anorthosite) and an as much as 2·8 km thick Upper Zone, composed of ferrogabbro and ferrodiorite (Naldrett et al., 2012; Cawthorn, 2015). During recent years, important progress has been made in constraining the compositional parameters of the most primitive Bushveld Complex magmas; see discussions in Wilson (2012, 2015), Yudovskaya et al. (2013, 2015) and Maier et al. (2016). Instead of the so-called B1 magma (12–13·5 wt % MgO, ∼56 wt % SiO2), which has long been considered as one probable parent for the Lower Zone and the Lower Critical Zone (Barnes et al., 2010; Wilson, 2012), much more primitive parental magmas, containing 19 wt % MgO (Wilson, 2015) and 18·7 wt% MgO (Maier et al., 2016) have been proposed. This is due to the inability of the B1 magma: (i) to produce the thick sequences of dunites and harzburgites in the Lower Zone, and (ii) a strong misbalance of Cr2O3 between the B1 composition and the extraordinarily high average Cr concentrations in the Bushveld rock sequences (Eales & Costin, 2012). Another contradiction is that the B1 magma was not capable of crystallizing olivine with Mg#>91 as observed in the Lower Zone of the Northern Limb and in the chilled sequence of the eastern Bushveld Complex (Wilson, 2012; Yudovskaya et al., 2013). These observations resulted in the conclusion that the B1 magma appears to be derived from a more mafic komatiitic parent; the composition of this magma can be approximated by the UM liquid from the Schwerin fold chill in the eastern Bushveld Complex, see Wilson (2015). Two model compositions were used in the example calculations (Table 2). The first one (denoted as ‘UM’) represents the most primitive UM liquid approximated by the Schwerin fold chill from the Basal Ultramafic Sequence of the Bushveld Complex, see the CH12/7&14UM column in table 4 of Wilson (2015). Focusing on S immiscibility calculations instead of modeling spinel crystallization, we postulated that the Cr2O3 concentration was zero and S concentration was 368 ppm in the UM melt. In fact, Wilson (2015) gave a much higher value 1914 ± 200 ppm S (based on analytical data) in the parent magma, whereas COMAGMAT-5·2 calculated SCSS = 1600 ppm at the UM liquidus temperature of 1456°C. This means that the ultramafic chill is locally over-saturated with sulfides, which contradicts previous conclusions on the sulfide under-saturated character of the B1 magma (Barnes et al., 2010; Ariskin et al., 2013), as a probable crystallization product of the more primitive UM magma (Wilson, 2015). Thus, the reason to reduce the initial S concentration in the initial melt is to avoid over-saturation of the UM magma with sulfide and sulfide immiscibility occurring too early in the modeled system. The second composition (denoted as ‘UM + 56 wt % Ol-91·6’, column 4 in Table 2) represents an approximation of a primitive cumulate rock, which is very similar to the olivine pyroxenite listed in the 5th column of Table 2. This composition is calculated as a mixture of 44 wt % of the proposed UM melt (the UM composition in column 3) with 56 wt% of olivine mg#91·6. We consider mg#91·6 to be a valid approximation of the original olivine composition in the Bushveld cumulates, because petrological reconstructions of the Lower Zone of the Northern Limb, the chilled sequence of the eastern Bushveld Complex and the Basal Ultramafic Sequence of the Bushveld Complex, indicate that the most primitive olivine varies in the range of mg#91–92 (Wilson, 2012, 2015; Yudovskaya et al., 2013). Note that the COMAGMAT-5 calculations for the UM magma yield the same liquidus olivine composition: mg#91·6 (see below). Sulfide immiscibility in the ultramafic magma and related cumulate Two series of plots in Figure 3 display the sulfide saturation history during equilibrium crystallization of two related systems, including the UM magma (Fig. 3a) and the proposed olivine enriched cumulate mush ‘UM + 56 wt % Ol-91·6’ (Fig. 3b). Both crystallization trajectories were calculated at QFM + 0·5 and anhydrous conditions, assuming P = 1 atm. The purpose of these calculations was not to replicate a slightly different order of crystallization at a more reliable pressure of 1–2 kbar (Wilson, 2015; Yudovskaya et al., 2015). Instead, we focused on marked differences in the onset of sulfide immiscibility in these samples at 1 atm., despite the fact that both systems were initially S-undersaturated, having almost the same initial melt (including S concentration) and olivine mg#91·6 at very similar temperatures (1456–1465°C), see Figure 3. Fig. 3. View largeDownload slide Modeled sequences of equilibrium crystallization, proportions of the modeled sulfide liquid with respect to precipitated solids, and compositional characteristics of crystallizing melts for two UM systems approximating a parental UM magma (A) and a related Ol-rich cumulate (B). Compositions of the UM parent (A) and the proposed ‘UM+56 wt.% Ol-91·6’ olivine cumulate (B) are listed in Table 2. Calculations using the COMAGMAT-5·2 model were carried out at QFM+0·5, anhydrous conditions and P=1 atm. The dashed green line in Figure 3b denotes the initial magma temperature consistent with ∼19 wt % MgO in the melt, which is in equilibrium with olivine mg#91·6 (Table 3). Light pink lines represent the evolution of S concentration in the modeled melts. Thick pink lines correspond to the sulfide liquid saturation of the modeled melt, when the sulfide solubility SCSS is equal to the S concentration in the melts. Note, that high NiO concentrations in the ‘UM+56 wt % Ol-91·6’ virtual melts at low percentage of crystallization are due to a large mass of Ni-containing olivine initially added to the UM parental melt. Fig. 3. View largeDownload slide Modeled sequences of equilibrium crystallization, proportions of the modeled sulfide liquid with respect to precipitated solids, and compositional characteristics of crystallizing melts for two UM systems approximating a parental UM magma (A) and a related Ol-rich cumulate (B). Compositions of the UM parent (A) and the proposed ‘UM+56 wt.% Ol-91·6’ olivine cumulate (B) are listed in Table 2. Calculations using the COMAGMAT-5·2 model were carried out at QFM+0·5, anhydrous conditions and P=1 atm. The dashed green line in Figure 3b denotes the initial magma temperature consistent with ∼19 wt % MgO in the melt, which is in equilibrium with olivine mg#91·6 (Table 3). Light pink lines represent the evolution of S concentration in the modeled melts. Thick pink lines correspond to the sulfide liquid saturation of the modeled melt, when the sulfide solubility SCSS is equal to the S concentration in the melts. Note, that high NiO concentrations in the ‘UM+56 wt % Ol-91·6’ virtual melts at low percentage of crystallization are due to a large mass of Ni-containing olivine initially added to the UM parental melt. Modeling the UM crystallization yields: Ol (1456°C, mg#91·6) → Ol (1222°C, mg#84·7) + Opx (mg#87·0) → Ol (1185°C, mg#82·7) + Opx (mg#85·2) + Sulfide liquid → Ol (1178°C, mg#82·3) + Opx (mg#84·9) + Sulfide liquid + Pl (An72·7) → Ol (1152°C, mg#78·6) + Opx (mg#82·2) + Sulfide liquid + Pl (An66·5) + Aug (mg#82·4) → plus ilmenite (at 1116°C) and Ti–magnetite (at 1105°C). Changes in the crystallizing mineral assemblages are shown in Figure 3a in terms of the percent crystallized. The ‘UM + 56 wt % Ol-91·6’ cumulate melt yields a similar order of crystallization, only starting at much higher temperatures: Ol (1724°C, mg#95·8) → Ol (1289°C, mg#89·5) + sulfide liquid → Ol (1231°C, mg#89·0) + Opx (mg#90·6) + sulfide liquid → Ol (1194°C, mg#88·6) + Opx (mg#90·1) + sulfide liquid + Pl (An73·8) → Ol (1165°C, mg#87·8) + Opx (mg#89·7) + sulfide liquid + Pl (An65·8) + Aug (mg#89·7) → plus ilmenite (at 1134°C and lower), see Figure 3b. In the following interpretation, it is noteworthy that neither the calculated maximum liquidus temperature (1724°C), nor the modeled liquidus olivine (mg#95·8), represent realistic estimates of initial magma parameters. Instead, the results of calculations at 1724°C should be considered as virtual characteristics of the 100% melted olivine rich ‘model cumulate’. The same is true for other calculations in the fictitious range 1724°C to ∼1465°C, with the latter value corresponding to the proposed two-phase mixture of a magmatic melt (19 wt % MgO) and olivine mg#91·6. Thus, a petrological interpretation only considers results obtained below the proposed UM magma temperature (approximately ≤1460°C); see Figure 3a and the dashed line in Figure 3b. The main differences between the ‘UM’ and ‘UM + 56 wt % Ol-91·6’ calculations within the apparent range of magma crystallization temperatures (∼1460–1120°C) are recorded in the onset of modeled sulfide immiscibility. The first sulfide liquid was calculated to occur at 1185°C in the case of the parental UM magma and at the much higher 1289°C temperature for the modeled olivine cumulate. We consider these differences as one more demonstration of the coupled compositional effect of FeO and Ni on sulfide solubility (SCSS), because the olivine cumulate contained approximately three times more NiO than the pure UM melt (Table 3), see details in Discussion. Table 3: Parameters of equilibrium and fractional crystallization for the proposed UM magma and modeled cumulate ‘UM + 56 wt.% Ol-91·6’, see Table 2 Melt components, wt.% Modeled Ol cumulate ‘UM + 56 wt.% Ol-91·6’ Proposed UM magma Crystallization Equilibrium Equilibrium Fractional Bulk NiO in the modeled systems  NiO, wt.% 0·317 0·090 Parameters of magmatic melt in equilibrium with olivine mg#91·6  Temp, °C 1465 1456  MgO, wt.% 19·4 18·9  FeO, wt.% 10·39 10·14  NiO, wt.% 0·114 0·090  Sinit, wt.% 0·036 0·037  SCSS, wt.% 0·163 0·160  mg# in Ol,% 91·6 91·6  NiO in Ol, wt. % 0·479 0·393  Ol, wt.% 44·5 0 The onset of sulfide immiscibility Mineral assemblage Ol+Sulfide liquid Ol+OPx+Sulfide liquid Pl+Pig+Sulfide liquid  Temp, °C 1289 1185 1116  MgOmelt, wt.% 10·22 6·66 2·93  FeO, wt.% 7·08 8·34 11·40  NiOmelt, wt.% 0·039 0·015 0·002  S (in melt) = SCSS, wt.% 0·050 0·062 0·086  mg# in Ol, % 89·5 82·7 –  NiO in Ol, wt.% 0·444 0·290 –  Ni/(Ni+Fe) in sulfide, atomic 0·410 0·318 0·237 Melt components, wt.% Modeled Ol cumulate ‘UM + 56 wt.% Ol-91·6’ Proposed UM magma Crystallization Equilibrium Equilibrium Fractional Bulk NiO in the modeled systems  NiO, wt.% 0·317 0·090 Parameters of magmatic melt in equilibrium with olivine mg#91·6  Temp, °C 1465 1456  MgO, wt.% 19·4 18·9  FeO, wt.% 10·39 10·14  NiO, wt.% 0·114 0·090  Sinit, wt.% 0·036 0·037  SCSS, wt.% 0·163 0·160  mg# in Ol,% 91·6 91·6  NiO in Ol, wt. % 0·479 0·393  Ol, wt.% 44·5 0 The onset of sulfide immiscibility Mineral assemblage Ol+Sulfide liquid Ol+OPx+Sulfide liquid Pl+Pig+Sulfide liquid  Temp, °C 1289 1185 1116  MgOmelt, wt.% 10·22 6·66 2·93  FeO, wt.% 7·08 8·34 11·40  NiOmelt, wt.% 0·039 0·015 0·002  S (in melt) = SCSS, wt.% 0·050 0·062 0·086  mg# in Ol, % 89·5 82·7 –  NiO in Ol, wt.% 0·444 0·290 –  Ni/(Ni+Fe) in sulfide, atomic 0·410 0·318 0·237 Table 3: Parameters of equilibrium and fractional crystallization for the proposed UM magma and modeled cumulate ‘UM + 56 wt.% Ol-91·6’, see Table 2 Melt components, wt.% Modeled Ol cumulate ‘UM + 56 wt.% Ol-91·6’ Proposed UM magma Crystallization Equilibrium Equilibrium Fractional Bulk NiO in the modeled systems  NiO, wt.% 0·317 0·090 Parameters of magmatic melt in equilibrium with olivine mg#91·6  Temp, °C 1465 1456  MgO, wt.% 19·4 18·9  FeO, wt.% 10·39 10·14  NiO, wt.% 0·114 0·090  Sinit, wt.% 0·036 0·037  SCSS, wt.% 0·163 0·160  mg# in Ol,% 91·6 91·6  NiO in Ol, wt. % 0·479 0·393  Ol, wt.% 44·5 0 The onset of sulfide immiscibility Mineral assemblage Ol+Sulfide liquid Ol+OPx+Sulfide liquid Pl+Pig+Sulfide liquid  Temp, °C 1289 1185 1116  MgOmelt, wt.% 10·22 6·66 2·93  FeO, wt.% 7·08 8·34 11·40  NiOmelt, wt.% 0·039 0·015 0·002  S (in melt) = SCSS, wt.% 0·050 0·062 0·086  mg# in Ol, % 89·5 82·7 –  NiO in Ol, wt.% 0·444 0·290 –  Ni/(Ni+Fe) in sulfide, atomic 0·410 0·318 0·237 Melt components, wt.% Modeled Ol cumulate ‘UM + 56 wt.% Ol-91·6’ Proposed UM magma Crystallization Equilibrium Equilibrium Fractional Bulk NiO in the modeled systems  NiO, wt.% 0·317 0·090 Parameters of magmatic melt in equilibrium with olivine mg#91·6  Temp, °C 1465 1456  MgO, wt.% 19·4 18·9  FeO, wt.% 10·39 10·14  NiO, wt.% 0·114 0·090  Sinit, wt.% 0·036 0·037  SCSS, wt.% 0·163 0·160  mg# in Ol,% 91·6 91·6  NiO in Ol, wt. % 0·479 0·393  Ol, wt.% 44·5 0 The onset of sulfide immiscibility Mineral assemblage Ol+Sulfide liquid Ol+OPx+Sulfide liquid Pl+Pig+Sulfide liquid  Temp, °C 1289 1185 1116  MgOmelt, wt.% 10·22 6·66 2·93  FeO, wt.% 7·08 8·34 11·40  NiOmelt, wt.% 0·039 0·015 0·002  S (in melt) = SCSS, wt.% 0·050 0·062 0·086  mg# in Ol, % 89·5 82·7 –  NiO in Ol, wt.% 0·444 0·290 –  Ni/(Ni+Fe) in sulfide, atomic 0·410 0·318 0·237 Fractional vs. equilibrium crystallization The crystallization sequence shown in Figure 4 represents the results of modeling perfect fractionation of the same UM magma shown in Figure 3a. During the early stages in the temperature range ∼1460–1200°C, modeling equilibrium (Fig. 3a) and fractional (Fig. 4) crystallization results in the same order of the two first minerals to crystallize, olivine → orthopyroxene, but the late stages differ considerably. First, note the much later sulfide immiscibility, which is observed at ∼70°C lower than the temperature of sulfide precipitation onset during equilibrium crystallization (Table 3). As far as orthopyroxene during fractionation was replaced by pigeonite, sulfide saturation is contemporaneous with the crystallization of a Pl + Pig gabbroic assemblage, followed by Pl + Pig + Aug + Mt (Fig. 4). Other sulfide saturation parameters during UM parental magma fractional crystallization are given in Table 3. Fig. 4. View largeDownload slide The order of fractional crystallization of the UM parental magma (Table 2) and comparison of S vs. SCSS trends and NiO in the melt in the case of fractional and equilibrium crystallization. Calculations by the COMAGMAT-5·2 model were carried at QFM+0·5, anhydrous conditions, and P=1 atm. The solid gray line and the solid black line represent changes in S and NiO contents in the melt in the case of equilibrium and fractional crystallization, respectively. The dotted lines display the evolution of the SCSS values for both crystallization trajectories. The onset of sulfide-silicate immiscibility is consistent with the intersection of the S and SCSS trends for each type of crystallization. The light and bold pink lines denote intervals of crystallization under sulfide saturated conditions. Fig. 4. View largeDownload slide The order of fractional crystallization of the UM parental magma (Table 2) and comparison of S vs. SCSS trends and NiO in the melt in the case of fractional and equilibrium crystallization. Calculations by the COMAGMAT-5·2 model were carried at QFM+0·5, anhydrous conditions, and P=1 atm. The solid gray line and the solid black line represent changes in S and NiO contents in the melt in the case of equilibrium and fractional crystallization, respectively. The dotted lines display the evolution of the SCSS values for both crystallization trajectories. The onset of sulfide-silicate immiscibility is consistent with the intersection of the S and SCSS trends for each type of crystallization. The light and bold pink lines denote intervals of crystallization under sulfide saturated conditions. Two pairs of curves to the right in Figure 4 display changes in the modeled S concentrations and SCSS in the melt, allowing insights into the differences in the onset of sulfide–silicate immiscibility between fractional and equilibrium crystallization. To interpret the modeled relationships, remember that our approach to simulating sulfide saturation is similar to that proposed by Li & Ripley (2005) and utilized by other workers (e.g. Barnes, 2007; Barnes et al., 2010; Ariskin et al., 2016). The evolution of melt S content in the case of equilibrium (the solid grey line) and fractional (the solid black line) crystallization was calculated assuming an initial S content of 368 ppm in the UM parent (Table 2). The sulfide solubility values (SCSS, dotted lines) were calculated as a function of varying temperature and melt composition, including Ni concentrations in the derivative melts (Ariskin et al., 2013). The point of S and SCSS lines intersection (i.e. S = SCSS) shows the onset of sulfide-silicate immiscibility, thus giving rise to precipitation of a Fe–Ni sulfide liquid, so that the S content in the melt follows the SCSS trend (Fig. 4). The modeled relationships demonstrate that higher sulfide solubility in the case of the fractional crystallization of the UM parent (compare SCSS values at the same temperature) should be correlated with an increase in FeO contents and much lower NiO in the melt, due to fractionation of olivine and low-Ca pyroxenes (Opx/Pig). Finally, this gives rise to much later sulfide saturation in the modeled system as compared to the course of equilibrium crystallization of the same magmatic melt (compare columns 2 & 3 in Table 3). DISCUSSION Problems of thermodynamics and the use of COMAGMAT-5 G-minimization vs. the mass action law As is evident from the thermodynamic considerations discussed here (see Eqs. (1–12)), despite its semi-empirical character, COMAGMAT was designed to follow fundamental thermodynamic principles. The main difference between its basic algorithm and those used in MELTS/pMELTS is that the algorithms used in the latter models directly minimize the isobaric potential at a given temperature (Ghiorso, 1987, 1994), whereas design of the COMAGMAT models allows the equilibrium state to be searched for at a pre-defined crystallization degree. In this case, the equilibrium problem is approached by means of the iterative solution of nonlinear equations consistent with the law of mass action (Eqs. (4, 5)). It is noteworthy that methods of calculating the products of chemical reactions based on the mass action law were developed in chemistry long before Gibbs formulated his thermodynamic theory. Thus, there is no principal conflict between these two methods, because the constants of heterogeneous equilibria are linked to the changes in the free energy through the fundamental equation ln K=-ΔG/RgT. In this sense, it would be possible to transform the energy functions used by Ghiorso and his followers to the ‘language of equilibrium constants’ describing mineral–melt equilibria. Their values (more accurately lnK vs. 1/T dependences) could then be used in the COMAGMAT algorithm to solve the equilibrium problem in the same manner as is done using the empirically calibrated mineral–melt geothermometers. Approaching the effect of low pressures Because all mineral–melt geothermometers used in the COMAGMAT-5 model were calibrated on experimental data obtained at 1 atm., formally, this model should be used only for ‘atmospheric’ calculations. However, based on our experience in modeling crystallization at elevated pressures (Ariskin et al., 1990; Ariskin & Barmina, 2004), we have slightly extended the possible range of pressures, up to about 1–2 kbar. There are two reasons for such approximations. First, the standard deviations of olivine and plagioclase thermometric calculations (±10–15°C) are two to three times higher than the effect of pressure on their liquidus temperatures (∼5°C/kbar). For pyroxenes this effect is known to be ∼10–12°C/kbar, which is also within the accuracy of the proposed pyroxene–melt geothermometers. Second, an increase in total pressure results in a very similar increase in the liquidus temperatures for olivine and plagioclase; therefore, the chemical evolution of the modeled melts crystallizing at slightly elevated pressures barely deviates from that obtained at 1 atm., at least up to the onset of pyroxene crystallization (Ariskin et al., 1990). Admitedly in the case of modeling tholeiitic or other systems, which crystallize olivine and plagioclase prior to crystallizing clinopyroxene, the effect of small pressure changes on the modeled crystallization of troctolitic assemblages should be negligible. This supports the application of COMAGMAT-5 to shallow magma chambers occurring at depths of no more 6–7 km, particularly for Ol±Pl saturated magmas similar to MORB (Ariskin et al., 2016). The effect of Ni on SCSS Using the sulfide COMAGMAT-5, it was shown (Ariskin et al., 2013) that despite relatively low concentrations, Ni has a pronounced lowering effect on S solubility, causing a significant temperature increase in the onset of sulfide immiscibility in melts with similar major element compositions. Despite the fact that this hypothesis has found some support in natural observations (e.g. Barnes et al., 2013), it continues to be debated. Recently Fortin et al. (2015) presented a new FeS–sulfide solubility equation (his so-called ‘Model-B’), which reproduced S contents both for their calibration dataset and for 65 natural Ni-containing sulfide-saturated MORB glasses within about 5%, even without considering the effect of Ni. It was concluded that ‘there is no need for a more complicated multi-species FeNiS model to calculate the SCSS, even in high-Ni, natural samples’ (Fortin et al., 2015, p. 113). In our study, we present an additional test of the Ni-effect on a magnesian Ni-rich sulfide-saturated glass dredged from the the Southern Atlantic. These results clearly demonstrate the strong effect of Ni on the modeled sulfide solubility (Table 1). This outcome seems to be of fundamental importance, particularly when the SCSS calculations are applied to the history of sulfide saturation in picritic magmas and mafic to ultramafic layered intrusions. This is because of the higher Ni contents in the high-temperature melt and the possible presence of a large amount of Ni-rich olivine in the parental magmas, so that the bulk composition of the ultramafic system has to control the sulfide immiscibility process during both magma crystallization and solidification of intercumulus melt in igneous cumulates (Table 3). CONCLUSIONS A detailed description of the thermodynamic framework and main options of the first ‘sulfide version’ of the COMAGMAT-5 magma crystallization model is presented. The new features include: (i) complete recalibration of silicate mineral– melt geothermometers; (ii) the capability to calculate the Ni content of the modeled femic minerals and sulfides; and (iii) incorporation of a recently proposed Fe–Ni sulfide solubility model (Ariskin et al., 2013) into a general algorithm of modeling crystallization for S-containing magmas. Due to updating of the basic COMAGMAT algorithm, the new model allows new geothermometers to be used for simulations of co-crystallization Aug±Pig±Opx, including peritectic relations between these phases and olivine. Combining these calculations with the proposed SCSS model, COMAGMAT-5 can be used to simulate equilibrium and fractional crystallization of S-saturated and S-undersaturated mafic to UM magmas, including changes in sulfide-silicate (±Fe–Ti oxides) proportions for multiply-saturated mineral assemblages. The COMAGMAT-5 model (ver. 5·2) is available for downloading from http://geo.web.ru/∼ariskin/soft.html-id=comagmat.htm. Following Ariskin et al. (2013), the strong effect of Ni on SCSS modeled in mafic to ultramafic systems must be taken into account, so that at higher Ni contents in the melt an essential decrease in the sulfide solubility should be expected. This statement was tested against a Ni-rich sulfide-saturated glass dredged from the southern Mid-Atlantic Ridge near the Bouvet Triple Junction (Kamenetsky et al., 2001, 2013). Applying COMAGMAT-5 to the BTJ glass, we obtained SCSS from 514 ppm (QFM) to 710 ppm (QFM + 1), close to the 600 ppm observed (Kamenetsky et al., 2013). The use of virtual ‘FeS’ solubility models results in much higher sulfide solubility, mostly > 1000 ppm (Li & Ripley, 2009; Fortin et al., 2015), even assuming the strong dampening effect of low H2O contents on SCSS postulated in Liu et al. (2007). Example calculations using compositional proxies of a UM parental melt and a related olivine cumulate from the Basal Ultramafic Sequence of the Bushveld Complex (see UM in Table 2; Wilson, 2015) demonstrate a high temperature interval in the modeled trajectories where sulfide–silicate immiscibility could occur. The highest 1289°C temperature (Ol + sulfide) was observed in crystallizing olivine cumulate, an intermediate 1185°C temperature (Ol + Opx + sulfide) was calculated in the case of equilibrium UM magma crystallization, and the lowest 1116°C temperature (Pl + Pig + sulfide) was obtained as a result of perfect fractionation of the same melt (Table 3). The modeled effect seems to be of general genetic significance, particularly for the study of the formation of sulfide-saturated volcanic suites and fertile layered intrusions. ACKNOWLEDGEMENTS We gratefully acknowledge thoughtful reviews by Sébastien Jégo, Guil Gualda, Gokce Ustunisik, Paula Antoshechkina and two anonymous reviewers, as well as constructive comments from editor Wendy Bohrson. We thank Leonid Danyushevsky for inspiring discussions, Jun-Ichi Kimura for his help with MELTS/pMELTS calculations, and Candace S. O'Connor for careful editing of the final version of the manuscript. FUNDING We acknowledge support of AngloAmerican, BHP Billiton, Votorantim Metais (through AMIRA project P962, 2007–2010) and the Australian Research Council funding to the CODES (Hobart, Australia) at the initial stage of development of the COMAGMAT model. Final calibrations, testing and example calculations using COMAGMAT-5 are supported by the Russian Science Foundation (project RSF 16–17-10129). SUPPLEMENTARY DATA Supplementary data for this paper are available at Journal of Petrology online. REFERENCES Almeev R. R. , Ariskin A. A. , Kimura J.-I. , Barmina G. S. ( 2013 ). The role of polybaric crystallization in genesis of andesitic magmas: phase equilibria simulations of Bezymianny volcanic subseries . Journal of Volcanology and Geothermal Research 263 , 182 – 192 . Google Scholar CrossRef Search ADS Ariskin A. A. ( 1999 ). Phase equilibria modeling in igneous petrology: use of COMAGMAT model for simulating fractionation of ferro-basaltic magmas and the genesis of high-alumina basalt . Journal of Volcanology and Geothermal Researches 90 , 115 – 162 . Google Scholar CrossRef Search ADS Ariskin A. A. ( 2003 ). The compositional evolution of differentiated liquids from the Skaergaard layered series as determined by geochemical thermometry . Russian Journal of Earth Sciences 5 , 1 – 29 . Google Scholar CrossRef Search ADS Ariskin A. A. , Barmina G. S. ( 1990 ). Equilibria thermometry between plagioclases and basalt or andesite magmas . Geochemistry International 27 , 129 – 134 . Ariskin A. A. , Barmina G. S. ( 2004 ). COMAGMAT: development of a magma crystallization model and its petrologic applications . Geochemistry International 42 (Suppl 1) , S1 – S157 . Ariskin A. A. , Yaroshevsky A. A. ( 2006 ). Crystallization differentiation of intrusive magmatic melt: development of a Convection-Accumulation model . Geochemistry International 44 , 72 – 93 . Google Scholar CrossRef Search ADS Ariskin A. A. , Frenkel M. Y. , Tsekhonya T. I. ( 1990 ). High-pressure fractional crystallization of tholeiitic magmas . Geochemistry International 27 , 10 – 20 . Ariskin A. A. , Barmina G. S. , Bychkov K. A. , Danyushevsky L. V. ( 2009 ). Parental magmas of mafic layered intrusions: using an updated COMAGMAT model for calculations of sulfide-silicate cotectics during their crystallization . Northwestern Geology 42 , 1 – 3 . Ariskin A. A. , Barmina G. S. , Frenkel M. Ya. , Nielsen R. L. ( 1993 ). COMAGMAT: a Fortran program to model magma differentiation processes . Computers and Geosciences 19 , 1155 – 1170 . Google Scholar CrossRef Search ADS Ariskin A. A. , Danyushevsky L. V. , Bychkov K. A. , McNeill A. W. , Barmina G. S. , Nikolaev G. S. ( 2013 ). Modeling solubility of Fe-Ni sulfides in basaltic magmas: the effect of Ni in the melt . Economic Geology 108 , 1983 – 2003 . Google Scholar CrossRef Search ADS Ariskin A. A. , Kislov E. V. , Danyushevsky L. V. , Nikolaev G. S. , Fiorentini M. L. , Gilbert S. , Goemann K. , Malyshev A. ( 2016 ). Cu-Ni-PGE fertility of the Yoko-Dovyren layered massif (Northern Transbaikalia, Russia): thermodynamic modeling of sulfide compositions in low mineralized dunites based on quantitative sulfide mineralogy . Mineralium Deposita 51 , 993 – 1011 . Google Scholar CrossRef Search ADS Barnes S. J. ( 2007 ). Cotectic precipitation of olivine and sulfide liquid from komatiite magma and the origin of komatiite-hosted disseminated nickel sulfide mineralization at Mount Keith and Yakabindie, Western Australia . Economic Geology 106 , 298 – 304 . Barnes S.-J. , Lightfoot P. C. ( 2005 ). Formation of magmatic nickel sulfide deposits and processes affecting their copper and Platinum Group Element contents . Economic Geology 100th Anniversary Volume, 179 – 213 . Barnes S.-J. , Maier W. D. , Curl E. A. ( 2010 ). Composition of the marginal rocks and sills of the Rustenburg Layered Suite, Bushveld Complex, South Africa: implications for the formation of the Platinum-Group element deposits . Economic Geology 105 , 1491 – 1511 . Google Scholar CrossRef Search ADS Barnes S. J. , Godel B. , Gürer D. , Brenan J. M. , Robertson J. , Paterson D. ( 2013 ). Sulfide-olivine Fe-Ni exchange and the origin of anomalously Ni rich magmatic sulfides . Economic Geology 108 , 1971 – 1982 . Google Scholar CrossRef Search ADS Berman R. G. ( 1988 ). Internally-consistent thermodynamic data for minerals in the system Na2O-K2O-CaO-MgO-FeO-Fe2O3-Al2O3-SiO2-TiO2-H2O-CO2 . Journal of Petrology 29 , 445 – 522 . Google Scholar CrossRef Search ADS Campbell I. H. , Naldrett A. J. ( 1979 ). The influence of silicate: sulfide ratios on the geochemistry of magmatic sulfides . Economic Geology 74 , 1503 – 1506 . Google Scholar CrossRef Search ADS Cawthorn R. G. ( 2015 ). The Bushveld Complex, South Africa. In: Charlier B. , Namur O. , Latypov R. , Tegner C. (eds) Layered Intrusions . Berlin : Springer , pp. 517 – 587 . Google Scholar CrossRef Search ADS Danyushevsky L. V. , Plechov P. ( 2011 ). Petrolog3: integrated software for modeling crystallization processes . Geochemistry, Geophysics, Geosystems 12 , Q07021 . Google Scholar CrossRef Search ADS Eales H. V. , Costin G. ( 2012 ). Crustally contaminated komatiite: primary source of the chromitites and Marginal, Lower, and Critical zone magmas in a staging chamber beneath the Bushveld Complex . Economic Geology 107 , 645 – 665 . Google Scholar CrossRef Search ADS Falloon T. J. , Danyushevsky L. V. , Ariskin A. , Green D. H. , Ford C. E. ( 2007 ). The application of olivine geothermometry to infer crystallization temperatures of parental liquids: Implications for the temperature of MORB magmas . Chemical Geology 241 , 207 – 233 . Google Scholar CrossRef Search ADS Ford C. , Russell D. , Craven J. , Fisk M. ( 1983 ). Olivine-liquid equilibria: temperature, pressure and composition dependence of the crystal/liquid cation partition coefficients for Mg, Fe2+, Ca and Mn . Journal of Petrology 24 , 256 – 265 . Google Scholar CrossRef Search ADS Fortin M.-A. , Riddle J. , Desjardins-Langlais Y. , Baker D. R. ( 2015 ). The effect of water on the sulfur concentration at sulfide saturation (SCSS) in natural melts . Geochimica et Cosmochimica Acta 160 , 100 – 116 . Google Scholar CrossRef Search ADS Frenkel M. Y. , Ariskin A. A. ( 1984 ). A computer algorithm for equilibration in a crystallizing basalt magma . Geochemistry International 21 , 63 – 73 . Frenkel M. Y. , Yaroshevsky A. A. , Ariskin A. A. , Barmina G. S. , Koptev-Dvornikov E. V. , Kireev B. S. ( 1989 ). Convective-cumulative model simulating the formation process of stratified intrusions. In: Bonin B. , Didier J. , Le Fort P. , Propach G. , Puga E. , Vistelius A. B. (eds) Magma-Crust Interactions and Evolution . Athens-Greece : Theophrastus Publications SA , pp. 3 – 88 . Ghiorso M. S. ( 1987 ). Modeling magmatic systems: thermodynamic relations. In: Carmichael I. S. E. , Eugster H. P. (eds) MSA Reviews in Mineralogy 17: Thermodynamic Modeling of Geological Materials: Minerals, Fluids and Melts, 443–465. Ghiorso M. S. ( 1994 ). Algorithms for the estimation of phase stability in heterogeneous thermodynamic systems . Geochimica et Cosmochimica Acta 58 , 5489 – 5501 . Google Scholar CrossRef Search ADS Ghiorso M. S. ( 1997 ). Thermodynamic models of igneous processes . Annual Review of Earth and Planetary Sciences 25 , 221 – 241 . Google Scholar CrossRef Search ADS Ghiorso M. S. , Sack R. O. ( 1995 ). Chemical mass transfer in magmatic processes: IV. A revised and internally consistent thermodynamic model for the interpolation and extrapolation of liquid–solid equilibria in magmatic systems at elevated temperatures and pressures . Contributions to Mineralogy and Petrology 119 , 197 – 212 . Google Scholar CrossRef Search ADS Ghiorso M. S. , Hirschmann M. M. , Reiners P. W. , Kress V. C. III . ( 2002 ). The pMELTS: a revision of MELTS for improved calculation of phase relations and major element partitioning related to partial melting of the mantle to 3 GPa . Geochemistry, Geophysics, Geosystems 3 , doi:10.1029/2001GC000217. Gongalsky B. I. , Krivolutskaya N. A. , Ariskin A. A. , Nikolaev G. S. ( 2016 ). The Chineysky gabbronorite-anorthosite layered massif (Northern Transbaikalia, Russia): its structure, Fe-Ti-V and Cu-PGE deposits, and parental magma composition . Mineralium Deposita 51 , 1013 – 1034 . Google Scholar CrossRef Search ADS Grove T. L. , Elkins-Tanton L. T. , Parman S. W. , Chatterjee N. , Muntener O. , Gaetani G. A. ( 2003 ). Fractional crystallization and mantle-melting controls on calc-alkaline differentiation trends . Contributions to Mineralogy and Petrology 145 , 515 – 533 . Google Scholar CrossRef Search ADS Gualda G. A. R. , Ghiorso M. S. ( 2015 ). MELTS_Excel: a Microsoft Excel-based MELTS interface for research and teaching of magma properties and evolution . Geochemistry, Geophysics, Geosystems 16 , 315 – 324 . Google Scholar CrossRef Search ADS Gualda G. A. R. , Ghiorso M. S. , Lemons R. V. , Carley T. L. ( 2012 ). Rhyolite-MELTS: a modified calibration of MELTS optimized for silicarich, fluid-bearing magmatic systems . Journal of Petrology 53 , 875 – 890 . Google Scholar CrossRef Search ADS Jugo P. J. ( 2009 ). Sulfur content at sulfide saturation in oxidized magmas . Geology 37 , 415 – 418 . Google Scholar CrossRef Search ADS Kamenetsky V. S. , Maas R. , Sushchevskaya N. M. , Norman M. D. , Cartwright I. , Peyve A. A. ( 2001 ). Remnants of Gondwanan continental lithosphere in oceanic upper mantle: evidence from the South Atlantic Ridge . Geology 29 , 243 – 246 . Google Scholar CrossRef Search ADS Kamenetsky V. S. , Maas R. , Fonseca R. O. C. , Ballhaus C. , Heuser A. , Brauns M. , Norman M. D. , Woodhead J. D. , Rodemann T. , Kuzmin D. V. , Bonatti E. ( 2013 ). Noble metals potential of sulfide-saturated melts from the subcontinental lithosphere . Geology 41 , 575 – 578 . Google Scholar CrossRef Search ADS Karpov I. K. , Chudnenko K. V. , Kulik D. A. , Avchenko O. V. , Bychinskii V. A. ( 2001 ). Minimization of Gibbs free energy in geochemical systems by convex programming . Geochemistry International 39 , 1108 – 1119 . Li C. , Ripley E. M. ( 2005 ). Empirical equations to predict the sulfur content of mafic magmas at sulfide saturation and applications to magmatic sulfide deposits . Mineralium Deposita 40 , 218 – 230 . Google Scholar CrossRef Search ADS Li C. , Ripley E. M. ( 2009 ). Sulfur contents at sulfide-liquid or anhydrite saturation in silicate melts: empirical equations and example applications . Economic Geology 104 , 405 – 412 . Google Scholar CrossRef Search ADS Liu Y. , Samaha N.-T. , Baker D. R. ( 2007 ). Sulfur concentration at sulfide saturation (SCSS) in magmatic silicate melts . Geochimica et Cosmochimica Acta 71 , 1783 – 1799 . Google Scholar CrossRef Search ADS Maier W. D. , Barnes S.-J. , Karykowski B. T. ( 2016 ). A chilled margin of komatiite and Mg-rich basaltic andesite in the western Bushveld Complex, South Africa . Contributions to Mineralogy and Petrology 171 , 57 . Google Scholar CrossRef Search ADS Mathez E. A. ( 1976 ). Sulfur solubility and magmatic sulfides in submarine basaltic glass . Journal of Geophysical Research 81 , 4269 – 4275 . Google Scholar CrossRef Search ADS McBirney A. R. ( 1996 ). The Skaergaard intrusion. In: Cawthorn, R. G. (ed.) Developments in Petrology , vol. 15: Layered Intrusions , pp. 147 – 180 . Naldrett A. J. , Wilson A. , Kinnaird J. , Yudovskaya M. , Chunnett G. ( 2012 ). The origin of chromitites and related PGE mineralization in the Bushveld Complex: new mineralogical and petrological constraints . Mineralium Deposita 47 , 209 – 232 . Google Scholar CrossRef Search ADS Nielsen R. L. ( 1990 ). Simulation of igneous differentiation processes. In: Nicholls J. , Russell J. K. (eds) Modern Methods of Igneous Petrology: Understanding Magmatic Processes. Reviews in Mineralogy 24 , 63 – 105 . Nielsen R. L. , Dungan M. A. ( 1983 ). Low-pressure mineral-melt equilibria in natural anhydrous mafic systems . Contributions to Mineralogy and Petrology 84 , 310 – 326 . Google Scholar CrossRef Search ADS Nikolaev G. S. , Ariskin A. A. , Barmina G. S. , Nazarov M. A. , Almeev R. R. ( 2016 ). Test of the Ballhaus–Berry–Green Ol–Opx–Sp oxybarometer and calibration of a new equation for estimating the redox state of melts saturated with olivine and spinel . Geochemistry International 54 , 301 – 320 . Google Scholar CrossRef Search ADS SACS (South African Committee for Stratigraphy) . ( 1980 ). Kent L. E. (Compiler) Stratigraphy of South Africa . Vol. 8 . Pretoria : Geological Survey of South Africa, Handbook , p. 690 . Shvarov Y. V. ( 1981 ). A general equilibrium criterion for an isobaric-isothermal model of a chemical system . Geochemistry International 18 , 38 – 45 . Wilson A. H. ( 2012 ). A chill sequence to the Bushveld Complex: insight into the first stage of emplacement and implications for the parental magmas . Journal of Petrology 53 , 1123 – 1168 . Google Scholar CrossRef Search ADS Wilson A. H. ( 2015 ). The earliest stages of emplacement of the Eastern Bushveld Complex: development of the lower zone, marginal zone and basal ultramafic sequence . Journal of Petrology 56 , 347 – 388 . Google Scholar CrossRef Search ADS Yudovskaya M. A. , Naldrett A. J. , Woolfe J. A. S. , Costin G. , Kinnaird J. A. ( 2015 ). Reverse compositional zoning in the Uitkomst chromitites as an indication of crystallization in a magmatic conduit . Journal of Petrology 56 , 2373 – 2394 . Google Scholar CrossRef Search ADS Yudovskaya M. A. , Kinnaird J. A. , Sobolev A. V. , Kuzmin D. V. , McDonald I. , Wilson A. H. ( 2013 ). Petrogenesis of the Lower Zone olivine-rich cumulates beneath the Platreef and their correlation with recognized occurrences in the Bushveld Complex . Economic Geology 108 , 1923 – 1952 . Google Scholar CrossRef Search ADS © The Author(s) 2018. Published by Oxford University Press. All rights reserved. For permissions, please e-mail: journals.permissions@oup.com This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/about_us/legal/notices) http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Journal of Petrology Oxford University Press

The COMAGMAT-5: Modeling the Effect of Fe–Ni Sulfide Immiscibility in Crystallizing Magmas and Cumulates

Loading next page...
 
/lp/ou_press/the-comagmat-5-modeling-the-effect-of-fe-ni-sulfide-immiscibility-in-lQqqXuxG7Z
Publisher
Oxford University Press
Copyright
© The Author(s) 2018. Published by Oxford University Press. All rights reserved. For permissions, please e-mail: journals.permissions@oup.com
ISSN
0022-3530
eISSN
1460-2415
D.O.I.
10.1093/petrology/egy026
Publisher site
See Article on Publisher Site

Abstract

ABSTRACT Details of the calibration and testing of the first ‘sulfide version’ of the COMAGMAT magma crystallization model (version 5·2, 2012–2014) are presented. The model’s updated empirical basis includes new mineral–melt geothermometers for olivine, plagioclase, high-Ca pyroxene, pigeonite, and orthopyroxene (calibrated at 1 atm. pressure), which are combined with a recently proposed Fe–Ni sulfide solubility model. This allows COMAGMAT-5 to be used for calculations of equilibrium and fractional crystallization of S-saturated and S-undersaturated magmas, including changes in the Fe/Ni ratio in silicate melts, femic minerals, and coexisting sulfides, as well as sulfide-silicate (±Fe–Ti oxides) proportions for multiply-saturated mineral assemblages. Based on our experience in tests of experimental data and modeling crystallization of mafic magmas, the possible range of application of COMAGMAT-5 may be extended up to 1–2 kbar pressure. The new model suggests a strong dependence of sulfide liquid immiscibility on the Ni content of the melt, as the increase of Ni is shown to decrease sulfide solubility, thus stabilizing the sulfide liquid in crystallizing mineral assemblages. This effect was tested for a Ni-rich sulfide-saturated glass dredged from the the southern Mid-Atlantic Ridge. The sulfide COMAGMAT was found to accurately predict the low S content observed in the glass (∼600 ppm); other existing ‘FeS’ solubility models yield higher values, mostly > 1000 ppm. In addition, using an experimentally studied high magnesia andesite composition, the results of calculations with two versions of COMAGMAT (3·72 & 5·2) were compared with those produced by MELTS family models. Application examples of COMAGMAT-5 include the modeling Fe–Ni sulfide saturation during equilibrium and fractional crystallization of ultramafic systems, approximating the most primitive magmas and cumulates from the Bushveld Complex in South Africa. INTRODUCTION The thermodynamic aspect of the genetic reconstructions of sulfide saturation history in mafic to ultramafic systems requires evaluating sulfide immiscibility timing and sulfide liquid composition at different stages of magmatic evolution. This involves constraints on the phase and chemical compositions of the crystallizing magma and, or, cumulates, particularly focusing on the equilibrium between the residual silicate melt and immiscible sulfide liquid. Some of these issues have been addressed in the ‘sulfide version’ of the COMAGMAT-5 magma crystallization model (Ariskin et al., 2009, 2013). This version combines an updated sulfide solubility model with newly-calibrated silicate mineral-melt geothermometers in a single algorithm designed to assess the effect of the Fe/Ni sulfide immiscibility on the compositional evolution of silicate melts, Fe–Mg minerals, and coexisting sulfides. Another feature of the new COMAGMAT is its ability to calculate sulfide-silicate (±Fe–Ti oxides) proportions for multiple-saturated mineral assemblages, and their effect on rock/mineral geochemistry during post-cumulus solidification of cumulate piles. Using these unique capabilities, the sulfide model has been used recently to quantify the geochemical evolution of an immiscible sulfide liquid which originated in olivine cumulates from the Yoko-Dovyren layered intrusion (Ariskin et al., 2016), as well as for quantifying temperature-compositional parameters of a parental magma calculated for the Chiney Complex in Southern Siberia, Russia (Gongalsky et al., 2016). In this study, we provide details on the calibration and other capabilities of the most recent COMAGMAT-5·2·2 model (http://geo.web.ru/∼ariskin/soft.html-id=comagmat.htm) and propose a methodology for its application to sulfide-saturated and sulfide-undersaturated systems. MAIN FEATURES OF THE EARLIER COMAGMAT MODELS COMAGMAT (Ariskin et al., 1993; Ariskin, 1999; Ariskin & Barmina, 2004) and MELTS (Ghiorso & Sack, 1995; Ghiorso et al., 2002; Gualda et al., 2012) are two families of popular thermodynamic programs, developed to model crystallization of mafic to ultramafic and silicic magmas. Both series of models are internally consistent and have some limitations when applied to volcanic systems and layered intrusions (Ghiorso, 1997; Ariskin & Barmina, 2004; Gualda & Ghiorso, 2015). Initially, COMAGMAT-3 was restricted to modeling low-pressure tholeiitic systems, such as mid-ocean ridge basalt (MORB) and large igneous province (LIP) magmas that follow the following order of crystallization: Ol → Pl → Aug → Pig (or Opx) → Fe–Ti oxides (Ariskin et al., 1993). Later it was updated to model crystallization of hydrous calc-alkaline magmas, starting from Ol + Aug cotectics at elevated pressures; see applications of COMAGMAT-3·5 to high-Mg magmas of the Klyuchevskoi volcano (Ariskin, 1999; Ariskin & Barmina, 2004) and to recent calculations for the parental magmas of the Bezymianny volcano in Kamchatka (Almeev et al., 2013). An important feature of COMAGMAT is that it allows for modeling crystallization of pigeonite and orthopyroxene, which is important for modeling the late crystallization stages of silica-oversaturated magmas. All previous versions of COMAGMAT (3·0–3·72, http://geo.web.ru/∼ariskin/soft.html-id=comagmat.htm) allowed for an alternative choice of ‘Pig’ or ‘Opx’ crystallization to model tholeiitic and calc-alkaline magmas. Thermodynamic fundamentals of COMAGMAT Incremental simulation of silicate melt equilibrium crystallization involves calculating mineral-melt equilibria in an evolving heterogeneous system, which consists of a residual melt (l) and crystallizing minerals (1 ≤ j ≤ M): nil=nibulk−∑j=1M∑r=1R(j)νi(r)jnrj, (1) where i is an oxide component (1 ≤ i ≤ N), r refers to end-members of a given j-mineral solid solution (1 ≤ r ≤ R(j)), nibulk,nil,nrj are the amounts of i-component in the system, melt, and crystallizing end-members of j-mineral, respectively, and νi(r)j represents the number of moles of i-component in one formula unit of an end-member. At a constant pressure, the mass-balance constraint (1) underlies calculations of the Gibbs free energy of a partially crystallized system at a set of independent (mainly P–T) parameters: G=∑i=1Nnilμil+∑j=1M∑r=1R(j)nrjμrj, (2) where μil and μrj are the chemical potentials of i-component in the melt and r-end-member of j-mineral. There are two main approaches to the search for the minimum Gibbs free energy at each step of crystallization (Frenkel & Ariskin, 1984). The first approach allows the isobaric potential at a given temperature to be minimized. This can be achieved through a number of algorithms based on nonlinear mathematical programming (e.g. Shvarov, 1981; Ghiorso, 1994; Karpov et al., 2001) utilizing an internally-consistent set of thermodynamic properties for the chosen components (e.g. Berman, 1988), thus ensuring that the model is applicable to a range of compositions and P–T conditions. An alternative approach is based on minimizing G at a given extent of crystallization (F), i.e. at fixed proportions of melt and crystals, searching for the equilibrium temperature (Frenkel & Ariskin, 1984; Ariskin et al., 1993). An advantage of this method is that it can account for the strong nonlinear dependence of T vs. F for multiply-saturated systems. Construction of crystallization trajectories as a function of the crystallization degree Fk=Fk−1+ΔF (where ΔF is an increment at each k-step of modeling) via small increments (e.g. ΔF ≤ 0·01) results in fast and accurate identification of inflection points which correspond to changes in crystallizing mineral assemblages. This capability is important for modeling the evolution of multi-mineral cotectics (e.g. Ol + Pl + Cpx → Ol + Pl + Cpx + Opx  → Ol + Pl + Cpx + Opx + Mt) characterized by a small dependence of their crystallization temperature on F. Another advantage of this approach is that it allows a finite-difference computational algorithm to be implemented when constructing magma differentiation models, which combine thermal and dynamical constraints of magma chambers with calculated crystallization sequences (Frenkel et al., 1989; Ariskin et al., 1993; Ariskin & Yaroshevsky, 2006). Basic principles of the COMAGMAT algorithm In a general form, the magma crystallization process can be described as a sequence of interdependent crystallization reactions producing end-members of rock-forming minerals EMrj (e.g. Fo, An, En) from the melt components Li (e.g. SiO2, MgO, NaAlO2 or Na2SiO3): ∑i=1Nνi(r)jLi=EMrj, (3) where νi(r)j are the stoichiometric coefficients of reaction (3), which correspond to the end-member components assumed in Eq. (1). Values of the equilibrium constants for the crystallization reactions (3) are given by the law of mass action: Krj=arj/∏i=1Naiνi(r)j, (4) where arj are the activities of end-members in the crystallizing minerals and ai are the activities of the melt components. In fact, the choice of thermodynamic components in the melt is arbitrary, as it depends on the accepted melt activity model. The lack of a comprehensive thermodynamic theory for multicomponent silicate melts leads to utilization of simplified descriptions of silicate melts, which involve a range of components unlikely to be physically present (e.g. NaAlO2; Nielsen & Dungan, 1983; or Fe2O3 and Na2SiO3; Ghiorso & Sack, 1995). The semi-empirical character of thermodynamic models for silicate melts is a feature of all current models of mineral–melt equilibria. The initial COMAGMAT model employed an ideal melt mixing model modified after Nielsen’s ‘two-lattice model’ (Ariskin et al., 1993). Assuming also ideal mixing for end-member components in mineral solutions, the link between the equilibrium constant (4) and the equilibrium temperature may be simplified to: ln Krj=ln xrj/∏i=1Nxiνi(r)j=−ΔG(3)/RgT=Arj/Trj+Brj, (5) where xrj are the concentrations of end-members in the crystallizing minerals, xi are the concentrations of the postulated melt components, and Arj>0 is the Arrhenius term representing the normalized enthalpy of r-phase melting, which is equivalent to the negative change in enthalpy of reaction (3), Arj=−ΔH(3)rj/Rg ( Rgis the gas constant). Modeling a k-step of equilibrium crystallization, Fk=Fk−1+ΔF initially involves the use of equation (1) to calculate a new melt composition ni(k)l by forming new amounts of the crystallizing mineral end-members nr(k)j (the total amount of new minerals is ΔF), using the same proportions as existed at the previous step Fk−1. This procedure can be described as a transition from the initial equilibrium state Fk−1 to a new non-equilibrated state due to a change in the melt composition following equation (1). Based on equations (1–5), a search for the minimum of the free energy function (2) for the new metastable system can be performed at the given crystallization degree Fk using a relationship described by Frenkel & Ariskin (1984): ∂G/∂nrj=−RgT ln Krj(o)+RgT ln Krj, (6) where Krj(o) is the reference state for pure endmembers r (i.e. forsterite or anorthite at their melting points) and Krj=Kr(k)j is calculated from (5) at each k-step of crystallization. Equation (6) is given to demonstrate that the decrease in the free energy of a metastable melt–mineral system at given conditions (P and Fk) can occur via iterative cycles of crystallization ( ∂nrj > 0) and/or dissolution ( ∂nrj < 0) of a number of end-members, leading to the equilibrium values of Kr(k)j following equation (4). Equation (6) forms the basis of the original COMAGMAT algorithm for free-energy minimization (Ariskin et al., 1993). The role of pseudo-liquidus mineral–melt temperatures Substituting K from Eq. (5), equation (6) can be re-written as ∂G/∂nrj=RgArj(T/Tj−1), (7) where Arj>0 represent parameters of corresponding mineral–melt geothermometers, T is the target (initially unknown) equilibrium temperature at step k, and Tj is the temperature of the current, notionally metastable state of each j-mineral ( Trj=Tj) with the current melt composition. Eq. (7) is used iteratively to determine the order and proportions of crystallization, including peritectic reactions at each step. For the completed derivation of expressions (6) and (7), see equations (25–45) in Ariskin & Barmina (2004). The COMAGMAT algorithm includes a special subroutine designed for calculating and comparing Tj (the so-called ‘pseudo-liquidus temperatures’, see Ariskin et al. (1993) and Danyushevsky & Plechov (2011)). Metastable mineral–melt temperatures Tj during iterations are calculated for each end-member component of every mineral present on the liquidus (i.e. nrj  > 0) following stoichiometry constraints: ∑r=1R(j)∏i=1Nxiνi(r)j exp ⁡(Arj/Tj+Brj)=1, (8) which require that the concentrations of end-members in each mineral sum to 1 ( ∑r=1R(j)xrj=1), where xrj=∏i=1Nxiνi(r)j exp ⁡(Arj/Tj+Brj) (9) The target temperature T must be within the range between the maximum and minimum of all pseudo-liquidus mineral temperatures Tj for the minerals currently present on the liquidus ( nj> 0 →nrj> 0). In this case, minimization in G is equivalent to minimization of the difference Tjmax⁡−Tj(nj>0)min⁡<εT, (10) where εT is the chosen computational accuracy during modelling. Detailed flowcharts of the main COMAGMAT algorithm describing its thermodynamic background are given by Ariskin & Barmina (2004). DEVELOPMENT OF COMAGMAT-5 COMAGMAT-5 (v. 5·2) is a new generation of the COMAGMAT programs designed to model crystallization of five silicate minerals (Ol–Pl–Aug–Pig–Opx), two Fe–Ti oxides (ilmenite and Fe–Ti spinel), and immiscible Fe–Ni sulfide liquids. This code was essentially rewritten, using more sophisticated compiling programs, as compared to previous versions of COMAGMAT (Ariskin et al., 1993; Ariskin, 1999); however, it preserved the main features of the basic algorithm described above. It includes a set of updated 1-atm mineral–melt geothermometers for olivine, plagioclase, augite, pigeonite, and orthopyroxene. Several new subroutines have been added, including SULSAT, which estimates the solubility of Fe–Ni sulfides based on the model of Ariskin et al. (2013). COMAGMAT-5 is also capable of modeling simultaneous crystallization of three pyroxenes including augite, pigeonite, and orthopyroxene. New silicate mineral–melt geothermometers The term ‘geothermometer’ is applied to a system of thermodynamic equations (5) which describe the partitioning of major and some minor elements between rock-forming minerals and silicate melts over a range of P–T–fO2–H2O conditions (Ariskin et al., 1993). The geothermometers are calibrated on a large dataset of melting and crystallization experiments in natural and synthetic systems ranging from ultramafic and high-Mg basalts to ferrobasalts and dacites (i.e. the INFOREX experimental database). Now the INFOREX melting-experiment database provides information on 15,418 sub-liquidus runs from 442 experimental studies published in the literature, see the INFOREX Reference List in Supplementary Data Electronic Appendix 1. The available information includes pressure, temperature, fO2, run duration, containers, 17415 compositions of 40 minerals and 8404 melt (glass) compositions, as well as signatures of nominally anhydrous and hydrous experiments in a wide range of natural and synthetic systems (Ariskin & Barmina, 2004; Nikolaev et al., 2016). Ariskin & Barmina (2004) have demonstrated that this approach results in geothermometers with a general precision of 10–15°C for each silicate mineral; these geothermometers can reproduce liquidus mineral compositions with a precision of 1–3 mol% for olivine and pyroxenes, and 2–4 mol% for plagioclase. Two major constraints underlie this approach. First, all mineral–melt geothermometers must be calibrated with respect to the equilibrium constants (5): ln Krj=a/T+b1 ln cl1+b2 ln cl2+⋯+d, (11) where Krj characterizes the reaction of r-end-member formation from melt components (3); cl1,cl2,… are empirically selected melt structural-chemical parameters; and a, b1,b2,… and d are the regression coefficients. Second, the Krj values need to be calculated using the same melt component activity model for each end-member and mineral. The above constraints are derived from the fundamental thermodynamic relationships (Ghiorso, 1987) to ensure that the crystallization model represents an internally-consistent formulation. Ignoring these constraints by selecting an arbitrary system of mineral–melt partitioning equations (Danyushevsky & Plechov, 2011) may result in incorrect proportions of crystallized minerals and large errors in the calculated liquid lines of descent (Ghiorso, 1987; Ariskin & Barmina, 2004). The Krj values in Eq. (11) are calculated using a two-lattice melt components model (Ariskin & Barmina, 1990) modified after Nielsen & Dungan (1983) and Nielsen (1990). Use of melt structural-chemical parameters cl (e.g. Si/O, Al/O, Al/Si) significantly improves the fit to the experimental data (Ariskin, 1999). The calibration data set To update geothermometers for olivine, plagioclase, augite, and pigeonite, we used an expanded set of coexisting mineral–melt compositions. The set includes ∼270 runs conducted at 1 atm (thus nominally anhydrous) for > 48 hours at temperatures from 1050–1250°C over a wide range of oxygen fugacity (IW - NNO + 1). These experimental glasses contain 7 ≤ FeO ≤ 18 wt %, 45 ≤ SiO2 ≤ 60 wt %, and 2 ≤ Na2O + K2O ≤ 5 wt %. This experimental array includes ∼150 mineral–melt pairs for olivine, 187 pairs for plagioclase, 125 pairs for high-Ca Cpx, and 43 pairs for pigeonite. To calibrate the new orthopyroxene–melt geothermometer an extra set of 96 experiments was added, which covers higher temperatures (to 1300°C) and more Si-rich compositions (48 ≤ SiO2 ≤ 68 wt %, see Supplementary Data Electronic Appendix 1; supplementary data are available for downloading at http://www.petrology.oxfordjournals.org). The calibration dataset was limited to 1065°C ≤T ≤1300°C as it is rare for volcanic suite and layered intrusion parental magmas to have higher temperatures (e.g. Ariskin & Barmina, 2004). The lower temperature range is included in the calibration dataset because protracted magma fractionation observed in many differentiated intrusions leads to highly evolved compositions (e.g. in the Skaergaard intrusion, the Layered Series rocks are composed of mineral assemblages corresponding to temperatures ∼1115–1080°C; McBirney, 1996; Ariskin, 2003). Since COMAGMAT-3 was calibrated at higher temperatures (∼1125–1350°C), calculations at and below 1100°C often resulted in over-expansion of the Ca-clinopyroxene field at the expense of plagioclase and pigeonite. Moreover, plagioclase compositions calculated for the most evolved melts were too sodic and potassic. Methods of calibration and results Calibrations of olivine and plagioclase geothermometers were carried out in the form of Eq. (11) by a least-squares method, using a two-stage procedure of multiple regression analysis (Ariskin, 1999; Ariskin & Barmina, 2004). For high-Ca clinopyroxene, pigeonite, and orthopyroxene a modified method was used, which is based on multiple regression analysis for normalized values: ln (Krj)′=ln Krj−aPxj/T=b1 ln cl1+b2 ln cl2+⋯+d, (12) where the Arrhenius terms aPxj for enstatite (En), ferrosilite (Fs), wollastonite (Wo), and the Al-component are postulated to be the same for all pyroxenes. Detailed comments about both approaches are given in Supplementary Data Electronic Appendix 2 (Part 1). Note that in addition to Mg, Fe, Ca, and Mn end-members, geothermometers for olivine and pyroxenes include the Ni end-members NiSi0·5O2 and NiSiO3, leading to a total system of 22 end-member geothermometers (Table 2·1 and Part 2 in Supplementary Data Electronic Appendix 2). The equations for each end-member were then combined into five mineral–melt equilibria models for olivine, plagioclase, augite, pigeonite, and orthopyroxene (Ariskin & Barmina, 2004). The proposed mineral–melt models have uncertainties for T of 10–15°C and for mineral compositions of ∼0·7–2 mol% of the end-member concentrations in olivine, orthopyroxene and augite, and of 3–4 mol% for pigeonite and plagioclase. Table 1: Results of sulfide solubility calculations for the BTJ glass at P = 1 atm SCSS model Melt parameters QFM QFM + 0·5 QFM + 1 Ni, ppm H2O, ppm T, °C SCSS, ppm T, °C SCSS, ppm T, °C SCSS COMAGMAT-5 (Ariskin et al., 2013) 310 0 1258·2 514* 1256·9 566* 1255·3 710* 0 0 1256·4 1241* 1255·0 1368* 1253·3 1710* Li & Ripley, 2009 0 515 1257 1054 1257 1042 1257 1028 Liu et al., 2007 0 515 1257 1027 1257 1002 1257 973 Fortin et al., 2015 (Model B) 0 515 1257 1121 – – – – SCSS model Melt parameters QFM QFM + 0·5 QFM + 1 Ni, ppm H2O, ppm T, °C SCSS, ppm T, °C SCSS, ppm T, °C SCSS COMAGMAT-5 (Ariskin et al., 2013) 310 0 1258·2 514* 1256·9 566* 1255·3 710* 0 0 1256·4 1241* 1255·0 1368* 1253·3 1710* Li & Ripley, 2009 0 515 1257 1054 1257 1042 1257 1028 Liu et al., 2007 0 515 1257 1027 1257 1002 1257 973 Fortin et al., 2015 (Model B) 0 515 1257 1121 – – – – * COMAGMAT-5 calculates the total sulfur content as SCSS (Sulfur Content at Sulfide Saturation), by adding S6+ using the model of Jugo (2009). Other models calculate only S2- as SCSS. Note that Model B by Fortin et al. (2015) does not take into account the effect of fO2. The major element composition of BTJ is given in Kamenetsky et al. (2013). Table 1: Results of sulfide solubility calculations for the BTJ glass at P = 1 atm SCSS model Melt parameters QFM QFM + 0·5 QFM + 1 Ni, ppm H2O, ppm T, °C SCSS, ppm T, °C SCSS, ppm T, °C SCSS COMAGMAT-5 (Ariskin et al., 2013) 310 0 1258·2 514* 1256·9 566* 1255·3 710* 0 0 1256·4 1241* 1255·0 1368* 1253·3 1710* Li & Ripley, 2009 0 515 1257 1054 1257 1042 1257 1028 Liu et al., 2007 0 515 1257 1027 1257 1002 1257 973 Fortin et al., 2015 (Model B) 0 515 1257 1121 – – – – SCSS model Melt parameters QFM QFM + 0·5 QFM + 1 Ni, ppm H2O, ppm T, °C SCSS, ppm T, °C SCSS, ppm T, °C SCSS COMAGMAT-5 (Ariskin et al., 2013) 310 0 1258·2 514* 1256·9 566* 1255·3 710* 0 0 1256·4 1241* 1255·0 1368* 1253·3 1710* Li & Ripley, 2009 0 515 1257 1054 1257 1042 1257 1028 Liu et al., 2007 0 515 1257 1027 1257 1002 1257 973 Fortin et al., 2015 (Model B) 0 515 1257 1121 – – – – * COMAGMAT-5 calculates the total sulfur content as SCSS (Sulfur Content at Sulfide Saturation), by adding S6+ using the model of Jugo (2009). Other models calculate only S2- as SCSS. Note that Model B by Fortin et al. (2015) does not take into account the effect of fO2. The major element composition of BTJ is given in Kamenetsky et al. (2013). Table 2: Compositions of magmatic melts and cumulate rocks used for testing COMAGMAT and MELTS on experimental data and example calculations with COMAGMAT-5 Melt components, wt % HMA-85-41c Proposed Bushveld magmas Cumulate rocks B1* UM** UM + 56 wt % Ol-91·6 Ol pyr*** SiO2 57·79 56·09 53·38 46·29 46·03 TiO2 0·60 0·28 0·39 0·17 0·10 Al2O3 14·46 11·31 9·73 4·28 4·86 FeO* 5·74 9·17 10·19 9·18 9·21 MnO 0·11 0·17 0·17 0·16 0·16 MgO 9·14 13·58 19·00 36·39 35·13 CaO 8·17 6·34 5·36 2·36 3·04 Na2O 3·11 1·43 1·34 0·59 0·35 K2O 0·71 1·05 0·54 0·24 0·14 P2O5 0·15 0·07 0·08 0·04 0·02 NiO – 0·04 0·09 0·32 0·29 Cr2O3 – 0·25 – – 0·90 Total 99·98 99·78 100·27 100 100·23 S, ppm – 438**** 368**** 162**** – Melt components, wt % HMA-85-41c Proposed Bushveld magmas Cumulate rocks B1* UM** UM + 56 wt % Ol-91·6 Ol pyr*** SiO2 57·79 56·09 53·38 46·29 46·03 TiO2 0·60 0·28 0·39 0·17 0·10 Al2O3 14·46 11·31 9·73 4·28 4·86 FeO* 5·74 9·17 10·19 9·18 9·21 MnO 0·11 0·17 0·17 0·16 0·16 MgO 9·14 13·58 19·00 36·39 35·13 CaO 8·17 6·34 5·36 2·36 3·04 Na2O 3·11 1·43 1·34 0·59 0·35 K2O 0·71 1·05 0·54 0·24 0·14 P2O5 0·15 0·07 0·08 0·04 0·02 NiO – 0·04 0·09 0·32 0·29 Cr2O3 – 0·25 – – 0·90 Total 99·98 99·78 100·27 100 100·23 S, ppm – 438**** 368**** 162**** – * HMA-85–41c is a high-Mg andesite used in melting experiments by Grove et al. (2003). Proposed magmas and cumulate rocks represent mafic to ultramafic (UM) compositions assumed to be parental for the Lower Zone and the Basal Ultramafic Sequence of the Bushveld Complex (Wilson, 2012, 2015). See column 4 in table 3 of Wilson (2012), **UM liquid composition from Schwerin fold chill listed as CH12/7&14UM in column 3 of table 4 in Wilson (2015), ***Olivine pyroxenite listed as CH7/91 in table 3 of Wilson (2015). ****S concentrations used in the COMAGMAT-5 modeling: 438 ppm S in B1 is taken from Barnes et al. (2010), 368 ppm initial S concentration estimated based on allowing UM to crystallize 16% olivine, which produces S concentration of 438 ppm, similar to that seen in B1; S in ‘UM + 56% Ol-91·6’ corresponds to a mixture of 44 wt.% of the UM composition with 56 wt.% of olivine mg#91.6. Table 2: Compositions of magmatic melts and cumulate rocks used for testing COMAGMAT and MELTS on experimental data and example calculations with COMAGMAT-5 Melt components, wt % HMA-85-41c Proposed Bushveld magmas Cumulate rocks B1* UM** UM + 56 wt % Ol-91·6 Ol pyr*** SiO2 57·79 56·09 53·38 46·29 46·03 TiO2 0·60 0·28 0·39 0·17 0·10 Al2O3 14·46 11·31 9·73 4·28 4·86 FeO* 5·74 9·17 10·19 9·18 9·21 MnO 0·11 0·17 0·17 0·16 0·16 MgO 9·14 13·58 19·00 36·39 35·13 CaO 8·17 6·34 5·36 2·36 3·04 Na2O 3·11 1·43 1·34 0·59 0·35 K2O 0·71 1·05 0·54 0·24 0·14 P2O5 0·15 0·07 0·08 0·04 0·02 NiO – 0·04 0·09 0·32 0·29 Cr2O3 – 0·25 – – 0·90 Total 99·98 99·78 100·27 100 100·23 S, ppm – 438**** 368**** 162**** – Melt components, wt % HMA-85-41c Proposed Bushveld magmas Cumulate rocks B1* UM** UM + 56 wt % Ol-91·6 Ol pyr*** SiO2 57·79 56·09 53·38 46·29 46·03 TiO2 0·60 0·28 0·39 0·17 0·10 Al2O3 14·46 11·31 9·73 4·28 4·86 FeO* 5·74 9·17 10·19 9·18 9·21 MnO 0·11 0·17 0·17 0·16 0·16 MgO 9·14 13·58 19·00 36·39 35·13 CaO 8·17 6·34 5·36 2·36 3·04 Na2O 3·11 1·43 1·34 0·59 0·35 K2O 0·71 1·05 0·54 0·24 0·14 P2O5 0·15 0·07 0·08 0·04 0·02 NiO – 0·04 0·09 0·32 0·29 Cr2O3 – 0·25 – – 0·90 Total 99·98 99·78 100·27 100 100·23 S, ppm – 438**** 368**** 162**** – * HMA-85–41c is a high-Mg andesite used in melting experiments by Grove et al. (2003). Proposed magmas and cumulate rocks represent mafic to ultramafic (UM) compositions assumed to be parental for the Lower Zone and the Basal Ultramafic Sequence of the Bushveld Complex (Wilson, 2012, 2015). See column 4 in table 3 of Wilson (2012), **UM liquid composition from Schwerin fold chill listed as CH12/7&14UM in column 3 of table 4 in Wilson (2015), ***Olivine pyroxenite listed as CH7/91 in table 3 of Wilson (2015). ****S concentrations used in the COMAGMAT-5 modeling: 438 ppm S in B1 is taken from Barnes et al. (2010), 368 ppm initial S concentration estimated based on allowing UM to crystallize 16% olivine, which produces S concentration of 438 ppm, similar to that seen in B1; S in ‘UM + 56% Ol-91·6’ corresponds to a mixture of 44 wt.% of the UM composition with 56 wt.% of olivine mg#91.6. Figure 1 shows the fit of the new models to the experimental data used for calibration. The new olivine model was also tested at higher temperatures on a set of 66 experiments at 1220°C  < T <1350°C (9–16 wt % MgO in the melt, 83–89% Fo), which were not included in the calibration dataset (Fig. 1a), demonstrating that the new model for olivine is applicable to at least 1350°C. The complete set of calculated temperatures and compositions and experimental data is given in Supplementary Data Electronic Appendix 1. The new calibrations for plagioclase and clinopyroxene resulted in a better fit of the calculated mineral compositions to the experimental data than achieved by the previous version of COMAGMAT-3.5 (Fig. 1c, d). Fig. 1. View largeDownload slide Comparison of the accuracy of the COMAGMAT-5 and COMAGMAT-3.5 calculations for equilibrium temperatures and mineral compositions. The dashed lines represent equal values. The temperature–compositional calculations by COMAGMAT were conducted on the same experimental arrays as those extracted from the INFOREX database (Ariskin & Barmina, 2004) and used in the calibration of mineral–melt geothermometers for the COMAGMAT-5 model, see details in the text (‘The calibration data set’ section). To test the validity of the model at higher temperatures, an additional set of 66 experimental Ol–melt pairs were used (1220°C < T <1350°C, 9–16% MgO in the melt, 83–89% Fo, see crosses in Fig. 1a). Note that the previous version of COMAGMAT-3.5 (Ariskin, 1999) worked well for low-alkaline tholeiitic systems; however, it essentially over-estimated the composition of plagioclase modeled in evolved crystallization products of mildly-alkaline magmas. Fig. 1. View largeDownload slide Comparison of the accuracy of the COMAGMAT-5 and COMAGMAT-3.5 calculations for equilibrium temperatures and mineral compositions. The dashed lines represent equal values. The temperature–compositional calculations by COMAGMAT were conducted on the same experimental arrays as those extracted from the INFOREX database (Ariskin & Barmina, 2004) and used in the calibration of mineral–melt geothermometers for the COMAGMAT-5 model, see details in the text (‘The calibration data set’ section). To test the validity of the model at higher temperatures, an additional set of 66 experimental Ol–melt pairs were used (1220°C < T <1350°C, 9–16% MgO in the melt, 83–89% Fo, see crosses in Fig. 1a). Note that the previous version of COMAGMAT-3.5 (Ariskin, 1999) worked well for low-alkaline tholeiitic systems; however, it essentially over-estimated the composition of plagioclase modeled in evolved crystallization products of mildly-alkaline magmas. Incorporation of Fe–Ni sulfide solubility model Recently, an additional subroutine (SULSAT) has been developed for calculating Sulfur Contents at Sulfide Saturation (SCSS, Campbell & Naldrett, 1979) in mafic magmas, as a function of pressure, temperature, fO2, major components, and Ni content in the melt (Ariskin et al., 2013). In contrast to other SCSS models, which approximate the immiscible sulfides by a pure FeS pyrrhotite liquid, see comments and references in Supplementary Data Electronic Appendix 2 (Part 3), we argued for a strong effect of minor components (such as Ni and Cu) which may be present in the natural sulfide liquids in amounts up to 10–20 wt %. Our SCSS model proposes the existence of positively charged Fe–Ni sulfide complexes in the silicate melt, which are formed as a result of complexation reactions between the sulfide-forming ions (Fe2+, Ni2+, S2-) and (Fe, Ni)S species. The proposed mechanism of sulfide solubility was calibrated and tested against 290 nominally anhydrous experimental and natural glass compositions (both Ni-free and Ni-bearing), producing a good fit between the observed and calculated variations of S content vs. temperature and Fe content in the melt (Ariskin et al., 2013). Integration of the SULSAT model into COMAGMAT-5 Incorporation of the SULSAT subroutine into the COMAGMAT-5 model required the addition of Ni and S as the principal thermodynamic components. This resulted in several changes to the algorithm, which allowed for combining the Fe–Ni sulfide model with equilibrium temperature and composition modeling for silicate minerals and Fe–Ti oxides. For initially S-undersaturated compositions, the algorithm uses SCSS calculations to determine the onset of sulfide immiscibility. This is achieved by comparing the modeled S content in the melt (which behaves as an incompatible component until sulfide saturation is reached) with the calculated SCSS at each stage of crystallization (e.g. Barnes & Lightfoot, 2005; Li & Ripley, 2005; Barnes, 2007; Barnes et al., 2010). When the calculated SCSS is the same as the S content in the melt, it is taken to indicate that the system has become sulfide saturated (see figs 9 and 10 in Ariskin et al. (2013)). Starting from this point, mass-balance constraints are used to calculate both the proportions and the compositions of all phases, including Fe–Ni sulfide. The algorithm also takes into account possible S over-saturation at the onset of crystallization, thus allowing for modeling systems containing an excess of sulfide phases. To account for the effect of fO2 on S speciation in the melt, we followed the approach of Jugo (2009), which allows for calculating the proportion of S6+ in a silicate melt. The value of bulk SCSS becomes a function of fO2 and increases with fO2 for a given major element composition. This is due to the addition of some sulfate component (as S6+) to the sulfide species ((Fe, Ni)S-complexes) at more oxidized conditions. The main options and capabilities of COMAGMAT-5 There are three main options in COMAGMAT-5: (i) calculating sulfur solubility for a given temperature and melt composition; (ii) calculating equilibrium temperatures for silicate and oxide minerals and their compositions for a given melt composition; and (iii) modeling crystallization processes for a given magma composition. A detailed overview of the main options and the updated structure of COMAGMAT-5 is given as Supplementary Data Electronic Appendix 2 (Part 4); examples of output files are available as spreadsheets in Supplementary Data Electronic Appendix 1. The first option allows SCSS to be calculated for a number of melt compositions, if their temperatures and fO2 are known. The SCSS calculations can be performed at a range of pressures up to 1 GPa. The second option allows both the temperature of crystallization and the composition of the mineral on the liquidus at 1 atm pressure to be calculated for a given melt composition. This routine is useful for testing the accuracy of the thermometric and compositional calculations in COMAGMAT-5 using experimental data. The third option allows crystallization to be modeled for multiply-saturated magmatic systems, including five silicate minerals (olivine, plagioclase, augite, pigeonite, orthopyroxene), two Fe–Ti oxides (ilmenite, magnetite), and a Fe–Ni sulfide liquid. Calculations can be performed for the cases of: (1) equilibrium crystallization in a closed system; (2) perfect fractional crystallization (assuming complete separation of crystallizing solids from the melt); and (3) an intermediate case where a certain portion of crystallizing minerals are postulated to re-equilibrate with the residual melt, whereas the remainder of the minerals are removed from the system (Nielsen, 1990; Ariskin & Barmina, 2004; Danyushevsky & Plechov, 2011). Changes in the user-interface and reorganization of the input/output files made COMAGMAT-5 more flexible and easy to use (see examples of its interface in Supplementary Data Electronic Appendix 2, Part 4). Currently, COMAGMAT-5 is the only algorithm capable of modeling co-crystallization of Aug±Pig±Opx, including calculations of peritectic relationships between orthopyroxene, pigeonite, and olivine, as well as between magnetite and ilmenite in plagioclase-pyroxene saturated systems. As compared to previous versions of COMAGMAT-3.n (Ariskin et al., 1993; Ariskin, 1999), the updated model for orthopyroxene provides more accurate calculations of its melting/crystallization proportions relative to olivine. Two examples of testing the COMAGMAT-5 on natural and experimental data are given in the section below. VERIFICATION OF THE SULFIDE COMAGMAT Test using a high-Ni natural glass composition The SCSS model used in COMAGMAT-5 (Ariskin et al., 2013) predicts that Ni has a pronounced effect on decreasing SCSS, despite its relatively low concentrations in silicate melts. Here we present a demonstration of the effect of Ni on SCSS using a sulfide-saturated glass composition dredged from the southern Mid-Atlantic Ridge near the Bouvet Triple Junction, sample S18–60/1 (the BTJ glass, Kamenetsky et al., 2001, 2013). This is a high-Mg (8·2 wt % MgO, mg# 67·1) and high-Ni (310 ppm) andesite (57·3 wt % SiO2) which is depleted in Cu (33 ppm) and S (on average 600 ppm), compared to typical MORB glasses (Mathez, 1976; Ariskin et al., 2013). A distinctive feature of this sample is the presence of numerous spherical Ni-rich sulfide globules, which have been interpreted to indicate sulfide-saturation of the melt (Kamenetsky et al., 2013). A detailed examination of the largest globule (∼180 μm in diameter), combining X-ray mapping and microprobe data, allowed the authors to estimate the average Ni and Fe contents in the sulfide melt; this melt is characterized by a high Ni/(Ni + Fe) molar ratio of ∼0·40. Rare micro-phenocrysts of olivine (Fo86·9) suggest the BTJ melt has olivine on its liquidus, thus allowing its liquidus temperature to be estimated by using olivine-melt geothermometers. Calculations of SCSS for the BTJ glass Falloon et al. (2007) have demonstrated that the olivine-melt geothermometer of Ford et al. (1983) accurately estimates crystallization temperatures of MORB melts at low pressure. Calculations for the BTJ glass, using the model of Ford et al. (1983) and COMAGMAT-5 at P = 1 atm and QFM buffer, resulted in 1255·2°C (Fo86·3) and 1258·2°C (Fo85·9), respectively. The ‘Sulfur solubility’ option of COMAGMAT-5 (see Supplementary Data Electronic Appendix 2) was used to calculate the SCSS for the BTJ glass composition to compare the calculated values with the measured S content in the glass (Table 1). The calculations were performed at the calculated liquidus temperature at 1 atm pressure and fO2 ranging from QFM to QFM + 1. Relatively oxidized conditions were selected to account for the oxygen contents observed in the sulfide globules from BTJ: 0·33–0·55 wt % O (Kamenetsky et al., 2013). Assuming 310 ppm Ni in the melt, COMAGMAT-5 calculated bulk S solubility ranging from 514 ppm S at QFM (4·5 % S presented as sulfate species S6+), to 566 ppm S at QFM + 0·5 (12·5% S6+), and 710 ppm S at QFM + 1 (29·3% S6+), see Table 1. The value at QFM + 0·5 is close to the 600 ppm S measured in the BTJ glass. Similar calculations based on the Li & Ripley (2009) equation resulted in SCSS values around 1040 ppm. We also tested the Liu et al. (2007) equation, which postulates a strong ‘dropping effect’ of H2O on SCSS at H2O < 0·1 wt%, see figure 3 in Ariskin et al. (2013). Using the H2O content determined in the BTJ glass (515 ppm H2O), these calculations resulted in the SCSS ranging from 1027 to 973 ppm (Table 1). The use of Model-B from Fortin et al. (2015) at the same parameters resulted in ∼1121 ppm S. As a result, none of the SCSS models, except COMAGMAT-5, was found to replicate the low sulfide solubility observed in the high-Ni BTJ glass. If we assume that the BTJ glass is Ni-free, COMAGMAT-5 calculates an SCSS of >1240 ppm (Table 1), which is similar to other models that do not take Ni content in the silicate melt into account. Thus, the 1·5–2 times overestimation of the SCSS based on pure FeS models is clearly due to the high Ni content observed in the glass. Results of COMAGMAT-5 modeling at QFM + 0·5 indicate that this melt should be in equilibrium with a high Ni sulfide (Ni/(Ni + Fe)=0·370), which is consistent with the Ni/(Ni + Fe)∼0·408 estimated for the BTJ sulfide globule (Kamenetsky et al., 2013). Test using experimental data The second test was carried out using the ‘Modeling crystallization’ option of COMAGMAT-5 (Supplementary Data Electronic Appendix 2, Part 4). To test the validity of the phase equilibria calculations, we selected results of eight 1 atm. synthesis melting experiments carried out on the 85–41c primitive magnesian andesite from Mt. Shasta (Grove et al., 2003), see Table 2. Two simulations of the equilibrium crystallization for this melt were carried out. The first was done using the basic version of the COMAGMAT-5, whereas the second modeling was performed using a slight correction of the calculated temperatures (see captions to Fig. 2a). In addition, we used a previous version of COMAGMAT-3.n to calculate the crystallization trajectory for the same melt. All simulations were carried out at 1 atm. and QFM, similar to conditions of the experimental study (Grove et al., 2003). Fig. 2. View largeDownload slide Experimental and modeled sequences for equilibrium crystallization of the high-Mg andesite 85–41c (Grove et al., 2003) at P = 1 atm and QFM. The dashed lines represent experimental liquidus temperatures for olivine, orthopyroxene, and clinopyroxene & plagioclase. (A) Calculations using the COMAGMAT-5·2 model were conducted without corrections for the modeled mineral temperatures and using small temperature shifts with respect to the original calculated values (−10°C for olivine, -12°C for plagioclase, -3°C for high-Ca pyroxene, and +5°C for orthopyroxene), see option 4·5 in Supplementary Data Electronic Appendix 2 (Part 4). For comparison, results of simulations by COMAGMAT-3·72 are shown. Both models were downloaded from http://geo.web.ru/∼ariskin/soft.html-id=comagmat.htm. All crystallization trajectories were constructed with crystallization increments from 0·5% up to 80% crystallized. The composition of the HMA is listed in Table 2. (B) MELTS calculations (Ghiorso & Sack, 1995) were carried out using the Adiabat_1ph version 1·8 run by an EXCEL-based front-end using a shell command (developed by J.-I. Kimura, JAMSTEC, Japan), whereas pMELTS (Ghiorso et al., 2002) and RMELTS (Gualda et al., 2012) codes were downloaded from http://melts.ofm-research.org/unix.html. The MELTS crystallization trajectory was modeled with a 10°C temperature increment in the 1300–1000°C range; calculations using pMELTS and RMELTS were carried out in the 1265–1000°C range with a 5°C temperature increment. The calculations using pMELTS and RMELTS resulted in a wide field of spinel (with onset of crystallization temperatures of 1160°C and 1220°C, correspondingly), not shown in the Figure 2b plots. A complete series of EXCEL-files including results of the modeling is available upon request from the authors. Fig. 2. View largeDownload slide Experimental and modeled sequences for equilibrium crystallization of the high-Mg andesite 85–41c (Grove et al., 2003) at P = 1 atm and QFM. The dashed lines represent experimental liquidus temperatures for olivine, orthopyroxene, and clinopyroxene & plagioclase. (A) Calculations using the COMAGMAT-5·2 model were conducted without corrections for the modeled mineral temperatures and using small temperature shifts with respect to the original calculated values (−10°C for olivine, -12°C for plagioclase, -3°C for high-Ca pyroxene, and +5°C for orthopyroxene), see option 4·5 in Supplementary Data Electronic Appendix 2 (Part 4). For comparison, results of simulations by COMAGMAT-3·72 are shown. Both models were downloaded from http://geo.web.ru/∼ariskin/soft.html-id=comagmat.htm. All crystallization trajectories were constructed with crystallization increments from 0·5% up to 80% crystallized. The composition of the HMA is listed in Table 2. (B) MELTS calculations (Ghiorso & Sack, 1995) were carried out using the Adiabat_1ph version 1·8 run by an EXCEL-based front-end using a shell command (developed by J.-I. Kimura, JAMSTEC, Japan), whereas pMELTS (Ghiorso et al., 2002) and RMELTS (Gualda et al., 2012) codes were downloaded from http://melts.ofm-research.org/unix.html. The MELTS crystallization trajectory was modeled with a 10°C temperature increment in the 1300–1000°C range; calculations using pMELTS and RMELTS were carried out in the 1265–1000°C range with a 5°C temperature increment. The calculations using pMELTS and RMELTS resulted in a wide field of spinel (with onset of crystallization temperatures of 1160°C and 1220°C, correspondingly), not shown in the Figure 2b plots. A complete series of EXCEL-files including results of the modeling is available upon request from the authors. Comparisons of the COMAGMAT calculation results with those obtained in experiments are shown in Figure 2a. These data demonstrate that the basic version of COMAGMAT-5·2 is able to predict the major minerals crystallizing in the experimental range; however, it overestimates by ∼20°C the liquidus temperature for olivine and underestimates by ∼20°C that for orthopyroxene. The modeled temperatures for plagioclase and augite are accurate, if the range of 15°C between two experimental points where these phases are absent and present (Grove et al., 2003) is taken into account. At temperatures below 1120°C, COMAGMAT-5 predicts a peritectic reaction between orthopyroxene and pigeonite and the occurrence of ilmenite. COMAGMAT-3·72 demonstrates a similar sequence of crystallization; however, the modeled olivine liquidus temperature is slightly higher as compared to COMAGMAT-5·2, whereas the modeled onset of orthopyroxene crystallization is delayed by 15°C (Fig. 2a). Examples of calculations with the corrected mineral temperatures are given in Figure 2a (see plot labeled ‘COMAGMAT 5·2 corrected’) to demonstrate the COMAGMAT-5 capability for ‘fine tuning’ the model, in order to make the modeling more consistent with available experimental data or natural observations. It is possible to see that even small shifts of the modeled temperatures around ±10°C or less allow construction of the crystallization order exactly as observed in the experiments (Grove et al., 2003): Ol → Ol + Opx → Opx only → Opx + Pl + Aug. Note that it would be difficult (or perhaps impossible) to do this using crystallization programs, which calculate phase equilibria at predetermined temperatures. Figure 2b represents similar comparisons of the thermodynamic calculations on the same HMA (high magnesium andesite) composition, using three models of the MELTS family, including the original MELTS (Ghiorso & Sack, 1995), pMELTS (Ghiorso et al., 2002), and RMELTS (Gualda et al., 2012). Similar to COMAGMAT-3·72 and COMAGMAT-5·2, both MELTS and RMELTS are accurate enough to predict the experimental crystallization sequence, whereas pMELTS strongly underestimates crystallization temperatures for clinopyroxene and, particularly, orthopyroxene. Note that the pMELTS calculations were carried out at pressures far from the recommended 1–4 GPa range of pressures (Ghiorso et al., 2002). Generally, both series of plots in Figure 2a and b demonstrate that modern magma crystallization models should have a user-friendly capability to slightly change or to correct input thermodynamic parameters, allowing petrologists to make the results of their calculations more consistent with results obtained using experimental data. Data for the corrected COMAGMAT-5 (see captions to Fig. 2a) show that even insignificant temperature shifts using option 4·5 in COMAGMAT-5·2 (Supplementary Data Electronic Appendix 2, Part 4) may produce essential changes in the relative volumes of crystallized pyroxenes and olivine. The reasons for the strong effect of the modeled temperatures are discussed in Supplementary Data Electronic Appendix 2 (Part 5). APPLICATION OF COMAGMAT-5 TO THE BUSHVELD COMPLEX Below two examples of sulfide COMAGMAT use are given which demonstrate the genetic potential of the model, when it is applied to the sulfide saturation history of large layered intrusions and their parental magmas. We focus on modeling sulfide immiscibility in initially sulfide-undersaturated systems; these systems approximate the most primitive magmas and cumulates from the Basal Ultramafic Sequence of the Bushveld Complex in South Africa (Wilson, 2015). In addition, essential differences in the sulfide saturation history during equilibrium and fractional crystallization of the same ultramafic (UM) magma are considered. The Bushveld Complex and its parental magmas The Bushveld Complex was formed within a stable cratonic shield ∼2·06 Ga ago, as a result of three distinctive events; it covers an area up to 120 000 km2, including intrusions of mafic/ultramafic to granitic rocks into the Transvaal sedimentary sequences of the Kaapvaal craton (Naldrett et al., 2012). The first suite of rocks consists of early mafic sills, followed by felsic volcanic rocks forming the roof of a 7–8 km thick mafic to ultramafic cumulate succession, which represents the second intrusive suite, broadly referred to as the ‘Bushveld Complex’ or, more specifically, the so-called ‘Rustenburg Layered Suite’ (RLS; South African Committee for Stratigraphy, 1980). The radially inward dipping layered rocks cover an area of 40 000 km2 or even more (Cawthorn, 2015), thus comprising the world’s largest mafic/ultramafic intrusion with giant Ni–Cu–PGE, Cr, and Ti–V–Fe deposits (Eales & Costin, 2012). The third component of the Bushveld intrusive event is the Lebowa Granite Suite, which is possibly 2 km thick (Cawthorn, 2015). The detailed structure of the RLS has been examinated in five peripheral areas of the Bushveld Complex, including three (initially probably interconnected) basin-like lobes, the Eastern Bushveld, the Southeastern Bushveld and the Western Bushveld, as well as the Far Western Bushveld and the Northern Limb. A classical stratigraphic succession through the RLS succession is divided into five major zones, including the Marginal Zone composed of norite and minor pyroxenite (up to 800 m thick), the 800–1300 m thick ultramafic Lower Zone (mostly orthopyroxenite, harzburgite, and minor dunite), the Critical Zone marked by incoming chromitites in addition to orthopyroxenite (1·3–1·8 km thick), which is overlain by the ∼3·2 km thick Main Zone (norite, gabbros and minor anorthosite) and an as much as 2·8 km thick Upper Zone, composed of ferrogabbro and ferrodiorite (Naldrett et al., 2012; Cawthorn, 2015). During recent years, important progress has been made in constraining the compositional parameters of the most primitive Bushveld Complex magmas; see discussions in Wilson (2012, 2015), Yudovskaya et al. (2013, 2015) and Maier et al. (2016). Instead of the so-called B1 magma (12–13·5 wt % MgO, ∼56 wt % SiO2), which has long been considered as one probable parent for the Lower Zone and the Lower Critical Zone (Barnes et al., 2010; Wilson, 2012), much more primitive parental magmas, containing 19 wt % MgO (Wilson, 2015) and 18·7 wt% MgO (Maier et al., 2016) have been proposed. This is due to the inability of the B1 magma: (i) to produce the thick sequences of dunites and harzburgites in the Lower Zone, and (ii) a strong misbalance of Cr2O3 between the B1 composition and the extraordinarily high average Cr concentrations in the Bushveld rock sequences (Eales & Costin, 2012). Another contradiction is that the B1 magma was not capable of crystallizing olivine with Mg#>91 as observed in the Lower Zone of the Northern Limb and in the chilled sequence of the eastern Bushveld Complex (Wilson, 2012; Yudovskaya et al., 2013). These observations resulted in the conclusion that the B1 magma appears to be derived from a more mafic komatiitic parent; the composition of this magma can be approximated by the UM liquid from the Schwerin fold chill in the eastern Bushveld Complex, see Wilson (2015). Two model compositions were used in the example calculations (Table 2). The first one (denoted as ‘UM’) represents the most primitive UM liquid approximated by the Schwerin fold chill from the Basal Ultramafic Sequence of the Bushveld Complex, see the CH12/7&14UM column in table 4 of Wilson (2015). Focusing on S immiscibility calculations instead of modeling spinel crystallization, we postulated that the Cr2O3 concentration was zero and S concentration was 368 ppm in the UM melt. In fact, Wilson (2015) gave a much higher value 1914 ± 200 ppm S (based on analytical data) in the parent magma, whereas COMAGMAT-5·2 calculated SCSS = 1600 ppm at the UM liquidus temperature of 1456°C. This means that the ultramafic chill is locally over-saturated with sulfides, which contradicts previous conclusions on the sulfide under-saturated character of the B1 magma (Barnes et al., 2010; Ariskin et al., 2013), as a probable crystallization product of the more primitive UM magma (Wilson, 2015). Thus, the reason to reduce the initial S concentration in the initial melt is to avoid over-saturation of the UM magma with sulfide and sulfide immiscibility occurring too early in the modeled system. The second composition (denoted as ‘UM + 56 wt % Ol-91·6’, column 4 in Table 2) represents an approximation of a primitive cumulate rock, which is very similar to the olivine pyroxenite listed in the 5th column of Table 2. This composition is calculated as a mixture of 44 wt % of the proposed UM melt (the UM composition in column 3) with 56 wt% of olivine mg#91·6. We consider mg#91·6 to be a valid approximation of the original olivine composition in the Bushveld cumulates, because petrological reconstructions of the Lower Zone of the Northern Limb, the chilled sequence of the eastern Bushveld Complex and the Basal Ultramafic Sequence of the Bushveld Complex, indicate that the most primitive olivine varies in the range of mg#91–92 (Wilson, 2012, 2015; Yudovskaya et al., 2013). Note that the COMAGMAT-5 calculations for the UM magma yield the same liquidus olivine composition: mg#91·6 (see below). Sulfide immiscibility in the ultramafic magma and related cumulate Two series of plots in Figure 3 display the sulfide saturation history during equilibrium crystallization of two related systems, including the UM magma (Fig. 3a) and the proposed olivine enriched cumulate mush ‘UM + 56 wt % Ol-91·6’ (Fig. 3b). Both crystallization trajectories were calculated at QFM + 0·5 and anhydrous conditions, assuming P = 1 atm. The purpose of these calculations was not to replicate a slightly different order of crystallization at a more reliable pressure of 1–2 kbar (Wilson, 2015; Yudovskaya et al., 2015). Instead, we focused on marked differences in the onset of sulfide immiscibility in these samples at 1 atm., despite the fact that both systems were initially S-undersaturated, having almost the same initial melt (including S concentration) and olivine mg#91·6 at very similar temperatures (1456–1465°C), see Figure 3. Fig. 3. View largeDownload slide Modeled sequences of equilibrium crystallization, proportions of the modeled sulfide liquid with respect to precipitated solids, and compositional characteristics of crystallizing melts for two UM systems approximating a parental UM magma (A) and a related Ol-rich cumulate (B). Compositions of the UM parent (A) and the proposed ‘UM+56 wt.% Ol-91·6’ olivine cumulate (B) are listed in Table 2. Calculations using the COMAGMAT-5·2 model were carried out at QFM+0·5, anhydrous conditions and P=1 atm. The dashed green line in Figure 3b denotes the initial magma temperature consistent with ∼19 wt % MgO in the melt, which is in equilibrium with olivine mg#91·6 (Table 3). Light pink lines represent the evolution of S concentration in the modeled melts. Thick pink lines correspond to the sulfide liquid saturation of the modeled melt, when the sulfide solubility SCSS is equal to the S concentration in the melts. Note, that high NiO concentrations in the ‘UM+56 wt % Ol-91·6’ virtual melts at low percentage of crystallization are due to a large mass of Ni-containing olivine initially added to the UM parental melt. Fig. 3. View largeDownload slide Modeled sequences of equilibrium crystallization, proportions of the modeled sulfide liquid with respect to precipitated solids, and compositional characteristics of crystallizing melts for two UM systems approximating a parental UM magma (A) and a related Ol-rich cumulate (B). Compositions of the UM parent (A) and the proposed ‘UM+56 wt.% Ol-91·6’ olivine cumulate (B) are listed in Table 2. Calculations using the COMAGMAT-5·2 model were carried out at QFM+0·5, anhydrous conditions and P=1 atm. The dashed green line in Figure 3b denotes the initial magma temperature consistent with ∼19 wt % MgO in the melt, which is in equilibrium with olivine mg#91·6 (Table 3). Light pink lines represent the evolution of S concentration in the modeled melts. Thick pink lines correspond to the sulfide liquid saturation of the modeled melt, when the sulfide solubility SCSS is equal to the S concentration in the melts. Note, that high NiO concentrations in the ‘UM+56 wt % Ol-91·6’ virtual melts at low percentage of crystallization are due to a large mass of Ni-containing olivine initially added to the UM parental melt. Modeling the UM crystallization yields: Ol (1456°C, mg#91·6) → Ol (1222°C, mg#84·7) + Opx (mg#87·0) → Ol (1185°C, mg#82·7) + Opx (mg#85·2) + Sulfide liquid → Ol (1178°C, mg#82·3) + Opx (mg#84·9) + Sulfide liquid + Pl (An72·7) → Ol (1152°C, mg#78·6) + Opx (mg#82·2) + Sulfide liquid + Pl (An66·5) + Aug (mg#82·4) → plus ilmenite (at 1116°C) and Ti–magnetite (at 1105°C). Changes in the crystallizing mineral assemblages are shown in Figure 3a in terms of the percent crystallized. The ‘UM + 56 wt % Ol-91·6’ cumulate melt yields a similar order of crystallization, only starting at much higher temperatures: Ol (1724°C, mg#95·8) → Ol (1289°C, mg#89·5) + sulfide liquid → Ol (1231°C, mg#89·0) + Opx (mg#90·6) + sulfide liquid → Ol (1194°C, mg#88·6) + Opx (mg#90·1) + sulfide liquid + Pl (An73·8) → Ol (1165°C, mg#87·8) + Opx (mg#89·7) + sulfide liquid + Pl (An65·8) + Aug (mg#89·7) → plus ilmenite (at 1134°C and lower), see Figure 3b. In the following interpretation, it is noteworthy that neither the calculated maximum liquidus temperature (1724°C), nor the modeled liquidus olivine (mg#95·8), represent realistic estimates of initial magma parameters. Instead, the results of calculations at 1724°C should be considered as virtual characteristics of the 100% melted olivine rich ‘model cumulate’. The same is true for other calculations in the fictitious range 1724°C to ∼1465°C, with the latter value corresponding to the proposed two-phase mixture of a magmatic melt (19 wt % MgO) and olivine mg#91·6. Thus, a petrological interpretation only considers results obtained below the proposed UM magma temperature (approximately ≤1460°C); see Figure 3a and the dashed line in Figure 3b. The main differences between the ‘UM’ and ‘UM + 56 wt % Ol-91·6’ calculations within the apparent range of magma crystallization temperatures (∼1460–1120°C) are recorded in the onset of modeled sulfide immiscibility. The first sulfide liquid was calculated to occur at 1185°C in the case of the parental UM magma and at the much higher 1289°C temperature for the modeled olivine cumulate. We consider these differences as one more demonstration of the coupled compositional effect of FeO and Ni on sulfide solubility (SCSS), because the olivine cumulate contained approximately three times more NiO than the pure UM melt (Table 3), see details in Discussion. Table 3: Parameters of equilibrium and fractional crystallization for the proposed UM magma and modeled cumulate ‘UM + 56 wt.% Ol-91·6’, see Table 2 Melt components, wt.% Modeled Ol cumulate ‘UM + 56 wt.% Ol-91·6’ Proposed UM magma Crystallization Equilibrium Equilibrium Fractional Bulk NiO in the modeled systems  NiO, wt.% 0·317 0·090 Parameters of magmatic melt in equilibrium with olivine mg#91·6  Temp, °C 1465 1456  MgO, wt.% 19·4 18·9  FeO, wt.% 10·39 10·14  NiO, wt.% 0·114 0·090  Sinit, wt.% 0·036 0·037  SCSS, wt.% 0·163 0·160  mg# in Ol,% 91·6 91·6  NiO in Ol, wt. % 0·479 0·393  Ol, wt.% 44·5 0 The onset of sulfide immiscibility Mineral assemblage Ol+Sulfide liquid Ol+OPx+Sulfide liquid Pl+Pig+Sulfide liquid  Temp, °C 1289 1185 1116  MgOmelt, wt.% 10·22 6·66 2·93  FeO, wt.% 7·08 8·34 11·40  NiOmelt, wt.% 0·039 0·015 0·002  S (in melt) = SCSS, wt.% 0·050 0·062 0·086  mg# in Ol, % 89·5 82·7 –  NiO in Ol, wt.% 0·444 0·290 –  Ni/(Ni+Fe) in sulfide, atomic 0·410 0·318 0·237 Melt components, wt.% Modeled Ol cumulate ‘UM + 56 wt.% Ol-91·6’ Proposed UM magma Crystallization Equilibrium Equilibrium Fractional Bulk NiO in the modeled systems  NiO, wt.% 0·317 0·090 Parameters of magmatic melt in equilibrium with olivine mg#91·6  Temp, °C 1465 1456  MgO, wt.% 19·4 18·9  FeO, wt.% 10·39 10·14  NiO, wt.% 0·114 0·090  Sinit, wt.% 0·036 0·037  SCSS, wt.% 0·163 0·160  mg# in Ol,% 91·6 91·6  NiO in Ol, wt. % 0·479 0·393  Ol, wt.% 44·5 0 The onset of sulfide immiscibility Mineral assemblage Ol+Sulfide liquid Ol+OPx+Sulfide liquid Pl+Pig+Sulfide liquid  Temp, °C 1289 1185 1116  MgOmelt, wt.% 10·22 6·66 2·93  FeO, wt.% 7·08 8·34 11·40  NiOmelt, wt.% 0·039 0·015 0·002  S (in melt) = SCSS, wt.% 0·050 0·062 0·086  mg# in Ol, % 89·5 82·7 –  NiO in Ol, wt.% 0·444 0·290 –  Ni/(Ni+Fe) in sulfide, atomic 0·410 0·318 0·237 Table 3: Parameters of equilibrium and fractional crystallization for the proposed UM magma and modeled cumulate ‘UM + 56 wt.% Ol-91·6’, see Table 2 Melt components, wt.% Modeled Ol cumulate ‘UM + 56 wt.% Ol-91·6’ Proposed UM magma Crystallization Equilibrium Equilibrium Fractional Bulk NiO in the modeled systems  NiO, wt.% 0·317 0·090 Parameters of magmatic melt in equilibrium with olivine mg#91·6  Temp, °C 1465 1456  MgO, wt.% 19·4 18·9  FeO, wt.% 10·39 10·14  NiO, wt.% 0·114 0·090  Sinit, wt.% 0·036 0·037  SCSS, wt.% 0·163 0·160  mg# in Ol,% 91·6 91·6  NiO in Ol, wt. % 0·479 0·393  Ol, wt.% 44·5 0 The onset of sulfide immiscibility Mineral assemblage Ol+Sulfide liquid Ol+OPx+Sulfide liquid Pl+Pig+Sulfide liquid  Temp, °C 1289 1185 1116  MgOmelt, wt.% 10·22 6·66 2·93  FeO, wt.% 7·08 8·34 11·40  NiOmelt, wt.% 0·039 0·015 0·002  S (in melt) = SCSS, wt.% 0·050 0·062 0·086  mg# in Ol, % 89·5 82·7 –  NiO in Ol, wt.% 0·444 0·290 –  Ni/(Ni+Fe) in sulfide, atomic 0·410 0·318 0·237 Melt components, wt.% Modeled Ol cumulate ‘UM + 56 wt.% Ol-91·6’ Proposed UM magma Crystallization Equilibrium Equilibrium Fractional Bulk NiO in the modeled systems  NiO, wt.% 0·317 0·090 Parameters of magmatic melt in equilibrium with olivine mg#91·6  Temp, °C 1465 1456  MgO, wt.% 19·4 18·9  FeO, wt.% 10·39 10·14  NiO, wt.% 0·114 0·090  Sinit, wt.% 0·036 0·037  SCSS, wt.% 0·163 0·160  mg# in Ol,% 91·6 91·6  NiO in Ol, wt. % 0·479 0·393  Ol, wt.% 44·5 0 The onset of sulfide immiscibility Mineral assemblage Ol+Sulfide liquid Ol+OPx+Sulfide liquid Pl+Pig+Sulfide liquid  Temp, °C 1289 1185 1116  MgOmelt, wt.% 10·22 6·66 2·93  FeO, wt.% 7·08 8·34 11·40  NiOmelt, wt.% 0·039 0·015 0·002  S (in melt) = SCSS, wt.% 0·050 0·062 0·086  mg# in Ol, % 89·5 82·7 –  NiO in Ol, wt.% 0·444 0·290 –  Ni/(Ni+Fe) in sulfide, atomic 0·410 0·318 0·237 Fractional vs. equilibrium crystallization The crystallization sequence shown in Figure 4 represents the results of modeling perfect fractionation of the same UM magma shown in Figure 3a. During the early stages in the temperature range ∼1460–1200°C, modeling equilibrium (Fig. 3a) and fractional (Fig. 4) crystallization results in the same order of the two first minerals to crystallize, olivine → orthopyroxene, but the late stages differ considerably. First, note the much later sulfide immiscibility, which is observed at ∼70°C lower than the temperature of sulfide precipitation onset during equilibrium crystallization (Table 3). As far as orthopyroxene during fractionation was replaced by pigeonite, sulfide saturation is contemporaneous with the crystallization of a Pl + Pig gabbroic assemblage, followed by Pl + Pig + Aug + Mt (Fig. 4). Other sulfide saturation parameters during UM parental magma fractional crystallization are given in Table 3. Fig. 4. View largeDownload slide The order of fractional crystallization of the UM parental magma (Table 2) and comparison of S vs. SCSS trends and NiO in the melt in the case of fractional and equilibrium crystallization. Calculations by the COMAGMAT-5·2 model were carried at QFM+0·5, anhydrous conditions, and P=1 atm. The solid gray line and the solid black line represent changes in S and NiO contents in the melt in the case of equilibrium and fractional crystallization, respectively. The dotted lines display the evolution of the SCSS values for both crystallization trajectories. The onset of sulfide-silicate immiscibility is consistent with the intersection of the S and SCSS trends for each type of crystallization. The light and bold pink lines denote intervals of crystallization under sulfide saturated conditions. Fig. 4. View largeDownload slide The order of fractional crystallization of the UM parental magma (Table 2) and comparison of S vs. SCSS trends and NiO in the melt in the case of fractional and equilibrium crystallization. Calculations by the COMAGMAT-5·2 model were carried at QFM+0·5, anhydrous conditions, and P=1 atm. The solid gray line and the solid black line represent changes in S and NiO contents in the melt in the case of equilibrium and fractional crystallization, respectively. The dotted lines display the evolution of the SCSS values for both crystallization trajectories. The onset of sulfide-silicate immiscibility is consistent with the intersection of the S and SCSS trends for each type of crystallization. The light and bold pink lines denote intervals of crystallization under sulfide saturated conditions. Two pairs of curves to the right in Figure 4 display changes in the modeled S concentrations and SCSS in the melt, allowing insights into the differences in the onset of sulfide–silicate immiscibility between fractional and equilibrium crystallization. To interpret the modeled relationships, remember that our approach to simulating sulfide saturation is similar to that proposed by Li & Ripley (2005) and utilized by other workers (e.g. Barnes, 2007; Barnes et al., 2010; Ariskin et al., 2016). The evolution of melt S content in the case of equilibrium (the solid grey line) and fractional (the solid black line) crystallization was calculated assuming an initial S content of 368 ppm in the UM parent (Table 2). The sulfide solubility values (SCSS, dotted lines) were calculated as a function of varying temperature and melt composition, including Ni concentrations in the derivative melts (Ariskin et al., 2013). The point of S and SCSS lines intersection (i.e. S = SCSS) shows the onset of sulfide-silicate immiscibility, thus giving rise to precipitation of a Fe–Ni sulfide liquid, so that the S content in the melt follows the SCSS trend (Fig. 4). The modeled relationships demonstrate that higher sulfide solubility in the case of the fractional crystallization of the UM parent (compare SCSS values at the same temperature) should be correlated with an increase in FeO contents and much lower NiO in the melt, due to fractionation of olivine and low-Ca pyroxenes (Opx/Pig). Finally, this gives rise to much later sulfide saturation in the modeled system as compared to the course of equilibrium crystallization of the same magmatic melt (compare columns 2 & 3 in Table 3). DISCUSSION Problems of thermodynamics and the use of COMAGMAT-5 G-minimization vs. the mass action law As is evident from the thermodynamic considerations discussed here (see Eqs. (1–12)), despite its semi-empirical character, COMAGMAT was designed to follow fundamental thermodynamic principles. The main difference between its basic algorithm and those used in MELTS/pMELTS is that the algorithms used in the latter models directly minimize the isobaric potential at a given temperature (Ghiorso, 1987, 1994), whereas design of the COMAGMAT models allows the equilibrium state to be searched for at a pre-defined crystallization degree. In this case, the equilibrium problem is approached by means of the iterative solution of nonlinear equations consistent with the law of mass action (Eqs. (4, 5)). It is noteworthy that methods of calculating the products of chemical reactions based on the mass action law were developed in chemistry long before Gibbs formulated his thermodynamic theory. Thus, there is no principal conflict between these two methods, because the constants of heterogeneous equilibria are linked to the changes in the free energy through the fundamental equation ln K=-ΔG/RgT. In this sense, it would be possible to transform the energy functions used by Ghiorso and his followers to the ‘language of equilibrium constants’ describing mineral–melt equilibria. Their values (more accurately lnK vs. 1/T dependences) could then be used in the COMAGMAT algorithm to solve the equilibrium problem in the same manner as is done using the empirically calibrated mineral–melt geothermometers. Approaching the effect of low pressures Because all mineral–melt geothermometers used in the COMAGMAT-5 model were calibrated on experimental data obtained at 1 atm., formally, this model should be used only for ‘atmospheric’ calculations. However, based on our experience in modeling crystallization at elevated pressures (Ariskin et al., 1990; Ariskin & Barmina, 2004), we have slightly extended the possible range of pressures, up to about 1–2 kbar. There are two reasons for such approximations. First, the standard deviations of olivine and plagioclase thermometric calculations (±10–15°C) are two to three times higher than the effect of pressure on their liquidus temperatures (∼5°C/kbar). For pyroxenes this effect is known to be ∼10–12°C/kbar, which is also within the accuracy of the proposed pyroxene–melt geothermometers. Second, an increase in total pressure results in a very similar increase in the liquidus temperatures for olivine and plagioclase; therefore, the chemical evolution of the modeled melts crystallizing at slightly elevated pressures barely deviates from that obtained at 1 atm., at least up to the onset of pyroxene crystallization (Ariskin et al., 1990). Admitedly in the case of modeling tholeiitic or other systems, which crystallize olivine and plagioclase prior to crystallizing clinopyroxene, the effect of small pressure changes on the modeled crystallization of troctolitic assemblages should be negligible. This supports the application of COMAGMAT-5 to shallow magma chambers occurring at depths of no more 6–7 km, particularly for Ol±Pl saturated magmas similar to MORB (Ariskin et al., 2016). The effect of Ni on SCSS Using the sulfide COMAGMAT-5, it was shown (Ariskin et al., 2013) that despite relatively low concentrations, Ni has a pronounced lowering effect on S solubility, causing a significant temperature increase in the onset of sulfide immiscibility in melts with similar major element compositions. Despite the fact that this hypothesis has found some support in natural observations (e.g. Barnes et al., 2013), it continues to be debated. Recently Fortin et al. (2015) presented a new FeS–sulfide solubility equation (his so-called ‘Model-B’), which reproduced S contents both for their calibration dataset and for 65 natural Ni-containing sulfide-saturated MORB glasses within about 5%, even without considering the effect of Ni. It was concluded that ‘there is no need for a more complicated multi-species FeNiS model to calculate the SCSS, even in high-Ni, natural samples’ (Fortin et al., 2015, p. 113). In our study, we present an additional test of the Ni-effect on a magnesian Ni-rich sulfide-saturated glass dredged from the the Southern Atlantic. These results clearly demonstrate the strong effect of Ni on the modeled sulfide solubility (Table 1). This outcome seems to be of fundamental importance, particularly when the SCSS calculations are applied to the history of sulfide saturation in picritic magmas and mafic to ultramafic layered intrusions. This is because of the higher Ni contents in the high-temperature melt and the possible presence of a large amount of Ni-rich olivine in the parental magmas, so that the bulk composition of the ultramafic system has to control the sulfide immiscibility process during both magma crystallization and solidification of intercumulus melt in igneous cumulates (Table 3). CONCLUSIONS A detailed description of the thermodynamic framework and main options of the first ‘sulfide version’ of the COMAGMAT-5 magma crystallization model is presented. The new features include: (i) complete recalibration of silicate mineral– melt geothermometers; (ii) the capability to calculate the Ni content of the modeled femic minerals and sulfides; and (iii) incorporation of a recently proposed Fe–Ni sulfide solubility model (Ariskin et al., 2013) into a general algorithm of modeling crystallization for S-containing magmas. Due to updating of the basic COMAGMAT algorithm, the new model allows new geothermometers to be used for simulations of co-crystallization Aug±Pig±Opx, including peritectic relations between these phases and olivine. Combining these calculations with the proposed SCSS model, COMAGMAT-5 can be used to simulate equilibrium and fractional crystallization of S-saturated and S-undersaturated mafic to UM magmas, including changes in sulfide-silicate (±Fe–Ti oxides) proportions for multiply-saturated mineral assemblages. The COMAGMAT-5 model (ver. 5·2) is available for downloading from http://geo.web.ru/∼ariskin/soft.html-id=comagmat.htm. Following Ariskin et al. (2013), the strong effect of Ni on SCSS modeled in mafic to ultramafic systems must be taken into account, so that at higher Ni contents in the melt an essential decrease in the sulfide solubility should be expected. This statement was tested against a Ni-rich sulfide-saturated glass dredged from the southern Mid-Atlantic Ridge near the Bouvet Triple Junction (Kamenetsky et al., 2001, 2013). Applying COMAGMAT-5 to the BTJ glass, we obtained SCSS from 514 ppm (QFM) to 710 ppm (QFM + 1), close to the 600 ppm observed (Kamenetsky et al., 2013). The use of virtual ‘FeS’ solubility models results in much higher sulfide solubility, mostly > 1000 ppm (Li & Ripley, 2009; Fortin et al., 2015), even assuming the strong dampening effect of low H2O contents on SCSS postulated in Liu et al. (2007). Example calculations using compositional proxies of a UM parental melt and a related olivine cumulate from the Basal Ultramafic Sequence of the Bushveld Complex (see UM in Table 2; Wilson, 2015) demonstrate a high temperature interval in the modeled trajectories where sulfide–silicate immiscibility could occur. The highest 1289°C temperature (Ol + sulfide) was observed in crystallizing olivine cumulate, an intermediate 1185°C temperature (Ol + Opx + sulfide) was calculated in the case of equilibrium UM magma crystallization, and the lowest 1116°C temperature (Pl + Pig + sulfide) was obtained as a result of perfect fractionation of the same melt (Table 3). The modeled effect seems to be of general genetic significance, particularly for the study of the formation of sulfide-saturated volcanic suites and fertile layered intrusions. ACKNOWLEDGEMENTS We gratefully acknowledge thoughtful reviews by Sébastien Jégo, Guil Gualda, Gokce Ustunisik, Paula Antoshechkina and two anonymous reviewers, as well as constructive comments from editor Wendy Bohrson. We thank Leonid Danyushevsky for inspiring discussions, Jun-Ichi Kimura for his help with MELTS/pMELTS calculations, and Candace S. O'Connor for careful editing of the final version of the manuscript. FUNDING We acknowledge support of AngloAmerican, BHP Billiton, Votorantim Metais (through AMIRA project P962, 2007–2010) and the Australian Research Council funding to the CODES (Hobart, Australia) at the initial stage of development of the COMAGMAT model. Final calibrations, testing and example calculations using COMAGMAT-5 are supported by the Russian Science Foundation (project RSF 16–17-10129). SUPPLEMENTARY DATA Supplementary data for this paper are available at Journal of Petrology online. REFERENCES Almeev R. R. , Ariskin A. A. , Kimura J.-I. , Barmina G. S. ( 2013 ). The role of polybaric crystallization in genesis of andesitic magmas: phase equilibria simulations of Bezymianny volcanic subseries . Journal of Volcanology and Geothermal Research 263 , 182 – 192 . Google Scholar CrossRef Search ADS Ariskin A. A. ( 1999 ). Phase equilibria modeling in igneous petrology: use of COMAGMAT model for simulating fractionation of ferro-basaltic magmas and the genesis of high-alumina basalt . Journal of Volcanology and Geothermal Researches 90 , 115 – 162 . Google Scholar CrossRef Search ADS Ariskin A. A. ( 2003 ). The compositional evolution of differentiated liquids from the Skaergaard layered series as determined by geochemical thermometry . Russian Journal of Earth Sciences 5 , 1 – 29 . Google Scholar CrossRef Search ADS Ariskin A. A. , Barmina G. S. ( 1990 ). Equilibria thermometry between plagioclases and basalt or andesite magmas . Geochemistry International 27 , 129 – 134 . Ariskin A. A. , Barmina G. S. ( 2004 ). COMAGMAT: development of a magma crystallization model and its petrologic applications . Geochemistry International 42 (Suppl 1) , S1 – S157 . Ariskin A. A. , Yaroshevsky A. A. ( 2006 ). Crystallization differentiation of intrusive magmatic melt: development of a Convection-Accumulation model . Geochemistry International 44 , 72 – 93 . Google Scholar CrossRef Search ADS Ariskin A. A. , Frenkel M. Y. , Tsekhonya T. I. ( 1990 ). High-pressure fractional crystallization of tholeiitic magmas . Geochemistry International 27 , 10 – 20 . Ariskin A. A. , Barmina G. S. , Bychkov K. A. , Danyushevsky L. V. ( 2009 ). Parental magmas of mafic layered intrusions: using an updated COMAGMAT model for calculations of sulfide-silicate cotectics during their crystallization . Northwestern Geology 42 , 1 – 3 . Ariskin A. A. , Barmina G. S. , Frenkel M. Ya. , Nielsen R. L. ( 1993 ). COMAGMAT: a Fortran program to model magma differentiation processes . Computers and Geosciences 19 , 1155 – 1170 . Google Scholar CrossRef Search ADS Ariskin A. A. , Danyushevsky L. V. , Bychkov K. A. , McNeill A. W. , Barmina G. S. , Nikolaev G. S. ( 2013 ). Modeling solubility of Fe-Ni sulfides in basaltic magmas: the effect of Ni in the melt . Economic Geology 108 , 1983 – 2003 . Google Scholar CrossRef Search ADS Ariskin A. A. , Kislov E. V. , Danyushevsky L. V. , Nikolaev G. S. , Fiorentini M. L. , Gilbert S. , Goemann K. , Malyshev A. ( 2016 ). Cu-Ni-PGE fertility of the Yoko-Dovyren layered massif (Northern Transbaikalia, Russia): thermodynamic modeling of sulfide compositions in low mineralized dunites based on quantitative sulfide mineralogy . Mineralium Deposita 51 , 993 – 1011 . Google Scholar CrossRef Search ADS Barnes S. J. ( 2007 ). Cotectic precipitation of olivine and sulfide liquid from komatiite magma and the origin of komatiite-hosted disseminated nickel sulfide mineralization at Mount Keith and Yakabindie, Western Australia . Economic Geology 106 , 298 – 304 . Barnes S.-J. , Lightfoot P. C. ( 2005 ). Formation of magmatic nickel sulfide deposits and processes affecting their copper and Platinum Group Element contents . Economic Geology 100th Anniversary Volume, 179 – 213 . Barnes S.-J. , Maier W. D. , Curl E. A. ( 2010 ). Composition of the marginal rocks and sills of the Rustenburg Layered Suite, Bushveld Complex, South Africa: implications for the formation of the Platinum-Group element deposits . Economic Geology 105 , 1491 – 1511 . Google Scholar CrossRef Search ADS Barnes S. J. , Godel B. , Gürer D. , Brenan J. M. , Robertson J. , Paterson D. ( 2013 ). Sulfide-olivine Fe-Ni exchange and the origin of anomalously Ni rich magmatic sulfides . Economic Geology 108 , 1971 – 1982 . Google Scholar CrossRef Search ADS Berman R. G. ( 1988 ). Internally-consistent thermodynamic data for minerals in the system Na2O-K2O-CaO-MgO-FeO-Fe2O3-Al2O3-SiO2-TiO2-H2O-CO2 . Journal of Petrology 29 , 445 – 522 . Google Scholar CrossRef Search ADS Campbell I. H. , Naldrett A. J. ( 1979 ). The influence of silicate: sulfide ratios on the geochemistry of magmatic sulfides . Economic Geology 74 , 1503 – 1506 . Google Scholar CrossRef Search ADS Cawthorn R. G. ( 2015 ). The Bushveld Complex, South Africa. In: Charlier B. , Namur O. , Latypov R. , Tegner C. (eds) Layered Intrusions . Berlin : Springer , pp. 517 – 587 . Google Scholar CrossRef Search ADS Danyushevsky L. V. , Plechov P. ( 2011 ). Petrolog3: integrated software for modeling crystallization processes . Geochemistry, Geophysics, Geosystems 12 , Q07021 . Google Scholar CrossRef Search ADS Eales H. V. , Costin G. ( 2012 ). Crustally contaminated komatiite: primary source of the chromitites and Marginal, Lower, and Critical zone magmas in a staging chamber beneath the Bushveld Complex . Economic Geology 107 , 645 – 665 . Google Scholar CrossRef Search ADS Falloon T. J. , Danyushevsky L. V. , Ariskin A. , Green D. H. , Ford C. E. ( 2007 ). The application of olivine geothermometry to infer crystallization temperatures of parental liquids: Implications for the temperature of MORB magmas . Chemical Geology 241 , 207 – 233 . Google Scholar CrossRef Search ADS Ford C. , Russell D. , Craven J. , Fisk M. ( 1983 ). Olivine-liquid equilibria: temperature, pressure and composition dependence of the crystal/liquid cation partition coefficients for Mg, Fe2+, Ca and Mn . Journal of Petrology 24 , 256 – 265 . Google Scholar CrossRef Search ADS Fortin M.-A. , Riddle J. , Desjardins-Langlais Y. , Baker D. R. ( 2015 ). The effect of water on the sulfur concentration at sulfide saturation (SCSS) in natural melts . Geochimica et Cosmochimica Acta 160 , 100 – 116 . Google Scholar CrossRef Search ADS Frenkel M. Y. , Ariskin A. A. ( 1984 ). A computer algorithm for equilibration in a crystallizing basalt magma . Geochemistry International 21 , 63 – 73 . Frenkel M. Y. , Yaroshevsky A. A. , Ariskin A. A. , Barmina G. S. , Koptev-Dvornikov E. V. , Kireev B. S. ( 1989 ). Convective-cumulative model simulating the formation process of stratified intrusions. In: Bonin B. , Didier J. , Le Fort P. , Propach G. , Puga E. , Vistelius A. B. (eds) Magma-Crust Interactions and Evolution . Athens-Greece : Theophrastus Publications SA , pp. 3 – 88 . Ghiorso M. S. ( 1987 ). Modeling magmatic systems: thermodynamic relations. In: Carmichael I. S. E. , Eugster H. P. (eds) MSA Reviews in Mineralogy 17: Thermodynamic Modeling of Geological Materials: Minerals, Fluids and Melts, 443–465. Ghiorso M. S. ( 1994 ). Algorithms for the estimation of phase stability in heterogeneous thermodynamic systems . Geochimica et Cosmochimica Acta 58 , 5489 – 5501 . Google Scholar CrossRef Search ADS Ghiorso M. S. ( 1997 ). Thermodynamic models of igneous processes . Annual Review of Earth and Planetary Sciences 25 , 221 – 241 . Google Scholar CrossRef Search ADS Ghiorso M. S. , Sack R. O. ( 1995 ). Chemical mass transfer in magmatic processes: IV. A revised and internally consistent thermodynamic model for the interpolation and extrapolation of liquid–solid equilibria in magmatic systems at elevated temperatures and pressures . Contributions to Mineralogy and Petrology 119 , 197 – 212 . Google Scholar CrossRef Search ADS Ghiorso M. S. , Hirschmann M. M. , Reiners P. W. , Kress V. C. III . ( 2002 ). The pMELTS: a revision of MELTS for improved calculation of phase relations and major element partitioning related to partial melting of the mantle to 3 GPa . Geochemistry, Geophysics, Geosystems 3 , doi:10.1029/2001GC000217. Gongalsky B. I. , Krivolutskaya N. A. , Ariskin A. A. , Nikolaev G. S. ( 2016 ). The Chineysky gabbronorite-anorthosite layered massif (Northern Transbaikalia, Russia): its structure, Fe-Ti-V and Cu-PGE deposits, and parental magma composition . Mineralium Deposita 51 , 1013 – 1034 . Google Scholar CrossRef Search ADS Grove T. L. , Elkins-Tanton L. T. , Parman S. W. , Chatterjee N. , Muntener O. , Gaetani G. A. ( 2003 ). Fractional crystallization and mantle-melting controls on calc-alkaline differentiation trends . Contributions to Mineralogy and Petrology 145 , 515 – 533 . Google Scholar CrossRef Search ADS Gualda G. A. R. , Ghiorso M. S. ( 2015 ). MELTS_Excel: a Microsoft Excel-based MELTS interface for research and teaching of magma properties and evolution . Geochemistry, Geophysics, Geosystems 16 , 315 – 324 . Google Scholar CrossRef Search ADS Gualda G. A. R. , Ghiorso M. S. , Lemons R. V. , Carley T. L. ( 2012 ). Rhyolite-MELTS: a modified calibration of MELTS optimized for silicarich, fluid-bearing magmatic systems . Journal of Petrology 53 , 875 – 890 . Google Scholar CrossRef Search ADS Jugo P. J. ( 2009 ). Sulfur content at sulfide saturation in oxidized magmas . Geology 37 , 415 – 418 . Google Scholar CrossRef Search ADS Kamenetsky V. S. , Maas R. , Sushchevskaya N. M. , Norman M. D. , Cartwright I. , Peyve A. A. ( 2001 ). Remnants of Gondwanan continental lithosphere in oceanic upper mantle: evidence from the South Atlantic Ridge . Geology 29 , 243 – 246 . Google Scholar CrossRef Search ADS Kamenetsky V. S. , Maas R. , Fonseca R. O. C. , Ballhaus C. , Heuser A. , Brauns M. , Norman M. D. , Woodhead J. D. , Rodemann T. , Kuzmin D. V. , Bonatti E. ( 2013 ). Noble metals potential of sulfide-saturated melts from the subcontinental lithosphere . Geology 41 , 575 – 578 . Google Scholar CrossRef Search ADS Karpov I. K. , Chudnenko K. V. , Kulik D. A. , Avchenko O. V. , Bychinskii V. A. ( 2001 ). Minimization of Gibbs free energy in geochemical systems by convex programming . Geochemistry International 39 , 1108 – 1119 . Li C. , Ripley E. M. ( 2005 ). Empirical equations to predict the sulfur content of mafic magmas at sulfide saturation and applications to magmatic sulfide deposits . Mineralium Deposita 40 , 218 – 230 . Google Scholar CrossRef Search ADS Li C. , Ripley E. M. ( 2009 ). Sulfur contents at sulfide-liquid or anhydrite saturation in silicate melts: empirical equations and example applications . Economic Geology 104 , 405 – 412 . Google Scholar CrossRef Search ADS Liu Y. , Samaha N.-T. , Baker D. R. ( 2007 ). Sulfur concentration at sulfide saturation (SCSS) in magmatic silicate melts . Geochimica et Cosmochimica Acta 71 , 1783 – 1799 . Google Scholar CrossRef Search ADS Maier W. D. , Barnes S.-J. , Karykowski B. T. ( 2016 ). A chilled margin of komatiite and Mg-rich basaltic andesite in the western Bushveld Complex, South Africa . Contributions to Mineralogy and Petrology 171 , 57 . Google Scholar CrossRef Search ADS Mathez E. A. ( 1976 ). Sulfur solubility and magmatic sulfides in submarine basaltic glass . Journal of Geophysical Research 81 , 4269 – 4275 . Google Scholar CrossRef Search ADS McBirney A. R. ( 1996 ). The Skaergaard intrusion. In: Cawthorn, R. G. (ed.) Developments in Petrology , vol. 15: Layered Intrusions , pp. 147 – 180 . Naldrett A. J. , Wilson A. , Kinnaird J. , Yudovskaya M. , Chunnett G. ( 2012 ). The origin of chromitites and related PGE mineralization in the Bushveld Complex: new mineralogical and petrological constraints . Mineralium Deposita 47 , 209 – 232 . Google Scholar CrossRef Search ADS Nielsen R. L. ( 1990 ). Simulation of igneous differentiation processes. In: Nicholls J. , Russell J. K. (eds) Modern Methods of Igneous Petrology: Understanding Magmatic Processes. Reviews in Mineralogy 24 , 63 – 105 . Nielsen R. L. , Dungan M. A. ( 1983 ). Low-pressure mineral-melt equilibria in natural anhydrous mafic systems . Contributions to Mineralogy and Petrology 84 , 310 – 326 . Google Scholar CrossRef Search ADS Nikolaev G. S. , Ariskin A. A. , Barmina G. S. , Nazarov M. A. , Almeev R. R. ( 2016 ). Test of the Ballhaus–Berry–Green Ol–Opx–Sp oxybarometer and calibration of a new equation for estimating the redox state of melts saturated with olivine and spinel . Geochemistry International 54 , 301 – 320 . Google Scholar CrossRef Search ADS SACS (South African Committee for Stratigraphy) . ( 1980 ). Kent L. E. (Compiler) Stratigraphy of South Africa . Vol. 8 . Pretoria : Geological Survey of South Africa, Handbook , p. 690 . Shvarov Y. V. ( 1981 ). A general equilibrium criterion for an isobaric-isothermal model of a chemical system . Geochemistry International 18 , 38 – 45 . Wilson A. H. ( 2012 ). A chill sequence to the Bushveld Complex: insight into the first stage of emplacement and implications for the parental magmas . Journal of Petrology 53 , 1123 – 1168 . Google Scholar CrossRef Search ADS Wilson A. H. ( 2015 ). The earliest stages of emplacement of the Eastern Bushveld Complex: development of the lower zone, marginal zone and basal ultramafic sequence . Journal of Petrology 56 , 347 – 388 . Google Scholar CrossRef Search ADS Yudovskaya M. A. , Naldrett A. J. , Woolfe J. A. S. , Costin G. , Kinnaird J. A. ( 2015 ). Reverse compositional zoning in the Uitkomst chromitites as an indication of crystallization in a magmatic conduit . Journal of Petrology 56 , 2373 – 2394 . Google Scholar CrossRef Search ADS Yudovskaya M. A. , Kinnaird J. A. , Sobolev A. V. , Kuzmin D. V. , McDonald I. , Wilson A. H. ( 2013 ). Petrogenesis of the Lower Zone olivine-rich cumulates beneath the Platreef and their correlation with recognized occurrences in the Bushveld Complex . Economic Geology 108 , 1923 – 1952 . Google Scholar CrossRef Search ADS © The Author(s) 2018. Published by Oxford University Press. All rights reserved. For permissions, please e-mail: journals.permissions@oup.com This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/about_us/legal/notices)

Journal

Journal of PetrologyOxford University Press

Published: Mar 1, 2018

There are no references for this article.

You’re reading a free preview. Subscribe to read the entire article.


DeepDyve is your
personal research library

It’s your single place to instantly
discover and read the research
that matters to you.

Enjoy affordable access to
over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month

Explore the DeepDyve Library

Search

Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly

Organize

Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place.

Access

Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals.

Your journals are on DeepDyve

Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more.

All the latest content is available, no embargo periods.

See the journals in your area

DeepDyve

Freelancer

DeepDyve

Pro

Price

FREE

$49/month
$360/year

Save searches from
Google Scholar,
PubMed

Create lists to
organize your research

Export lists, citations

Read DeepDyve articles

Abstract access only

Unlimited access to over
18 million full-text articles

Print

20 pages / month

PDF Discount

20% off