Pyroxene Geothermometer¶
Using opx and cpx solution models from Sack and Ghiorso (1994a, b, c):¶
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.¶
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.¶
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