ENKI

Fe-Ti oxide geothermobarometer

Constructed based on Fe-Ti exchange for oxide pairs, using the methodology of Ghiorso and Evans (2008)

Required Python code to load the phase library …

import numpy as np
from scipy import optimize as optim
import thermoengine as thermo

Get access to a thermodynamic database (by default, the Berman (1988)/MELTS database).

modelDB = thermo.model.Database()

Create a Python reference to the Spinel (“Mag”) and Rhombehedral oxide (“Ilm”) solution phase class.

Mag = modelDB.get_phase('SplS')
Ilm = modelDB.get_phase('Rhom')

Optional - Obtain some properties of the selected Oxide solution

Name, formulas of endmembers, molecular weights of endmembers, abbreviation, number of endmember components, names of endmembers

print (Mag.props['phase_name'])
print ('Num of endmembers: ', Mag.props['endmember_num'])
print (Mag.props['formula'])
print (Mag.props['endmember_name'])
print ()
print (Ilm.props['phase_name'])
print ('Num of endmembers: ', Ilm.props['endmember_num'])
print (Ilm.props['formula'])
print (Ilm.props['endmember_name'])
Spinel
Num of endmembers:  5
['FeCr2O4' 'FeAl2O4' 'Fe3O4' 'MgAl2O4' 'Fe2TiO4']
['chromite' 'hercynite' 'magnetite' 'spinel' 'ulvospinel']

Ilmenite ss
Num of endmembers:  5
['MgTiO3' 'Fe2O3' 'FeTiO3' 'MnTiO3' 'Al2O3']
['geikielite' 'hematite' 'ilmenite' 'pyrophanite' 'corundum']

Step 1 - Input compositions of oxide pairs

case = 1
  • case 0: Brad’s composition

  • case 1: Oxide pair from website geothermometer

if case == 0:
    Mag_mol_oxides = thermo.chem.format_mol_oxide_comp(
        {'SiO2':0.0, 'TiO2':4.15, 'Al2O3':2.83, 'Fe2O3':57.89, 'FeO':32.48,
        'MnO':0.35, 'MgO':1.0, 'CaO':0.0, 'Na2O':0.0, 'Cr2O3':0.07},convert_grams_to_moles=True)
    Ilm_mol_oxides = thermo.chem.format_mol_oxide_comp(
        {'SiO2':0, 'TiO2':37.27, 'Al2O3':0.29, 'Fe2O3':29.62,'FeO':29.19,
        'MnO':0.39, 'MgO':2.0, 'CaO':0.0, 'Na2O':0.0, 'Cr2O3':0.01},convert_grams_to_moles=True)
else: # example from website, gives 724 °C, 1.61 delta NNO 0.894 a TiO2
    Mag_mol_oxides = thermo.chem.format_mol_oxide_comp(
        {'SiO2':0.0, 'TiO2':4.35, 'Al2O3':1.94, 'Fe2O3':0.00, 'Cr2O3':0.18, 'FeO':86.34,
        'MnO':0.44, 'MgO':1.2, 'CaO':0.0, 'Na2O':0.0},convert_grams_to_moles=True)
    Ilm_mol_oxides = thermo.chem.format_mol_oxide_comp(
        {'SiO2':0.0, 'TiO2':28.73, 'Al2O3':0.35, 'Fe2O3':0.00, 'Cr2O3':0.0, 'FeO':65.98,
        'MnO':0.23, 'MgO':1.02, 'CaO':0.0, 'Na2O':0.0},convert_grams_to_moles=True)

Step 2 - Convert analytical composition to moles of endmember components

Note that a method - test_endmember_comp() - is called to test the validity of the projected composition
Also note that the “intrinsic” conversion routines now take care of FeO-Fe2O3 conversion to balance the cation-anion ratio of the phase. Input compositions of Fe2O3 and/or FeO are adjusted to balance charge and maintain phase stoichiometry.
def validate_endmember_comp(moles_end, phase):
    print(phase.props['phase_name'])
    sum = 0.0
    for i in range(0,phase.props['endmember_num']):
        print ("mole number of {0:10.10s} = {1:13.8f}".format(
            phase.props['endmember_name'][i], moles_end[i]))
        sum += moles_end[i]
    if not phase.test_endmember_comp(moles_end):
        print ("Calculated composition is infeasible!")

Mag_moles_end = Mag.calc_endmember_comp(
    mol_oxide_comp=Mag_mol_oxides, method='intrinsic', normalize=True)
validate_endmember_comp(Mag_moles_end, Mag)
print()
Ilm_moles_end = Ilm.calc_endmember_comp(
    mol_oxide_comp=Ilm_mol_oxides, method='intrinsic', normalize=True)
validate_endmember_comp(Ilm_moles_end, Ilm)
Spinel
mole number of chromite   =    0.00266617
mole number of hercynite  =   -0.02419364
mole number of magnetite  =    0.83193038
mole number of spinel     =    0.06702845
mole number of ulvospinel =    0.12256865

Ilmenite ss
mole number of geikielite =    0.03853892
mole number of hematite   =    0.44719307
mole number of ilmenite   =    0.50410315
mole number of pyrophanit =    0.00493747
mole number of corundum   =    0.00522739

Implement a Fe-Ti oxide geothermometer.

Consider Fe-Ti exchange between oxides

Rhom(Ilm) + Spinel(Mag) = Spinel (Ulv) + Rhom(Hm)

FeTiO3 + Fe3O4 = Fe2TiO4 + Fe2O3

Fix the pressure at 2000 bars

P = 2000.0

Corection terms from Ghiorso and Evans (2008) that modify the MELTS models

Correction terms for ulvospinel derived in Ghiorso and Evans (2008)

def UlvCorr(t, correctReaction=True):
    tr = 298.15
    h = - 162.0 + 284.5
    s = 0.0
    if correctReaction:
        h += 2039.106175
        s +=    1.247790
    l1 = - 0.039452*np.sqrt(4.184)
    l2 = 7.54197e-5*np.sqrt(4.184)
    h = h + 0.5*l1*l1*(t*t-tr*tr) + (2.0/3.0)*l1*l2*(t*t*t - tr*tr*tr) + 0.25*l2*l2*(t*t*t*t - tr*tr*tr*tr)
    s = s + l1*l1*(t - tr) + l1*l2*(t*t - tr*tr) + (1.0/3.0)*l2*l2*(t*t*t - tr*tr*tr)
    return h - t*s

Ghiorso and Evans (2008) used the Vinet integral; MELTS uses the Berman integral

We must substract the latter from computed chemical potentials and add in the former.

def BermanVint(t, p, v0, v1, v2, v3, v4):
    pr = 1.0
    tr = 298.15
    return v0*((v1/2.0-v2)*(p*p-pr*pr)+v2*(p*p*p-pr*pr*pr)/3.0 + (1.0-v1+v2+v3*(t-tr)+v4*(t-tr)*(t-tr))*(p-pr))
def VinetVint(t, p, v0, alpha, K, Kp):
    eta = 3.0*(Kp-1.0)/2.0
    x   = 1.0
    x0  = 1.0
    pr  = 1.0
    tr  = 298.15

    iter = 0
    while True:
        fn = x*x*(p/10000.0) - 3.0*K*(1.0-x)*np.exp(eta*(1.0-x)) - x*x*alpha*K*(t-tr)
        dfn = 2.0*x*(p/10000.0) + 3.0*K*(1.0+eta*(1.0-x))*np.exp(eta*(1.0-x)) - 2.0*alpha*K*(t-tr)
        x = x - fn/dfn
        iter += 1
        if ((iter > 500) or (fn*fn < 1.0e-15)):
            break
    # print (iter, x)

    iter = 0
    while True:
        fn = x0*x0*(pr/10000.0) - 3.0*K*(1.0-x0)*np.exp(eta*(1.0-x0)) - x0*x0*alpha*K*(t-tr)
        dfn = 2.0*x0*(pr/10000.0) + 3.0*K*(1.0+eta*(1.0-x0))*np.exp(eta*(1.0-x0)) - 2.0*alpha*K*(t-tr)
        x0 = x0 - fn/dfn
        iter += 1
        if ((iter > 500) or (fn*fn < 1.0e-15)):
            break
    # print (iter, x0)

    a  = (9.0*v0*K/(eta*eta))*(1.0 - eta*(1.0-x))*np.exp(eta*(1.0-x))
    a += v0*(t-tr)*K*alpha*(x*x*x - 1.0) - 9.0*v0*K/(eta*eta)
    a -= (9.0*v0*K/(eta*eta))*(1.0 - eta*(1.0-x0))*np.exp(eta*(1.0-x0))
    a -= v0*(t-tr)*K*alpha*(x0*x0*x0 - 1.0) - 9.0*v0*K/(eta*eta)

    return -a*10000.0 + p*v0*x*x*x - pr*v0*x0*x0*x0

Berman integral for the reaction FeTiO3 + Fe3O4 = Fe2TiO4 + Fe2O3

def rBerVint(T, P):
    vIntBerMag = BermanVint(T, P, 4.452, -0.582E-6, 1.751E-12, 30.291E-6, 138.470E-10)
    vIntBerUlv = BermanVint(T, P, 4.682, 0.0, 0.0, 0.0, 0.0)
    vIntBerHem = BermanVint(T, P, 3.027, -0.479e-6, 0.304e-12, 38.310e-6, 1.650e-10)
    vIntBerIlm = BermanVint(T, P, 3.170, -0.584e-6, 1.230e-12, 27.248e-6, 29.968e-10)
    return vIntBerUlv + vIntBerHem - vIntBerMag -vIntBerIlm

Vinet integral for the reaction FeTiO3 + Fe3O4 = Fe2TiO4 + Fe2O3

def rVinetVint(T, P):
    vIntVinetMag = VinetVint(T, P, 4.452, 30.291E-6, 171.821, 9.3387)
    vIntVinetUlv = VinetVint(T, P, 4.682, 30.291E-6, 171.821, 9.3387)
    vIntVinetHem = VinetVint(T, P, 3.027, 38.310E-6, 208.768, 1.64992)
    vIntVinetIlm = VinetVint(T, P, 3.170, 27.248E-6, 171.233, 6.21289)
    return vIntVinetUlv + vIntVinetHem - vIntVinetMag - vIntVinetIlm

This method computes the free energy of the exchange reaction …

def deltaG(T, P, mag_mols, ilm_mols):
    muMag = Mag.chem_potential(T, P, mol=mag_mols, endmember=2)
    muUlv = Mag.chem_potential(T, P, mol=mag_mols, endmember=4) + UlvCorr(T)
    muIlm = Ilm.chem_potential(T, P, mol=ilm_mols, endmember=2)
    muHem = Ilm.chem_potential(T, P, mol=ilm_mols, endmember=1)
    deltaG = muUlv + muHem - muIlm - muMag - rBerVint(T, P) + rVinetVint(T, P)
    return deltaG

This next function is used by the minimizer to zero the free energy of the exchange reaction …

def boundary(P, Tlims, deltaG, mag_mols, ilm_mols):
    Afun = lambda T, P=P: deltaG(T, P, mag_mols, ilm_mols)
    Tbound = optim.brentq(Afun, Tlims[0], Tlims[1])
    return Tbound

Calculate the equilibrium temperature for this oxide pair …

Teq = boundary(P, [500.,2000.], deltaG, Mag_moles_end, Ilm_moles_end)
print('Equilibrium Temp = ', Teq-273.15, ' °C')
Equilibrium Temp =  724.072806179801  °C

Calculate Oxygen fugacity from the reaction

O2 + 4 Fe3O4 = 6 Fe2O3

Note that the properties of oxygen are defined here for consistency instead of using the built-in functions.
Also note that the chemical potentials of hematite and magnetite are adjusted to remove the Berman-type volume integrals and replace them with the Vinet-type volume integrals to be consistent with Ghiorso and Evans (2008)
def muO2(t, p):
    tr = 298.15
    hs = 23.10248*(t-tr) + 2.0*804.8876*(np.sqrt(t)-np.sqrt(tr)) - 1762835.0*(1.0/t-1.0/tr) \
       - 18172.91960*np.log(t/tr) + 0.5*0.002676*(t*t-tr*tr)
    ss = 205.15 + 23.10248*np.log(t/tr)  - 2.0*804.8876*(1.0/np.sqrt(t)-1.0/np.sqrt(tr)) \
       - 0.5*1762835.0*(1.0/(t*t)-1.0/(tr*tr)) + 18172.91960*(1.0/t-1.0/tr) + 0.002676*(t-tr)
    return hs - t*ss
def deltaNNO (T, P, mag_mols, ilm_mols):
    muHem  = Ilm.chem_potential(T, P, mol=ilm_mols, endmember=1)
    muHem -= BermanVint(T, P, 3.027, -0.479e-6, 0.304e-12, 38.310e-6, 1.650e-10)
    muHem += VinetVint(T, P, 3.027, 38.310E-6, 208.768, 1.64992)
    muMag =  Mag.chem_potential(T, P, mol=mag_mols, endmember=2)
    muMag -= BermanVint(T, P, 4.452, -0.582E-6, 1.751E-12, 30.291E-6, 138.470E-10)
    muMag += VinetVint(T, P, 4.452, 30.291E-6, 171.821, 9.3387)
    muOxy  = muO2(T, P)
    logfO2 = (6.0*muHem - 4.0*muMag -  muOxy)/(8.3144598*T)/np.log(10.0)
    return logfO2 - (-25018.7/T + 12.981 + 0.046*(P-1.0)/Teq -0.5117*np.log(T))

Calculate the equilibrium oxygen fugacity for this oxide pair …

print(deltaNNO(Teq, P, Mag_moles_end, Ilm_moles_end))
1.613240788544294

Calculate the temperature for Fe-Mg exchange

FeAl2O4 (hercynite) + MgTiO3 (geikielite) = MgAl2O4 (spinel) + FeTiO3 (ilmenite)

The method below is used by the minimizer to evaluate the free energy change of the Fe-Mg exchange reaction …

def deltaGfemg(T, P, mag_mols, ilm_mols):
    muSpn  = Mag.chem_potential(T, P, mol=mag_mols, endmember=3)
    muSpn -= BermanVint(T, P, 3.977, -0.489E-6, 0.0, 21.691E-6, 50.528E-10)
    muSpn += VinetVint(T, P, 3.977, 21.691E-6, 204.499, 4.0)

    muHer  = Mag.chem_potential(T, P, mol=mag_mols, endmember=1)
    muHer -= BermanVint(T, P, 0.973948*4.184, 0.0, 0.0, 0.0, 0.0)
    muHer += VinetVint(T, P, 0.973948*4.184, 21.691E-6, 204.499, 4.0)

    muIlm  = Ilm.chem_potential(T, P, mol=ilm_mols, endmember=2)
    muIlm -= BermanVint(T, P, 3.170, -0.584e-6, 1.230e-12, 27.248e-6, 29.968e-10)
    muIlm += VinetVint(T, P, 3.170, 27.248E-6, 171.233, 6.21289)

    muGei  = Ilm.chem_potential(T, P, mol=ilm_mols, endmember=0)
    muGei -= BermanVint(T, P, 3.086, -0.584e-6, 1.230e-12, 27.248e-6, 29.968e-10)
    muGei += VinetVint(T, P, 3.086, 27.2476341e-6, 171.240, 6.21527)

    deltaG = muSpn + muIlm - muHer - muGei
    return deltaG

Compute the Fe-Mg exchange temperature (if possible) …

Tlow  = deltaGfemg(500.0, P, Mag_moles_end, Ilm_moles_end)
Thigh = deltaGfemg(2000.0, P, Mag_moles_end, Ilm_moles_end)
if np.sign(Tlow) != np.sign(Thigh):
    Tfemg = boundary(P, [500.,2000.], deltaGfemg, Mag_moles_end, Ilm_moles_end)
    print('Fe-Mg Equilibrium Temp = ', Tfemg-273.15, ' °C')
else:
    print('No Fe-Mg equilibration temperature found.')
No Fe-Mg equilibration temperature found.

Calculate the activity of TiO2 relative to rutile saturation

Rut = modelDB.get_phase('Rt')
def aTiO2(T, P, mag_mols, ilm_mols):
    muUlv  = Mag.chem_potential(T, P, mol=mag_mols, endmember=4) + UlvCorr(T, correctReaction=False)
    muUlv -= BermanVint(T, P, 4.682, 0.0, 0.0, 0.0, 0.0)
    muUlv += VinetVint(T, P, 4.682, 30.291E-6, 171.821, 9.3387)
    muIlm  = Ilm.chem_potential(T, P, mol=ilm_mols, endmember=2)
    muIlm -= BermanVint(T, P, 3.170, -0.584e-6, 1.230e-12, 27.248e-6, 29.968e-10)
    muIlm += VinetVint(T, P, 3.170, 27.248E-6, 171.233, 6.21289)
    muRut = Rut.chem_potential(T, P)
    return np.exp(-(muRut+muUlv-2.0*muIlm)/(8.3143*T))
aTiO2(Teq, P, Mag_moles_end, Ilm_moles_end)
0.8938677107462636