A New Iterative Method to Evaluate Fractional Crystallization in Subalkalic Basalts

Yucheng Pan and Rodey Batiza
Department of Geology and Geophysics, SOEST
University of Hawaii, 1680 East-West Rd
Honolulu, HI 96822

Published previously in abstract form as:  Pan, Y., The application of Microsoft Excel "Solver" to VGP curricula: a handy iteration math tool, EOS Trans AGU, 78 (46), Fall Meeting suppl., F667, 1997;  Pan, Y. and Batiza, R., A new iterative method to evaluate fractional crystallization in subalkalic basalts, Penrose Conference Series: Evolution of Ocean Island volcanoes, 69, 1998.

Detailed instructions to perform the calculation described below are included in the Excel application spreadsheet "PB" which can be used on both Macintosh and Windows computers.

1. Introduction

Knowing the compositions and abundances of solid phases that crystallize from basalt during cooling is useful for estimating the conditions of crystallization, and the thermal structure of the crust and upper mantle. Further, such information can be used to correct basalt compositions to evaluate melting and other deeper processes. We describe a direct, iterative method for estimating the compositions and proportions of fractionating solids that can be used as a complement to the indirect least-squares methods and calibrated liquid line of descent models (e.g. Ghiorso and Sack, 1995; Weaver and Langmuir, 1990).

This method iteratively solves a multi-component, multiphase mass balance to achieve a "best fit". Our method was developed for use in subalkalic basalt systems at relatively high temperature, and considers only the components Si, Al, Ti, Na, and (Fe+Mg+Ca) in solid phases olivine (ol), plagioclase (pl), and high-Ca pyroxene (cpx) solid solutions. We emphasize that the technique is designed for use in the basalt range of compositions. Stoichiometric and mass balance considerations lead to a system of five equations with four unknowns (amounts of ol, pl, and cpx, and the anorthite mole fraction in plagioclase (An#)) which is over-determined. One advantage of this new method is that a variety of petrologically "reasonable" constraints can be added and used to constrain the optimal solution. In applying this method, we have used both Microsoft Excel version 5.0 (with "Solver"), and Matlab. We provide here an Excel spreadsheet "PB"that can be used to carry out the calculations outlined below.

2. Methods

Suppose olivine (ol), plagioclase (pl), and pyroxene (cpx) are the phases crystallized from a parental liquid with a moles pl, b moles cpx , and c moles ol. The parental liquid is normalized to 1 mole cations (in the form of SiO₂, TiO₂, AlO_1.5, FeOT (total iron as FeO), MnO, MgO, CaO, NaO_0.5, KO_0.5, and PO_2.5). Suppose the pl composition is Na_x Ca_1-xAl_2-xSi_2+xO₈(x to be determined) and the cpx composition is Na_n(Ca,Mg,Fe)_mAl_hSi_sTi_tO₆(amounts of the minor components Nan, Alh, and Tit are assumed and the major components (Ca,Mg,Fe)_mand Si_swill be determined). The total moles of cations in the crystallized solid would be moles because there are five moles of cations in one mole pl, four moles of cations in one mole cpx, and three moles of cations in one mole ol. Using the notation of []₀and [] to represent the mole concentrations of an element in the parent and daughter liquids respectively, we have 5 equations of the general form: daughter concentration = [total mass in parent liquid - mass in crystallized solid] / liquid left:

Eq. 1
Eq. 2
Eq. 3
Eq. 4
Eq. 5

There are five equations and four major unknowns (a, b, c, x), so five sets of solutions can be obtained by sequentially excluding one of the five equations. Instead of all 5, we only used 3 sets of solutions: without Na (excluding Eq. 1), without Si (excluding Eq. 3), and without Ti (excluding Eq. 5). These three are used because Na and Si have the largest analytical errors and Ti is special because it is used as an incompatible element in the calculation. Ti can be replaced by other incompatible element that can be analyzed with high precision (spreadsheet "PB"provides this option).

The above equations (1-5) are solved by substitution, and we have tried almost all possible orders of substitution for each set of the 3 sets of solutions. Certain orders of substitution produce equations with a squared term or close to zero denominators. The solutions we used in the spreadsheet "PB"do not have these problems. Because the final expressions are very long, we used a lot of symbols (like m1, h1 in the spreadsheet) to represent combinations of known parameters during substitutions to facilitate the algebra. The use of these symbols also helped to detect near-zero denominators. Because the steps used to solve the simultaneous equations are so numerous and the intermediate equations are so cumbersome, we do not include their full derivation. To satisfy potential users of the method that the solutions are in fact correct, we provide a "test" sheet for formulating and solving simple synthetic mass balance problems. By using the test sheet, readers can easily satisfy themselves that the calculations in the spreadsheet "PB" are valid.

3. Assumptions, Petrological Constraint and Optimization Strategies

The above phase calculations assume that all the daughter liquids are derived from a homogeneous parental liquid by subtraction of crystallization solids only and that the possible phases crystallized are ol, pl and cpx. Therefore only subalkalic basalts with more than ~5 weight percent MgO can be used in the calculation.

In general, we use idealized phase compositions and ignore minor mineralogic constituents that occur at levels of < 400 ppm in common crystallizing phases. Thus, we use ideal ol ([MgFe]₂SiO₄) and pl (Na_xCa_1-xAl_2-xSi_2+xO₈) composition. We use cpx composition Na_0.012 (Ca,Mg,Fe)_1.922 Al_0.12Si_1.934Ti_0.012O₆in the calculation before optimization because we have been interested in applying this technique to MORB and these values of the minor elements Na, Al, and Ti are appropriate for MORB compositions (e.g., Niu and Batiza, 1994). During the optimization carried out by "PB", the minor components Al, Na, Ti in cpx are optimized within given limits (from Na_0.004 (Ca,Mg,Fe)_1.974Al_0.04Si_1.978Ti_0.004O₆to Na_0.02 (Ca,Mg,Fe)_1.87Al_0.2Si_1.89Ti_0.02O₆). These limits, again, are appropriate for MORB lavas. Those users wishing to use the program for other basalt suites, may change these starting values and ranges for minor elements in cpx in the appropriate places in "PB".

In accordance with the assumption of a single parental liquid, we find that this method works best when the array of liquids to be modeled have chemistry that varies smoothly and monotonously. For natural lava suites showing scatter about such a trend, we first remove outliers and fit the trend to a 3d order polynomial. We then choose 6 liquid compositions with decreasing MgO from a well-defined part of LLD with a large number of samples and little scatter as input parent and daughter liquids for the spreadsheet.

In "PB", the parent-daughter calculations progress in a pair-wise progression to lower MgO, and the optimization considers all five pairs simultaneously in the best-fit calculation. We only use 6 liquid compositions as a compromise because addition of more liquid compositions decrease the speed of optimization significantly.

Optimization is needed for a number of reasons. First, the difference among the 3 sets of solutions (without Si, Na, and Ti) may reflect a violation of the assumption of a single parental liquid. PB minimizes the difference among the three independent solutions by optimizing the liquid compositions used. Secondly, we find that some solutions are petrologically unreasonable, for example showing an increase in calculated An# with decreasing temperature. By optimizing the liquid compositions used, we also require that the calculated An# of pl should decrease as crystallization proceeds and that the magnitude of this decrease (slope on a plot of An# vs. MgO) should be close (factor of 0.5 to 2, changeable in the spreadsheet "PB") to that from experiments. We used Panjasawatwong et al.’s (1995) empirical equation to calculate the An# of pl in equilibrium with a liquid to obtain the slope of the experimentally determined An# decrease against MgO. In their equation, we used a pressure of 1 kbar and temperature approximated by T(K) = 0.1936*MgO^3 - 3.6529*MgO^2 + 39.351*MgO + 1294.1 empirically calibrated with multiple runs of liquid LLDs using Weaver and Langmuir (1990). Users wishing to do so may change these values in "PB".

Using PB, it is also possible to calculate Ca content in cpx (see Appendix). We found that the calculated number of Ca atoms in one cpx formula (Ca# of cpx) was often over 1 (which of course is impossible). It turns out that the calculated Ca# of cpx is very sensitive to small changes in the liquid compositions (even smaller than analytical error) and to the assumed contents of minor elements in cpx. To remedy the problem of unreasonable Ca contents of cpx, and also to get around the problem of sensitivity of the Ca# calculation, we constrain Ca# to be below 0.85 only. Users wishing more stringent constraints can change PB accordingly. In addition to the An# in plagioclase, and the Ca# of cpx constraints, we impose a condition of no negative mineral amounts. In general, the goal of the optimization is to minimize the discrepancy among the 3 sets of solutions with minimum changes in the liquid compositions (from those initially input) and minimum changes in the solid phase proportions as crystallization proceeds. We impose the condition of minimum changes in the solid phase proportions to prevent petrologically unreasonable zigzag changes, and lack of convergence on a stable solution.

Results of PB do not always converge, and instabilities have in some cases been found. The set of optimization constraints given above seems to minimize unstable behavior and increases the chances of convergence, at least for the suites of MORB and oceanic island basalts we have used in experiments. If the solution fails to converge, the most likely reason is a poorly defined LLD (too few samples and/or large scatter). To aid convergence, the user may also relax the constraints for Ca# of cpx and An# of pl. Since this method is an iterative rather than an inverse method, we do not explore the solution space. The calculation procedure and detailed instructions on how to use it are in the spreadsheet "PB".

Appendix

Suppose the composition of cpx is Na_n Ca_xx (Mg,Fe)_m-xx Al_h Si_s Ti_t O₆ (xx to be determined). We have:

Eq. 6
Eq. 7
Eq. 8
Eq. 9
Eq. 10
Eq. 11

There are six equations and five unknowns (a, b, c, x, xx). So the 5 unknowns can be determined by any 5 of the 6 equations. To avoid having to choose between positive and negative square roots in the solution of an equation with a squared term, we devise the following solution.

Suppose in a Ca (abscissa) vs. Na diagram, the straight line is:

Eq. 12

In a Ti (abscissa) vs. Ca diagram, the straight line is:

Eq. 13

In an Al (abscissa) vs. Ca diagram, the straight line is:

Eq. 14

In a Mg+Fe (abscissa) vs. Ca diagram, the straight line is:

Eq.15

In an Si (abscissa) vs. Ca diagram, the straight line is:

Eq. 16

Plug equations 6-11 in the above 5 equations. Because we only use slopes and intercepts in this calculation, only the ratios of plagioclase, pyroxene and olivine are defined (1:(b/a):(c/a), essentially two unknowns) rather than their absolute amount (a, b, c, three unknowns). By reducing one unknown, we can avoid equations with a squared term.

References

Ghiorso, M. S. and R. O., Sack, 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, Contrib. Mineral. Petrol., 119, 197-212, 1995.

Niu, Y. and Batiza, R., Magmatic processes at a slow spreading ridge segment: 26 S Mid-Atlantic Ridge, J. Geophs. Res., 99: 19719-19740, 1994.

Panjasawatwong, Y., Danyushevsky, L.V., Crawford A.J., and Harris, K.L., An experimental study of the effects of melt composition on plagioclase- melt equilibria at 5 and 10 kbar: implications for the origin of magmatic high-An plagioclase, Contrib Mineral Petrol, 118: 420-432, 1995.

Weaver, J.S. and Langmuir, C.H, Calculation of phase equilibrium in mineral-melt systems, Computers & Geosciences, 16, 1-19, 1990.