This page was generated from docs/Examples/PhaseDiagrams/PhaseDiagram_example_MaunaLoa.ipynb. Interactive online version: .

Phase Diagram calculations at Mauna Loa

Given a set composition we can use ptt to assess the phase stability at different P-T conditions to construct a phase diagram.
These calculations can be performed with either the MELTS or Holland & Powell family of thermodynamic models (through MAGEMin).
Also use outputs to assess the variability in key geochemical parameters across P-T space (e.g., phase proportions or chemistry).

Before any calculations can be run users need to download the alphaMELTS for MATLAB files and/or install the MAGEMinCalc and juliacall packages. Please see the installation guide on ReadTheDocs.

[1]:

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
# If the alphaMELTS for Python files have not been added to your Python path (see installation guide) then use the two lines below to add
# the location of the alphaMELTS files here.
# import sys
# sys.path.append(r'C:\Users\penny\Box\Berkeley_new\MELTS_Installation\alphamelts_py\alphamelts-py-2.3.1-win64')


import petthermotools as ptt

print(ptt.__version__)
# If you are running in jupyter lab or notebook, remove the word widget.
%matplotlib widget

alphaMELTS for Python files successfully located.
0.3.23

[2]:

ptt.activate_petthermotools_env()

Julia environment detected.
Detected IPython. Loading juliacall extension. See https://juliapy.github.io/PythonCall.jl/stable/compat/#IPython

  Activating project at `~/.petthermotools_julia_env`
Precompiling MAGEMinCalc...
    600.5 ms  ? PythonCall
Info Given MAGEMinCalc was explicitly requested, output will be shown live
Using libMAGEMin.dylib from MAGEMin_jll
┌ Warning: Module PythonCall with build ID fafbfcfd-3501-17ed-0393-e686105ae5fb is missing from the cache.
│ This may mean PythonCall [6099a3de-0909-46bc-b1f4-468b9a2dfc0d] does not support precompilation but is imported by a module that does.
└ @ Base loading.jl:2541
   1191.0 ms  ? MAGEMinCalc
[ Info: Precompiling MAGEMinCalc [dd69d18c-a947-4316-a873-d8ea006ff1b2] (cache misses: wrong dep version loaded (2), mismatched flags (2))
┌ Warning: Module PythonCall with build ID fafbfcfd-3501-17ed-0393-e686105ae5fb is missing from the cache.
│ This may mean PythonCall [6099a3de-0909-46bc-b1f4-468b9a2dfc0d] does not support precompilation but is imported by a module that does.
└ @ Base loading.jl:2541
┌ Info: Skipping precompilation due to precompilable error. Importing MAGEMinCalc [dd69d18c-a947-4316-a873-d8ea006ff1b2].
└   exception = Error when precompiling module, potentially caused by a __precompile__(false) declaration in the module.

Using libMAGEMin.dylib from MAGEMin_jll
Using libMAGEMin.dylib from MAGEMin_jll
      From worker 5:    Using libMAGEMin.dylib from MAGEMin_jll
      From worker 2:    Using libMAGEMin.dylib from MAGEMin_jll
      From worker 4:    Using libMAGEMin.dylib from MAGEMin_jll
      From worker 3:    Using libMAGEMin.dylib from MAGEMin_jll

As in all examples we will import the key packages including ptt. The only difference here is that we will also use %matplotlib widget to make the phase diagram plots interactive (more explanation later).

The cell below can be used for users running MacOS to suppress the outputs (you need to run it twice). I’m still trying to find a better way to deadl with this - if you have suggestions please let me know.

[4]:

# Should be used to suppress outputs on MacOS - run twice - if on windows, wont do anything.
import platform
if platform.system() == "Darwin":
    import sys
    import os
    sys.stdout = open(os.devnull, 'w')
    sys.stderr = open(os.devnull, 'w')

Bulk composition

We are going to recreate the study of Wieser et al. (2025b), which places constraints on magma storage at Mauna Loa Volcano (Bull Volcanol, https:// doi. org/ 10. 1007/ s00445- 025- 01869-2). It has been suggested based on pMELTS modelling that the presence of high Mg# Opx indicates sub-Moho magma storage. We can assess if this is true using different thermodynamic models.

[5]:

# Estimated primary melt composition using reverse crystallization calculations in Petrolog3
Prim_MELT = {
    'SiO2_Liq': 48.5858,
    'TiO2_Liq': 1.4381,
    'Al2O3_Liq': 9.73,
    'FeOt_Liq': 10.144,        # 9.064 + 0.8998*1.1996
    'MgO_Liq': 20.0,
    'CaO_Liq': 7.2186,
    'Na2O_Liq': 1.7327,
    'K2O_Liq': 0.3227,
    'Cr2O3_Liq': 0.0702,
    'H2O_Liq': 0.2,
    'Fe3Fet_Liq': 0.106       # (0.8998*1.1996) / FeOt
}

We can now run ptt.phaseDiagram_calc which takes an array of pressure and temperature conditions and evaluates the equilibrium assemblage at each combination of values. The import parameters to tweak are:

The i_max parameters determines the maximum number of iterations that the code will loop through (a new loop is started if one calculation fails) unless every point either has a solution or hits a failure point. Smaller i_max, more likely to have holes in the phase diagram (as if one calculation kills a subprocess with 7 more calculations to run, you will get 7 holes until you have i_max>1). 15 is generally a good value to use.
The refine = 2 parameter helps refine the phase boundaries. Consider that your 25 pressure and temperature points creates the grid shown below. Initially, the phase boundary between olivine and liquid is put at the midpoint between the solved values. What the first refine step does is insert a new calculation anywhere two neighbouring points have a different phase assemble (the red dots). These are then solved, and the position of the phase boundary is obtained at a higher precision. To ensure that the matrix has a constant point spacing, new points are inserted between existing points with the same phase (e.g. the point between Ol - Ol is filled with olivine - in this diagram, these are the smaller circles in panel c). refine=2 would do this again after the first refinement stage.

[6]:

Res_pMELTS = ptt.phaseDiagram_calc(Model = "pMELTS", bulk = Prim_MELT,
                            P_bar = np.linspace(200.0,15000.0, 25),
                            T_C = np.linspace(1150.0, 1450.0, 25),
                            i_max = 15,
                            refine = 2)

[7]:

Res_pMELTS

[7]:

	T_C	P_bar	mass_g	H_J	S_J/K	V_cm^3	rho_kg/m^3	log10(fO2)	dVdP_cm^3/bar	SiO2_Cpx	...	H_J_Plag	S_J/K_Plag	Cp_J/(kg.K^2)_Plag	dCpdT_J/(kg.K^2)_Plag	dVdT_cm^3/K_Plag	dPdT_bar/K_Plag	d2VdT2_cm^3/K^2_Plag	d2VdTdP_cm^3/(bar.K)_Plag	d2VdP2_cm^3/bar^2_Plag	molwt_Plag
0	1150.0	200.00	100.0	-1.248859e+06	255.862662	35.541637	2813.601422	-8.208987	0.000353	NaN	...	-4055.510671	0.701873	1227.622241	0.000035	0.000003	NaN	8.456372e-10	0.0	7.925588e-13	273.791174
1	1150.0	354.17	100.0	-1.248532e+06	255.708136	35.507705	2816.297855	-8.175625	0.000360	NaN	...	-9315.357523	1.615936	1227.626545	0.000079	0.000006	NaN	1.940876e-09	0.0	1.835080e-12	273.722894
2	1150.0	508.33	100.0	-1.248206e+06	255.553609	35.473774	2818.994287	-8.142263	0.000368	NaN	...	-14575.204374	2.529999	1227.630849	0.000124	0.000009	NaN	3.036115e-09	0.0	2.877602e-12	273.654613
3	1150.0	662.50	100.0	-1.247879e+06	255.399083	35.439842	2821.690719	-8.108900	0.000376	NaN	...	-19835.051226	3.444063	1227.635154	0.000168	0.000013	NaN	4.131354e-09	0.0	3.920124e-12	273.586332
4	1150.0	816.67	100.0	-1.247552e+06	255.244556	35.405911	2824.387152	-8.075538	0.000383	NaN	...	-25094.898077	4.358126	1227.639458	0.000213	0.000016	NaN	5.226593e-09	0.0	4.962645e-12	273.518051
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
9404	1450.0	14383.33	100.0	-1.144643e+06	290.385564	34.672850	2884.100960	-5.168325	0.000004	NaN	...	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN
9405	1450.0	14537.50	100.0	-1.144694e+06	290.046490	34.619078	2888.601626	-5.153195	0.000019	NaN	...	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN
9406	1450.0	14691.67	100.0	-1.144745e+06	289.707416	34.565307	2893.102291	-5.138066	0.000033	NaN	...	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN
9407	1450.0	14845.83	100.0	-1.144797e+06	289.368342	34.511535	2897.602957	-5.122936	0.000048	NaN	...	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN
9408	1450.0	15000.00	100.0	-1.144848e+06	289.029268	34.457763	2902.103623	-5.107807	0.000063	NaN	...	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN

9409 rows × 219 columns

The result is a DataFrame with 97 equally spaced temperature and pressure values (so total size of 97*97). Why 97 you ask? Well initially its 25 in x and 25 in y. The refine then adds 24 points (one inbetween each point), so 49X49. The refine again adds another 48 points so (49+48) = 97 X 97. For cases where the phase assemblage of 2 neighboring points was identical the phase proportions and phase chemistry is estimated through simply linear interpolation between the two neighbors. To view the phase diagram results we can use the plot_phaseDiagram function in ptt. This will create a phase diagram labelled by number, with each number corresponding to a phase assemblage listed on the right. By using %matplotlib widget above you should be able to hover over the plot and see which number (and therefore phase assemblage) is present at each P-T condition.

[8]:

ptt.plot_phaseDiagram(Combined = Res_pMELTS)

[8]:

(<Figure size 1000x600 with 2 Axes>,
 array([<Axes: xlabel='Temperature ($\\degree$C)', ylabel='Pressure (bar)'>,
        <Axes: >], dtype=object))

Lets try Melts V.1.0.2

[9]:

Res_M12 = ptt.phaseDiagram_calc(Model = "MELTSv1.0.2", bulk = Prim_MELT,
                            P_bar = np.linspace(200.0,15000.0, 25),
                            T_C = np.linspace(1150.0, 1450.0, 25),
                            i_max = 15,
                            refine = 2)

[10]:

f, a = ptt.plot_phaseDiagram(Combined = Res_M12)

Calculations using MAGEMin

For calculations with MAGEMin it is important to note that this is one of the major purposes for which the original MAGEMin code was developed, so it will likely always be straightforward to create phase diagrams using the MAGEMinApp (https://github.com/ComputationalThermodynamics/MAGEMinApp.jl). The point of this code is not to try and replace MAGEMin, but to provide Python users with the ability to perform calcualtions with a range of different thermodynamic models. In the cells below we repeat the calculations with the Model kwarg is changed to Green2025 which will initiate the calculation in julia via MAGEMin

[11]:

Res_green = ptt.phaseDiagram_calc(Model = "Green2025", bulk = Prim_MELT,
                            P_bar = np.linspace(200.0,15000.0, 25),
                            T_C = np.linspace(1150.0, 1450.0, 25),
                            i_max = 15, refine=2)

[12]:

f, a = ptt.plot_phaseDiagram(Combined = Res_green)

[13]:

f, a = ptt.plot_phaseDiagram_multi({"Green 2025": Res_green,
                                    "MELTSv1.0.2": Res_M12,
                                    "pMELTS": Res_pMELTS})

Lets calculate Mg# as that was a big part of the original arguements

[14]:

Res_green['Mg#_Opx']=(Res_green['MgO_Opx']/40.3044)/((Res_green['MgO_Opx']/40.3044)+(Res_green['FeOt_Opx']/71.844))
Res_M12['Mg#_Opx']=(Res_M12['MgO_Opx']/40.3044)/((Res_M12['MgO_Opx']/40.3044)+(Res_M12['FeOt_Opx']/71.844))
Res_pMELTS['Mg#_Opx']=(Res_pMELTS['MgO_Opx']/40.3044)/((Res_pMELTS['MgO_Opx']/40.3044)+(Res_pMELTS['FeOt_Opx']/71.844))

In addition to the phase diagrams shown above it is relatively simple to extract additional information from the calculations. For example, it is trivial to extract the mass of liquid (mass fraction) at each point and reshape this array into a P-T grid so that we can visualize this data as contours.

The function make_grid_variables_from_phase_diagram takes your dataframe output, and puts it into a PT grid. This allows for easier plotting.

[15]:

X_green, Y_green, Z_green = ptt.make_grid_variables_from_phase_diagram(
    Results = Res_green,
    x_col = "T_C",
    y_col = "P_bar",
)

liq_green = Z_green["mass_g_Liq"]

[16]:

X_pMELTS, Y_pMELTS, Z_pMELTS = ptt.make_grid_variables_from_phase_diagram(
    Results = Res_pMELTS,
    x_col = "T_C",
    y_col = "P_bar",
)

[17]:

X_M12, Y_M12, Z_M12 = ptt.make_grid_variables_from_phase_diagram(
    Results = Res_M12,
    x_col = "T_C",
    y_col = "P_bar",
)

From here we can use the plt.contourf function to visualize the mass fraction of melt at each P-T value:

[18]:

import matplotlib.pyplot as plt

fig, (ax1, ax2, ax3) = plt.subplots(
    1, 3, figsize=(12, 4), sharex=True, sharey=True
)

# ------------------------------------------------
# Panel 1: Opx in Green
pcm1 = ax1.pcolormesh(
    X_green, Y_green, Z_green["mass_g_Opx"],
    shading="auto",
    cmap="viridis", vmin=0, vmax=10,
)
ax1.set_title("Green2025")
ax1.set_xlabel("T (°C)")
ax1.set_ylabel("P (bar)")
fig.colorbar(pcm1, ax=ax1, label="% Opx")

# ------------------------------------------------
# Panel 2: Opx in pMELTS

pcm2 = ax2.pcolormesh(
    X_pMELTS, Y_pMELTS, Z_pMELTS["mass_g_Opx"],
    shading="auto",
    cmap="viridis", vmin=0, vmax=10,
)
ax2.set_title("pMELTS")
ax2.set_xlabel("T (°C)")
fig.colorbar(pcm2, ax=ax2, label="% Opx")

# ------------------------------------------------
# Panel 3: Opx in v1.0.2

pcm3 = ax3.pcolormesh(
    X_M12, Y_M12, Z_M12["mass_g_Opx"],
    shading="auto",
    cmap="viridis", vmin=0, vmax=10,
)
ax3.set_title("MELTS 1.0.2")
ax3.set_xlabel("T (°C)")
fig.colorbar(pcm3, ax=ax3, label="% Opx")

plt.tight_layout()
plt.show()

We can also plot the P-T conditions where Ol and Opx are stable (the observed phase assemblage)

[19]:

Res_green_Opx_Ol=Res_green.loc[(Res_green['mass_g_Ol']>0) & (Res_green['mass_g_Opx']>0) ]
Res_green_Opx=Res_green.loc[(Res_green['mass_g_Opx']>0) ]


Res_pMELTS_Opx_Ol=Res_pMELTS.loc[(Res_pMELTS['mass_g_Ol']>0) & (Res_pMELTS['mass_g_Opx']>0 )]
Res_pMELTS_Opx=Res_pMELTS.loc[(Res_pMELTS['mass_g_Opx']>0 )]


Res_M12_Opx_Ol=Res_M12.loc[(Res_M12['mass_g_Ol']>0) & (Res_M12['mass_g_Opx']>0 )]
Res_M12_Opx=Res_M12.loc[ (Res_M12['mass_g_Opx']>0 )]

[20]:

fig, (ax1, ax2) = plt.subplots(
    1, 2, figsize=(12, 4), sharex=True, sharey=True
)

# ------------------------------------------------
# Panel 1: Opx + Ol in Green
ax1.set_title('Ol + Opx stable')

ax1.plot(Res_green_Opx_Ol['T_C'], Res_green_Opx_Ol['P_bar'], '.r', label='Green2025')

ax1.plot(Res_pMELTS_Opx_Ol['T_C'], Res_pMELTS_Opx_Ol['P_bar'], '.k', alpha=0.3, label='pMELTS')

ax1.plot(Res_M12_Opx_Ol['T_C'], Res_M12_Opx_Ol['P_bar'], '.c', alpha=0.1, label='1.0.2')

# Lets just do Opx stable
# Panel 2: Opx stable
ax2.set_title('Opx stable')

ax2.plot(Res_green_Opx['T_C'], Res_green_Opx['P_bar'], '.r')

ax2.plot(Res_pMELTS_Opx['T_C'], Res_pMELTS_Opx['P_bar'], '.k', alpha=0.2)
ax2.plot(Res_M12_Opx['T_C'], Res_M12_Opx['P_bar'], '.c', alpha=0.1)
ax1.legend()
ax1.set_xlabel('T (C)')
ax1.set_ylabel('P (bar)')

[20]:

Text(0, 0.5, 'P (bar)')

[21]:

fig, (ax1, ax2) = plt.subplots(
    1, 2, figsize=(12, 4), sharex=True, sharey=True
)

# ------------------------------------------------
# Panel 1: Opx + Ol in Green
ax1.set_title('Ol + Opx stable')

ax1.plot(Res_green_Opx_Ol['Mg#_Opx'], Res_green_Opx_Ol['P_bar'], '.r', label='Green2025')

ax1.plot(Res_pMELTS_Opx_Ol['Mg#_Opx'], Res_pMELTS_Opx_Ol['P_bar'], '.k', alpha=0.2, label='pMELTS')
#ax1.plot(Res_M12_Opx_Ol['Mg#_Opx'], Res_M12_Opx_Ol['P_bar'], '.c', alpha=0.2)
ax1.legend()
# Lets just do Opx stable
# Panel 2: Opx stable
ax2.set_title('Opx stable')

ax2.plot(Res_green_Opx['Mg#_Opx'], Res_green_Opx['P_bar'], '.r')

ax2.plot(Res_pMELTS_Opx['Mg#_Opx'], Res_pMELTS_Opx['P_bar'], '.k', alpha=0.2)
#ax2.plot(Res_M12_Opx['Mg#_Opx'], Res_M12_Opx['P_bar'], '.c', alpha=0.1)
ax1.legend()
ax1.set_xlabel('Mg#_Opx')
ax1.set_ylabel('P (bar)')

[21]:

Text(0, 0.5, 'P (bar)')

[22]:

fig, axes = plt.subplots(
    3, 2, figsize=(12, 10), sharex=True, sharey=True
)

# =================================================
# Row 1 — Green2025
# =================================================
axes[0, 0].set_title('Ol + Opx stable')
axes[0, 0].plot(
    Res_green_Opx_Ol['Mg#_Opx'],
    Res_green_Opx_Ol['P_bar'],
    'xr', ms=2, label='Green2025'
)
axes[0, 0].legend()

axes[0, 1].set_title('Opx stable')
axes[0, 1].plot(
    Res_green_Opx['Mg#_Opx'],
    Res_green_Opx['P_bar'],
    'xr', ms=2
)

# =================================================
# Row 2 — pMELTS
# =================================================
axes[1, 0].plot(
    Res_pMELTS_Opx_Ol['Mg#_Opx'],
    Res_pMELTS_Opx_Ol['P_bar'],
    '.k', alpha=0.2, label='pMELTS'
)
axes[1, 0].legend()

axes[1, 1].plot(
    Res_pMELTS_Opx['Mg#_Opx'],
    Res_pMELTS_Opx['P_bar'],
    '.k', alpha=0.2
)

# =================================================
# Row 3 — MELTS v1.2.0
# =================================================
axes[2, 0].plot(
    Res_M12_Opx_Ol['Mg#_Opx'],
    Res_M12_Opx_Ol['P_bar'],
    '.c', alpha=0.3, label='v1.0.2'
)
axes[2, 0].legend()

axes[2, 1].plot(
    Res_M12_Opx['Mg#_Opx'],
    Res_M12_Opx['P_bar'],
    '.c', alpha=0.3
)

# =================================================
# Labels
# =================================================
for ax in axes[:, 0]:
    ax.set_ylabel('P (bar)')

for ax in axes[-1, :]:
    ax.set_xlabel('Mg#$_{Opx}$')

plt.tight_layout()

[ ]: