MagmaForge: Find pMELTS Liquidus with buffered fO2

At present the standard liquidus finding routine does not work with pMELTS when you want to buffer the oxygen fugacity. We plan to fix this limitation soon, but for now this notebook provides a work around (though will take a little longer to calculate).

Initialization

Import necessary packages:

from thermoengine import magmaforge, redox, rockychem, model

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import root_scalar

Define a bulk composition:

# MORB
wt_oxides = pd.Series({
        'SiO2':  48.68,
        'TiO2':   1.01,
        'Al2O3': 17.64,
        'Fe2O3':  0.0,
        'Cr2O3':  0.0425,
        'FeO':    8.0,
        'MgO':    9.10,
        'CaO':   12.45,
        'Na2O':   2.65,
        'K2O':    0.03,
        'P2O5':   0.08,
        'H2O':    0.2
        })

Define some utility functions for solving for the fO2 buffer

def misfit(x, P, wt_oxides, logfO2):
    T = x
    logfO2_abs = redox.redox_buffer(T, P, buffer=logfO2[0]) + logfO2[1]
    mol_oxides = rockychem.convert_bulk_comp(wt_oxides=wt_oxides, to='mol_oxides', drop_missing=False)
    adj_mol_oxides = redox.adjust_liquid_redox_ratio(logfO2=logfO2_abs, T=T, P=P, mol_oxides=mol_oxides)
    adj_wt_oxides = rockychem.convert_bulk_comp(mol_oxides=adj_mol_oxides, to='wt_oxides', drop_missing=False)

    sys = magmaforge.System(comp=adj_wt_oxides,
                            T_liquidus=True,
                            P_bar=P,
                            # logfO2=logfO2,
                            database="MELTS_pMELTS")

    return sys.T_K - T

def find_liquidus(P_bar, logfO2):
    P = P_bar
    res = root_scalar(misfit, args=(P, wt_oxides, logfO2), bracket=(1000.0, 3000.0))
    assert res.converged, "Something went wrong, check your inputs!"

    return res.root - 273.15

Find the liquidus for a single pressure

Set the pressure and fO2 buffer, then run the calculation and print the output:

Tliq = find_liquidus(P_bar = 10000.0,
                     logfO2=('IW', 0.0)
                     )

print(f"The Liquidus Temperature is {Tliq :.1f}˚C")

The Liquidus Temperature is 1307.3˚C

Plot the liquidus over a range of Delta logfO2

Set the minimum and maximum Delta logfO2 values, and the buffer to use, then run the calculations:

P_bar = 10000.0
buffer = 'IW'
D_fO2_min = -2
D_fO2_max = 2
n = 5

liquidus = np.zeros(n)
logfO2s = np.linspace(-2, 2, n)

for i in range(n):
    liquidus[i] = find_liquidus(P_bar, (buffer, logfO2s[i]))

fig, ax = plt.subplots()

ax.plot(logfO2s, liquidus, marker='s')
ax.set_xlabel(f'$\Delta \log(fO_2)$ ({buffer})')
ax.set_ylabel('Liquidus Temperature (˚C)')

plt.show()

/workspaces/ThermoEngineLite/doc/examples/_2_magmaforge/plot_magmaforge_find_liquidus_buffered_pMELTS.py:109: DeprecationWarning: invalid escape sequence '\D'
  ax.set_xlabel(f'$\Delta \log(fO_2)$ ({buffer})')

Plot the liquidus over a pressure range

Set the minimum and maximum pressure and number of points and run calculations:

Pmin_bar = 10000.0
Pmax_bar = 30000.0
n = 5

liquidus = np.zeros(n)
pressures = np.linspace(Pmin_bar, Pmax_bar, n)
logfO2 = ('IW', 0.0)

for i in range(n):
    liquidus[i] = find_liquidus(pressures[i], logfO2)

fig, ax = plt.subplots()

ax.plot(liquidus, pressures, marker='o')
ax.invert_yaxis()
ax.set_ylabel('Pressure (bar)')
ax.set_xlabel('Temperature (˚C)')

plt.show()