MagmaForge: Equilibrium Decompression

Demonstrates how to run an equilibrium (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  = {
        'SiO2': 44.71,
        'TiO2': 0.130,
        'Al2O3': 3.98,
        'Cr2O3': 0.0,
        'FeO':  8.18,
        'MgO':  38.73,
        'CaO':  3.17,
        'Na2O':  0.13,
        'K2O':  0.006,
        'P2O5':  0.019,
        'MnO': 0.0,
        'CoO': 0.0,
        'NiO': 0.0
        }

P_GPa = 2.5
Fe3_tFe = 0.1

Find the solidus

Currently, ThermoEngine cannot equilibrate subsolidus, as a different equilibration routine is required when the liquid is not present. We plan to get subsolidus equilibration implemented soon.

sys = magmaforge.System(comp=DMM,
                        T_liquidus=True,
                        P_GPa=P_GPa,
                        Fe3_tFe = Fe3_tFe,
                        database='MELTS_pMELTS',
                        )

sys.crystallize(method='equil',
                T_step=5,
                fix_fO2=False,
                calc_args={'debug':0})

solidus_T = sys.T_C

print(f"Solidus found at {solidus_T :.1f}˚C")

Solidus found at 1477.0˚C

System and Calculations

Define the system using magmaforge:

sys = magmaforge.System(comp=DMM,
                        T_C = solidus_T,
                        P_GPa=P_GPa,
                        Fe3_tFe=Fe3_tFe,
                        database='MELTS_pMELTS'
                        )

Decompress system:

sys.melt_by_decompression(P_final_GPa=0.5,
                          method='equil',
                          P_step_GPa=0.1,
                          )

<thermoengine.magmaforge.system.System object at 0xffff8f621a20>

Display Phase Mass Table

phase_masses = sys.history.get_phase_mass_table(phases='present', index='P_GPa')
#pd.set_option('display.max_columns', None)
#pd.set_option('display.max_row', None)
print(phase_masses)

       Olivine  Clinopyroxene  Orthopyroxene    Spinel    Garnet     Liquid
2.5  55.649136      13.525118      20.913526  0.000000  8.233126   0.824802
2.4  55.525459      13.669084      21.168481  0.000000  7.716916   1.065769
2.3  55.348998      13.737738      21.494497  0.000000  7.136937   1.427537
2.2  55.114522      13.699536      21.910850  0.000000  6.474116   1.946684
2.1  54.830037      13.536460      22.426934  0.000000  5.703318   2.648959
2.0  53.870528      14.954336      24.513450  1.045378  1.935719   2.826297
1.9  53.254859      15.242768      25.616947  1.504237  0.000000   3.526898
1.8  53.086921      14.183910      25.734088  1.326037  0.000000   4.814752
1.7  53.002098      12.854121      25.815185  1.104180  0.000000   6.370123
1.6  53.048190      11.250733      25.793314  0.832738  0.000000   8.220734
1.5  53.281593       9.373145      25.581419  0.504215  0.000000  10.405336
1.4  53.770558       7.223384      25.070570  0.109607  0.000000  12.971589
1.3  54.321449       5.256324      24.347241  0.000000  0.000000  15.220694
1.2  55.046559       3.300331      23.405727  0.000000  0.000000  17.393091
1.1  56.041954       1.261026      22.169815  0.000000  0.000000  19.672913
1.0  57.288526       0.000000      20.146381  0.000000  0.000000  21.710802
0.9  58.761806       0.000000      17.010127  0.000000  0.000000  23.373774
0.8  60.528232       0.000000      13.430253  0.000000  0.000000  25.187223
0.7  62.615959       0.000000       9.370475  0.000000  0.000000  27.159275
0.6  65.055139       0.000000       4.784776  0.000000  0.000000  29.305793
0.5  67.652734       0.000000       0.000000  0.000000  0.000000  31.492975

Plot System Evolution

Automatic plotting functions for decompression melting coming soon!

oxides_to_plot = ['MgO', 'SiO2', 'Al2O3', 'FeO', 'Fe2O3', 'CaO']

pressures = sys.history.get_pressures('GPa')
temperatures = sys.history.get_temperatures('C')

phase_masses = sys.history.get_phase_mass_table(phases='present', index='P_GPa')

liq_comps = sys.history.get_phase_comp_table('Liquid')

phase_colors = {'Spinel': 'black',
                'Garnet': 'purple',
                'Olivine': 'yellowgreen',
                'Clinopyroxene' : 'green',
                'Orthopyroxene': 'darkolivegreen',}

# Create plot

fig, ax = plt.subplots(1, 3, sharey='row', gridspec_kw={'wspace': 0.1, 'width_ratios': [1,1,2]}, figsize=(8,4))

# Plot phase masses
for ph in phase_masses.columns:
    if ph != 'Liquid':
        ax[0].plot(phase_masses[ph], pressures,
                   color=phase_colors[ph],label=ph)

ax[0].plot(phase_masses['Liquid'], pressures,
           lw=2, c='red', label='Melt')

# Plot temperature
ax[1].plot(temperatures, pressures, c='k', lw=2)

# Plot liquid composition
for ox in oxides_to_plot:
  if liq_comps[ox].iloc[0] > 0:
    ax[2].plot(liq_comps[ox], pressures, label=ox)

# Format axes
ax[0].invert_yaxis()
ax[0].legend(ncol=2, loc='upper left', bbox_to_anchor=(0,1.3))
ax[2].legend(ncol=2, loc='upper left', bbox_to_anchor=(0, 1.3))
plt.subplots_adjust(top=0.8)

ax[0].set_ylabel('Pressure (GPa)')
ax[0].set_xlabel('Phase Mass')
ax[1].set_xlabel('Temperature (˚C)')
ax[2].set_xlabel('Liquid Comp (wt%)')

plt.show()

Save Output

All of the output can be saved to csv files for future use, for example, uncomment the following:

# phase_fracs.to_csv('MELTS_phase_fracs.csv')
# liq_comps.to_csv('MELTS_liquid_comps.csv')