MagmaForge: Mantle Melting Pseudosection

Demonstrates how to run a fractional (isentropic) decompression routine using MagmaForge.

Open this code in an executable MyBinder instance (MyBinder links may be slow to load– please be patient!):

Initialization

Import necessary packages:

from thermoengine import magmaforge

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

Define a bulk composition

DMM_comp_adj = {
    'SiO2': 44.7354,
    'TiO2': 0.130074,
    'Al2O3': 3.98226,
    'Fe2O3': 0.275456,
    'Cr2O3': 0.570323,
    'FeO': 8.01254,
    'MgO': 38.752,
    'CaO': 3.1718,
    'Na2O': 0.130074
}

Run crystallization at different P slices

state_histories = []
P_vals  = np.linspace(2,0.5,10)
for P_GPa in P_vals:
    # Start at liquidus
    sys = magmaforge.System(comp=DMM_comp_adj,
                        T_liquidus=True,
                        P_GPa=P_GPa,
                        database='MELTS_pMELTS',
                        )
    # Crystallize to solidus
    sys.crystallize(method='equil',
                    T_step=5,
                    fix_fO2=False,
                    calc_args={'debug':0},
                    #melt_frac_cutoff=0.001
                    )

    state_histories.append(sys.history)

/workspaces/ThermoEngineLite/thermoengine/thermoengine/redox.py:109: RuntimeWarning: invalid value encountered in divide
  ln_Fe_oxide_ratio = np.squeeze(np.log(mol_oxides[:, OXIDES=='Fe2O3']/mol_oxides[:, OXIDES=='FeO']))
/workspaces/ThermoEngineLite/thermoengine/thermoengine/redox.py:116: RuntimeWarning: invalid value encountered in divide
  mol_oxides /= mol_oxide_tot[:,np.newaxis]
/workspaces/ThermoEngineLite/thermoengine/thermoengine/redox.py:109: RuntimeWarning: invalid value encountered in divide
  ln_Fe_oxide_ratio = np.squeeze(np.log(mol_oxides[:, OXIDES=='Fe2O3']/mol_oxides[:, OXIDES=='FeO']))
/workspaces/ThermoEngineLite/thermoengine/thermoengine/redox.py:116: RuntimeWarning: invalid value encountered in divide
  mol_oxides /= mol_oxide_tot[:,np.newaxis]

Calculate Liquidus, Solidus, and Melt Isopleths

T_liquidus = np.zeros_like(P_vals)
T_solidus = np.zeros_like(P_vals)
F_vals = np.linspace(0.10,0.90,9)
T_F_array = np.zeros([len(F_vals), len(P_vals)])
for P_idx,state_history in enumerate(state_histories):
    T_liquidus[P_idx] = state_history.get_temperatures('C')[0]
    T_solidus[P_idx] = state_history.get_temperatures('C')[-1]

    F = state_history.get_melt_fractions()
    for F_idx, F_val in enumerate(F_vals):
        T_idx = np.abs(F - F_val).argmin()
        T_F_array[F_idx,P_idx] = state_history.get_temperatures('C')[T_idx]

Create Plot

fig, ax = plt.subplots(nrows=1, figsize=(4,4))

ax.plot(T_liquidus, P_vals, lw=2, color='r')
ax.plot(T_solidus, P_vals, lw=2, color='k')
ax.plot(T_F_array.T, P_vals, lw=1, ls=':', color='k')

ax.fill_betweenx(P_vals,1100,T_solidus,color='black', alpha=0.2)
ax.fill_betweenx(P_vals,1800,T_liquidus,color='red', alpha=0.2)

ax.set_xlabel('T (°C)')
ax.set_xlim([1450-273.15,2000-273.15])

ax.set_ylabel('P (GPa)')
ax.set_ylim([P_vals.max(), P_vals.min()])

ax.text(1550-273.15,1.5,'solid')
ax.text(1900-273.15,1.0,'liquid', color='r')
ax.text(1290,0.48,'melt fraction isopleths', color='k', fontsize=8)

ax.set_title('Partial Melting of DMM', pad=20)

plt.tight_layout()
plt.show()