This page was generated from docs/Examples/BarometryTests/MELTS_geobarometry_example3.ipynb. Interactive online version: .

Python Notebook Download

MELTS geobarometry part 3 - multi-phase saturation in mafic magmas

import numpy as np
import sys
import matplotlib.pyplot as plt

import petthermotools as ptt

ptt.__version__

alphaMELTS for Python files successfully located.

'0.3.23'

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...
    568.0 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
   1132.4 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 2:    Using libMAGEMin.dylib from MAGEMin_jll
      From worker 3:    Using libMAGEMin.dylib from MAGEMin_jll
      From worker 5:    Using libMAGEMin.dylib from MAGEMin_jll
      From worker 4:    Using libMAGEMin.dylib from MAGEMin_jll

To demonstrate how PetThermoTools can be used to estimate the pressure of storage for any magma saturated in multiple solid phases we’ll use the compsition of an olivine, clinopyroxene, and plagioclase saturated melt from experiment Y0174-12 of Neave et al. (2019). We’ll also perform these minimizations using the Green et al. (2025) thermodynamic model implemented through MAGEMin (see installation instructions for how to run MAGEMin calculations).

bulk = {'SiO2_Liq': 49.47,
       'TiO2_Liq': 1.09,
       'Al2O3_Liq': 15.58,
       'FeOt_Liq': 12.60,
       'MgO_Liq': 7.56,
       'CaO_Liq': 11.63,
       'Na2O_Liq': 1.84,
       'K2O_Liq': 0.13,
       'H2O_Liq': 1.0,
       'Fe3Fet_Liq': 0.18}

Before we run the barometry calculations there are several other things that need to be specified:

• The oxygen fugacity of the system. At the moment we cannot run these calculations at an fO2 buffer, instead the Fe redox state can be specified in the ``bulk`` variable above. This is due to a small bug that I’m trying to fix as soon as possible

• The water content (also set in the bulk variable above).

• The pressure range and number of steps.

• The phases of interest.

• The thermodynamic model we want to use for the calculations.

# fO2_offset = 0.8 # offset from FMQ buffer
P_bar = np.linspace(250.0,7500.0,24) #bars
phases = ['olivine1', 'plagioclase1', 'clinopyroxene1']
Model = "Green2025" # calculations performed using the rhyolite-MELTS v1.0.2 thermodynamic model as in Harmon et al. (2018)

Now we can use these parameters to calculate pressure. Note, with a natural glass data the H\(_2\)O content and oxygen fugacity may not be known, which leads to additional uncertainty in the calculation as these parameters have a notable influence on the saturation curves investigated here.

Notably, in the example below we’ll also specify a maximum temperaure interval on 50 degrees Celsius below the liquidus and a temperature step for each calculatio of 1 degree Celsius.

minimum, xstal = ptt.mineral_cosaturation(bulk = bulk,
                                   Model = Model,
                                   phases = phases,
                                   P_bar = P_bar,
                                   T_initial_C = 1150.0,
                                   dt_C = 1.0, # temperature interval for each calculation
                                   T_maxdrop_C = 50.0, # maximum temperature below the liquidus the code will search for 3 phase saturation
                                   timeout = 300,
                                   T_cut_C = 25.0)

Limiting calculations to 4 processes due to limitations in available memory
Computing 51 points... 100%|█████████████████████████████| Time: 0:00:09
Computing 51 points... 100%|█████████████████████████████| Time: 0:00:08
Computing 51 points... 100%|█████████████████████████████| Time: 0:00:09
Computing 51 points... 100%|█████████████████████████████| Time: 0:00:09
Computing 51 points... 100%|█████████████████████████████| Time: 0:00:08
Computing 51 points... 100%|█████████████████████████████| Time: 0:00:09
Computing 51 points... 100%|█████████████████████████████| Time: 0:00:09
Computing 51 points... 100%|█████████████████████████████| Time: 0:00:09
Computing 51 points... 100%|█████████████████████████████| Time: 0:00:08
Computing 51 points... 100%|█████████████████████████████| Time: 0:00:09
Computing 51 points... 100%|█████████████████████████████| Time: 0:00:09
Computing 51 points... 100%|█████████████████████████████| Time: 0:00:09
Computing 51 points... 100%|█████████████████████████████| Time: 0:00:10
Computing 51 points... 100%|█████████████████████████████| Time: 0:00:10
Computing 51 points... 100%|█████████████████████████████| Time: 0:00:10
Computing 51 points... 100%|█████████████████████████████| Time: 0:00:10
Computing 51 points... 100%|█████████████████████████████| Time: 0:00:09
Computing 51 points... 100%|█████████████████████████████| Time: 0:00:09
Computing 51 points... 100%|█████████████████████████████| Time: 0:00:09
Computing 51 points... 100%|█████████████████████████████| Time: 0:00:09
Computing 51 points... 100%|█████████████████████████████| Time: 0:00:09
Computing 51 points... 100%|█████████████████████████████| Time: 0:00:09
Computing 51 points... 100%|█████████████████████████████| Time: 0:00:09
Computing 51 points... 100%|█████████████████████████████| Time: 0:00:10

To investigate the results, we can use in-built plotting functions to look at the position of the saturation curves in P-T space:

ptt.plot_surfaces(Results = minimum, P_bar = P_bar, phases = phases)

(<Figure size 500x400 with 1 Axes>,
 <Axes: xlabel='P (bars)', ylabel='T ($\\degree$C)'>)

We can also investigate the temperature residual between the curves using the plottin function below.

ptt.residualT_plot(Results = minimum, P_bar = P_bar, phases = phases, interpolate = True, ylim = [0, 100])