ENKI

pMELTS

Versions of MELTS implemented are:
- MELTS v. 1.0.2 ➞ (rhyolite-MELTS, Gualda et al., 2012)
- MELTS v. 1.1.0 ➞ (rhyolite-MELTS + new CO2, works at the ternary minimum)
- MELTS v. 1.2.0 ➞ (rhyolite-MELTS + new H2O + new CO2)
- pMELTS v. 5.6.1

Initialize tools and packages that are required to execute this notebook.

from thermoengine import equilibrate
import matplotlib.pyplot as plt
import numpy as np
%matplotlib inline

Create a pMELTS v 5.6.1 instance.

Rhyolite-MELTS version 1.0.2 is the default model.

melts = equilibrate.MELTSmodel(version="5.6.1")

Optional: Generate some information about the implemented model.

oxides = melts.get_oxide_names()
phases = melts.get_phase_names()
#print (oxides)
#print (phases)

Required: Input initial composition of the system (liquid), in wt% or grams of oxides.

Early Bishop Tuff average melt inlusion composition

feasible = melts.set_bulk_composition({'SiO2':  45.47,
                                       'TiO2':   0.11,
                                       'Al2O3':  4.0,
                                       'Fe2O3':  0.585,
                                       'Cr2O3':  0.0,
                                       'FeO':    6.696,
                                       'MnO':    0.0,
                                       'MgO':   38.53,
                                       'NiO':    0.0,
                                       'CoO':    0.0,
                                       'CaO':    3.59,
                                       'Na2O':   0.31,
                                       'K2O':    0.0,
                                       'P2O5':   0.0,
                                       'H2O':    0.0})

Optional: Suppress phases that are not required in the simulation.

b = melts.get_phase_inclusion_status()
melts.set_phase_inclusion_status({'Actinolite':False, 'Aegirine':False, \
                                  'Aenigmatite':False, 'Akermanite':False, 'Andalusite':False, \
                                  'Anthophyllite':False, 'Apatite':False, 'Biotite':False, 'Chromite':False, \
                                  'Coesite':False, 'Corundum':False, 'Cristobalite':False, 'Cummingtonite':False, \
                                  'Fayalite':False, 'Forsterite':False, 'Gehlenite':False, 'Hematite':False, \
                                  'Hornblende':False, 'Ilmenite':False, 'Ilmenite ss':False, 'Kalsilite':False, \
                                  'Kalsilite ss':False, 'Kyanite':False, 'Leucite':False, 'Lime':False, \
                                  'Liquid Alloy':False, 'Magnetite':False, 'Melilite':False, 'Muscovite':False, \
                                  'Nepheline':False, 'Nepheline ss':False, 'OrthoOxide':False, 'Panunzite':False, \
                                  'Periclase':False, 'Perovskite':False, 'Phlogopite':False, 'Quartz':False, \
                                  'Rutile':False, 'Sanidine':False, 'Sillimanite':False, 'Solid Alloy':False, \
                                  'Sphene':False, 'Tridymite':False, 'Whitlockite':False})

a = melts.get_phase_inclusion_status()
for phase in b.keys():
    if b[phase] != a[phase]:
        print ("{0:<15s} Before: {1:<5s} After: {2:<5s}".format(phase, repr(b[phase]), repr(a[phase])))
Actinolite      Before: True  After: False
Aegirine        Before: True  After: False
Aenigmatite     Before: True  After: False
Akermanite      Before: True  After: False
Andalusite      Before: True  After: False
Anthophyllite   Before: True  After: False
Apatite         Before: True  After: False
Biotite         Before: True  After: False
Chromite        Before: True  After: False
Coesite         Before: True  After: False
Corundum        Before: True  After: False
Cristobalite    Before: True  After: False
Cummingtonite   Before: True  After: False
Fayalite        Before: True  After: False
Forsterite      Before: True  After: False
Gehlenite       Before: True  After: False
Hematite        Before: True  After: False
Hornblende      Before: True  After: False
Ilmenite        Before: True  After: False
Ilmenite ss     Before: True  After: False
Kalsilite       Before: True  After: False
Kalsilite ss    Before: True  After: False
Kyanite         Before: True  After: False
Leucite         Before: True  After: False
Lime            Before: True  After: False
Liquid Alloy    Before: True  After: False
Magnetite       Before: True  After: False
Melilite        Before: True  After: False
Muscovite       Before: True  After: False
Nepheline       Before: True  After: False
Nepheline ss    Before: True  After: False
OrthoOxide      Before: True  After: False
Panunzite       Before: True  After: False
Periclase       Before: True  After: False
Perovskite      Before: True  After: False
Phlogopite      Before: True  After: False
Quartz          Before: True  After: False
Rutile          Before: True  After: False
Sanidine        Before: True  After: False
Sillimanite     Before: True  After: False
Solid Alloy     Before: True  After: False
Sphene          Before: True  After: False
Tridymite       Before: True  After: False
Whitlockite     Before: True  After: False

Compute the equilibrium state at some specified T (°C) and P (MPa).

Print status of the calculation.

output = melts.equilibrate_tp(1300.0, 1000.0, initialize=True)
(status, t, p, xmlout) = output[0]
print (status, t, p)
success, Optimal residual norm. 1300.0 1000.0

Summary output of equilibrium state …

melts.output_summary(xmlout)
T (°C)      1300.00
P (MPa)     1000.00
Augite           15.6279 (g)  Na0.05Ca0.67Fe''0.10Mg0.97Fe'''0.04Ti0.01Al0.32Si1.84O6
Olivine          53.9240 (g)  (Ca0.00Mg0.91Fe''0.09Mn0.00Co0.00Ni0.00)2SiO4
Spinel            1.1339 (g)  Fe''0.12Mg0.88Fe'''0.10Al1.89Cr0.00Ti0.00O4
Liquid            4.0410 (g)  wt%:SiO2 48.94 TiO2  0.48 Al2O3 19.05 Fe2O3  1.96 Cr2O3  0.00 FeO  5.04 MnO  0.00 MgO
                                  11.73 NiO  0.00 CoO  0.00 CaO  8.39 Na2O  4.42 K2O  0.00 P2O5  0.00 H2O  0.00
Orthopyroxene    24.5641 (g)  Na0.01Ca0.06Fe''0.15Mg1.67Fe'''0.02Ti0.00Al0.22Si1.88O6

Compute the equilibrium state at some specified S (J/K) and P (MPa).

s = melts.get_property_of_phase(xmlout,'System', 'Entropy')
print ("{0:<20s} {1:13.6e} {2:<10s}".format('Entropy', s, melts.get_units_of_property('Entropy')))
output = melts.equilibrate_sp(s, p, initialize=True)
(status, t, p, xmlout) = output[0]
melts.output_summary(xmlout)
Entropy               2.505416e+02 J/K
T (°C)      1300.00
P (MPa)     1000.00
Augite           15.6279 (g)  Na0.05Ca0.67Fe''0.10Mg0.97Fe'''0.04Ti0.01Al0.32Si1.84O6
Olivine          53.9240 (g)  (Ca0.00Mg0.91Fe''0.09Mn0.00Co0.00Ni0.00)2SiO4
Spinel            1.1339 (g)  Fe''0.12Mg0.88Fe'''0.10Al1.89Cr0.00Ti0.00O4
Liquid            4.0410 (g)  wt%:SiO2 48.94 TiO2  0.48 Al2O3 19.05 Fe2O3  1.96 Cr2O3  0.00 FeO  5.04 MnO  0.00 MgO
                                  11.73 NiO  0.00 CoO  0.00 CaO  8.39 Na2O  4.42 K2O  0.00 P2O5  0.00 H2O  0.00
Orthopyroxene    24.5641 (g)  Na0.01Ca0.06Fe''0.15Mg1.67Fe'''0.02Ti0.00Al0.22Si1.88O6

Run the sequence of calculations along a S=constant, P path:

Output is sent to an Excel file and plotted in the notebook

number_of_steps = 20
s_increment_of_steps = 0.0
p_increment_of_steps = -50.0

plotPhases = ['Liquid', 'Olivine', 'Orthopyroxene', 'Augite', 'Plagioclase']
# matplotlib colors b : blue, g : green, r : red, c : cyan, m : magenta, y : yellow, k : black, w : white.
plotColors = [ 'ro', 'bo', 'go', 'co', 'mo']

wb = melts.start_excel_workbook_with_sheet_name(sheetName="Summary")
melts.update_excel_workbook(wb, xmlout)

n = len(plotPhases)
xPlot = np.zeros(number_of_steps)
yPlot = np.zeros((n, number_of_steps))
y2Plot = np.linspace(t, t-100.0, number_of_steps)
y3Plot = np.linspace(t, t-100.0, number_of_steps)
xPlot[0] = p
for i in range (0, n):
    yPlot[i][0] = melts.get_property_of_phase(xmlout, plotPhases[i])
y2Plot[0] = t

plt.ion()
fig = plt.figure()
ax = fig.add_subplot(111)
ax.set_xlim([min(p, p+p_increment_of_steps*number_of_steps), max(p, p+p_increment_of_steps*number_of_steps)])
ax.set_ylim([0., 100.])
ax2 = ax.twinx()
ax2.set_ylabel('T', color='k')
ax2.set_ylim([t-100, t])

graphs = []
for i in range (0, n):
    graphs.append(ax.plot(xPlot, yPlot[i], plotColors[i]))
graphs.append(ax2.plot(xPlot, y2Plot, 'k-'))
graphs.append(ax2.plot(xPlot, y3Plot, 'k--'))
handle = []
for (graph,) in graphs:
    handle.append(graph)
ax.legend(handle, plotPhases, loc='upper left')

for i in range (1, number_of_steps):
    output = melts.equilibrate_sp(s+s_increment_of_steps, p+p_increment_of_steps)
    (status, t, p, xmlout) = output[0]
    print ("{0:<30s} {1:8.2f} {2:8.2f}".format(status, t, p))
    xPlot[i] = p
    for j in range (0, n):
        yPlot[j][i] = melts.get_property_of_phase(xmlout, plotPhases[j])
    y2Plot[i] = t
    j = 0
    for (graph,) in graphs:
        graph.set_xdata(xPlot)
        if j < n:
            graph.set_ydata(yPlot[j])
        elif j == n:
            graph.set_ydata(y2Plot)
        j = j + 1
    fig.canvas.draw()
    melts.update_excel_workbook(wb, xmlout)

melts.write_excel_workbook(wb, "MELTSv102summary.xlsx")
success, Optimal residual norm.  1295.97   950.00
success, Optimal residual norm.  1291.41   900.00
success, Optimal residual norm.  1286.25   850.00
success, Optimal residual norm.  1281.49   800.00
success, Optimal residual norm.  1277.55   750.00
success, Optimal residual norm.  1273.18   700.00
success, Optimal residual norm.  1268.41   650.00
success, Optimal residual norm.  1263.21   600.00
success, Optimal residual norm.  1257.61   550.00
success, Optimal residual norm.  1251.57   500.00
success, Optimal residual norm.  1245.10   450.00
success, Optimal residual norm.  1238.75   400.00
success, Optimal residual norm.  1234.03   350.00
success, Optimal residual norm.  1229.33   300.00
success, Optimal residual norm.  1224.60   250.00
success, Optimal residual norm.  1219.79   200.00
success, Optimal residual norm.  1218.15   150.00
success, Optimal residual norm.  1216.97   100.00
success, Optimal residual norm.  1215.78    50.00
_images/pMELTS-v5.6.1-adiabatic_18_1.png