ENKI

Pyroxene Geothermometer

Using opx and cpx solution models from Sack and Ghiorso (1994a, b, c):

Sack RO, Ghiorso MS (1994a) Thermodynamics of multicomponent pyroxenes: I. Formulation of a general model. Contrib Mineral Petrol 116, 277-286
Sack RO, Ghiorso MS (1994b) Thermodynamics of multicomponent pyroxenes: II. Phase relations in the quadrilateral. Contrib Mineral Petrol 116, 287-300
Sack RO, Ghiorso MS (1994c) Thermodynamics of multicomponent pyroxenes: III. Calibration of Fe2+(Mg)-1, TiAl(MgSi)-1, TiFe3+(MgSi)-1, AlFe3+(MgSi)-1, NaAl(CaMg)-1, Al2(MgSi)-1 and Ca(Mg)-1 exchange reactions between pyroxenes and silicate melts. Contrib Mineral Petrol 118, 271-296

Initialize some required packages, and load the phase library.

from ctypes import cdll
from ctypes import util
from rubicon.objc import ObjCClass, objc_method
cdll.LoadLibrary(util.find_library('phaseobjc'))
/Users/ghiorso/anaconda3/lib/python3.7/site-packages/rubicon/objc/ctypes_patch.py:24: UserWarning: rubicon.objc.ctypes_patch has only been tested with Python 3.4 through 3.6. The current version is sys.version_info(major=3, minor=7, micro=6, releaselevel='final', serial=0). Most likely things will work properly, but you may experience crashes if Python's internals have changed significantly.
  .format(sys.version_info)
<CDLL '/usr/local/lib/libphaseobjc.dylib', handle 7f9754c053b0 at 0x7f977852c090>

Define some conversion functions that …

… take dictionaries of oxide names and oxides values and return molecular weights and arrays of molar concentrations.

def oxide_mw (formulas=["H2O"]):
    result = {}
    PhaseBase = ObjCClass('PhaseBase')
    for formula in formulas:
        obj = PhaseBase.alloc().init()
        obj.setPhaseFormula_(formula)
        result[formula]= obj.mw
    return result

import ctypes
def oxides_wts_to_element_moles (oxides={"H2O" : 100.0}):
    e = (ctypes.c_double*107)()
    ctypes.cast(e, ctypes.POINTER(ctypes.c_double))
    for i in range (0, 107):
        e[i] = 0.0
    PhaseBase = ObjCClass('PhaseBase')
    for formula, value in oxides.items():
        obj = PhaseBase.alloc().init()
        obj.setPhaseFormula_(formula)
        moles = value/obj.mw
        elements = obj.formulaAsElementArray
        for i in range (0, 107):
            coeff = elements.valueAtIndex_(i)
            if coeff != 0.0:
                e[i] += coeff*moles
    return e

def element_moles_to_pyx_moles(e):
    m = (ctypes.c_double*nc)()
    ctypes.cast(m, ctypes.POINTER(ctypes.c_double))
    p = 2000.0
    Na = 11
    Mg = 12
    Al = 13
    Si = 14
    Ca = 20
    Ti = 22
    Cr = 24
    Mn = 25
    Fe = 26
    sumcat  = e[Na] +     e[Mg] +     e[Al] +     e[Si] +     e[Ca] +     e[Ti] +     e[Cr] +     e[Mn] + e[Fe]
    sumchg  = e[Na] + 2.0*e[Mg] + 3.0*e[Al] + 4.0*e[Si] + 2.0*e[Ca] + 4.0*e[Ti] + 3.0*e[Cr] + 2.0*e[Mn]
    if e[Na]+e[Ca] > 0.25*sumcat:
        corrSi = 4.0*(e[Na]+e[Ca]) - sumcat
    else:
        corrSi = 0.0
    sumcat += corrSi;

    # catch low-P oxidized samples and acmites
    if (p < 1000.0) or (e[Na] > e[Al]):
        fe3 = 3.0*sumcat - sumchg - 2.0*e[Fe]
        fe2 = e[Fe] - fe3
        if fe3 < 0.01*e[Fe]:
            fe3 = 0.01*e[Fe]
            fe2 = 0.99*e[Fe]
        if fe2 < 0.01*e[Fe]:
            fe2 = 0.01*e[Fe]
            fe3 = 0.99*e[Fe]
    else:
        fe2 = e[Fe]
        fe3 = 0.0

    m[0] = -fe3/2.0 - fe2 - e[Mn] - e[Al]/2.0 - e[Cr]/2.0 + e[Ca] + e[Na]/2.0 - e[Ti]
    m[1] =  fe3/4.0 + fe2/2.0 + e[Mn]/2.0 + e[Al]/4.0 + e[Cr]/4.0 - e[Ca]/2.0 + e[Mg]/2.0 - e[Na]/4.0
    m[2] =  fe2 + e[Mn]
    m[3] = -fe3/2.0 + e[Al]/2.0 + e[Cr]/2.0 - e[Na]/2.0 + e[Ti]
    m[4] =  fe3/2.0 - e[Al]/2.0 - e[Cr]/2.0 + e[Na]/2.0 + e[Ti]
    m[5] =  fe3/2.0 + e[Al]/2.0 + e[Cr]/2.0 - e[Na]/2.0 - e[Ti]
    m[6] =  e[Na]
    return m

Test the oxide formula to molecular weight method.

print (oxide_mw(["Al2O3", "SiO2"]))
{'Al2O3': 101.96127999999999, 'SiO2': 60.0843}

Implement a two-pyroxene geothermometer.

Phase

SiO2

TiO2

Al2O3

FeO

MnO

MgO

CaO

Na2O

cpx

5 1.946 15385

0 .7111 53846

0.150 769231

12 .80 846 154

0. 556 923 077

12 .69 576 923

20 .63 307 692

0 .3811 53846

:m ath:sigma

0 .3252 45469

0 .1598 08058

0.025 443754

0. 305 754 049

0. 038 860 698

0. 250 984 829

0. 200 434 912

0 .0158 30837

opx

5 0.929 25926

0 .4255 55556

0.128 888889

28 .49 518 519

1. 103 703 704

18 .33 037 037

0 .97 962 963

0 .0251 85185

:m ath:sigma

0 .4603 25353

0 .1076 43762

0.023 912233

0. 493 233 993

0. 045 161 943

0. 257 241 005

0. 022 951 702

0 .0070 00203

Values in wt%. Averages and standard deviations computed from analyzed pyroxenes found in the late eruptive units.

Instantiate a clinopyroxene with the specified composition.

As an illustration of use, compute and print properties at 800 °C and 200 MPa. Properties are output as a Python dictionary.
Output the number of components, their names, and their formulas.
CpxBerman = ObjCClass('CpxBerman')
cpx = CpxBerman.alloc().init()
nc = cpx.numberOfSolutionComponents()
e = oxides_wts_to_element_moles ({'SiO2':51.94615385, 'TiO2':0.711153846, 'Al2O3':0.150769231, 'FeO':12.80846154,
                                  'MnO':0.556923077, 'MgO':12.69576923, 'CaO':20.63307692, 'Na2O':0.381153846})
mCpx = element_moles_to_pyx_moles(e)

if (cpx.testPermissibleValuesOfComponents_(mCpx) == 1):
    print ('Cpx composition is feasible')
else:
    print ('Cpx composition is infeasible')

t = 1073.15 # K
p = 2000.0  # bars
potential = cpx.getChemicalPotentialFromMolesOfComponents_andT_andP_(mCpx, t, p)

for i in range (0, nc):
    component = cpx.componentAtIndex_(i)
    print("{0:>20s}{1:15.2f}".format(component.phaseName, potential.valueAtIndex_(i)))
Cpx composition is feasible
            diopside    -3466683.65
      clinoenstatite    -3352376.72
        hedenbergite    -3150901.60
   alumino-buffonite    -3575864.14
           buffonite    -3151406.39
            essenite    -3221797.55
             jadeite    -3336129.22

Instantiate an orthopyroxene with the specified composition.

As an illustration of use, compute and print properties at 800 °C and 200 MPa. Properties are output as a Python dictionary.
Output the number of components, their names, and their formulas.
OpxBerman = ObjCClass('OpxBerman')
opx = OpxBerman.alloc().init()
nc = opx.numberOfSolutionComponents()
e = oxides_wts_to_element_moles ({'SiO2':50.92925926, 'TiO2':0.425555556, 'Al2O3':0.128888889, 'FeO':28.49518519,
                                  'MnO':1.103703704, 'MgO':18.33037037, 'CaO':0.97962963, 'Na2O':0.025185185})
mOpx = element_moles_to_pyx_moles(e)

if (opx.testPermissibleValuesOfComponents_(mOpx) == 1):
    print ('Opx composition is feasible')
else:
    print ('Opx composition is infeasible')

t = 1073.15 # K
p = 2000.0  # bars
potential = opx.getChemicalPotentialFromMolesOfComponents_andT_andP_(mOpx, t, p)

for i in range (0, nc):
    component = opx.componentAtIndex_(i)
    print("{0:>20s}{1:15.2f}".format(component.phaseName, potential.valueAtIndex_(i)))
Opx composition is feasible
            diopside    -3470369.82
      clinoenstatite    -3356457.02
        hedenbergite    -3153990.55
   alumino-buffonite    -3556103.19
           buffonite    -3144518.00
            essenite    -3487657.71
             jadeite    -3582560.72

Define an Fe-Mg exchange reaction between opx and cpx.

CaMgSi2O6 [cpx] + CaFeSi2O6 [opx] = CaMgSi2O6 [opx] +CaFeSi2O6 [cpx]

Note that the get_properties function for the class instance returns a Python dictionary. The chemical potential of the endmember components are retrieved from this dictionary by using the name of the component as a key. Otherwise, the other thermodynamic properties are extensive (mass dependent) quantities and pertain to the phase as a whole.

def deltaG(t, p):
    cpxPotentials = cpx.getChemicalPotentialFromMolesOfComponents_andT_andP_(mCpx, t, p)
    opxPotentials = opx.getChemicalPotentialFromMolesOfComponents_andT_andP_(mOpx, t, p)
    return opxPotentials.valueAtIndex_(0) + cpxPotentials.valueAtIndex_(2) - cpxPotentials.valueAtIndex_(0) - opxPotentials.valueAtIndex_(2)

The Gibbs free energy computed by the deltaG function defined above must be zero at equilibrium. In order to find this zero, we … ## . . . import a minimizer routine from SciPy called BrentQ. We will use BrentQ to find the temperature that zeroes the Gibbs free energy of a reaction within a specified range of values.

from scipy.optimize import brentq

Solve for the temperature that zeroes the exchange free energy.

Upper and lower bounds on T are specified by Tmin and Tmax (both in K). The pressure is specified in bars.

Tmin = 500.0
Tmax = 1500.0
p = 2000.0
print ('Equilibrium T (°C) = ', brentq(deltaG, Tmin, Tmax, args=(p)) - 273.15)
Equilibrium T (°C) =  695.2941439766115