This page was generated from
docs/Examples/BarometryTests/MELTS_geobarometry_example1.ipynb.
Interactive online version:
.
MELTS geobarometry part 1 - quartz - 2-feldspar saturated rhyolites
The position of the quartz - 2-feldspar eutectic in high-silca rhyolitic magmas is sensitive to pressure (and possibly other variables; fO\(_2\), H\(_2\)O, CO\(_2\)).
If MELTS accurately recreates phase stability (and the sensitivity to the relevant parameters), it could be used as a barometer in igneous systems by searching for the pressure at which the mineral saturation curves intersect.
Before any calculations can be run users need to download the alphaMELTS for MATLAB files. Please see the installation guide on ReadTheDocs.
[1]:
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import petthermotools as ptt
# 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'MELTS')
alphaMELTS for Python files successfully located.
[3]:
# used to suppress MELTS outputs in MacOS systems (run twice)
import platform
if platform.system() == "Darwin":
import sys
import os
sys.stdout = open(os.devnull, 'w')
sys.stderr = open(os.devnull, 'w')
To demonstrate how PetThermoTools can be used to evaluate the pressure of magma storage we will use the composition of a quartz-hosted melt inclusion from Anderson et al. (2000) to evaluate the method:
[4]:
bulk = {'SiO2_Liq': 77.6,
'TiO2_Liq': 0.09,
'Al2O3_Liq': 12.3,
'FeOt_Liq': 0.65,
'MgO_Liq': 0.02,
'CaO_Liq': 0.41,
'Na2O_Liq': 4.49,
'K2O_Liq': 4.69}
In this case, we evaluate the saturation temperature of quartz, plagioclase, and k-feldspar at 32 pressures between 250 and 5000 bars. The model runs an equilibrium crystallisation calculation at each pressure (these calculations are run in parallel),with the calculations starting at the liquidus and continuing until either: (i) all three phases are saturated; or (ii) the temperature is more than 25\(^o\)C below the liquidus (T_maxdrop_C).
The function will then calculate the maximum temperature separating the saturation of the 3 phases and use a spline fit to determine the pressure of the minimum offset. This pressure is then taken forward as the pressure of storage/extraction.
[5]:
P_bar = np.linspace(250, 5000, 64)
phases = ['quartz1', 'plagioclase1', 'alkali-feldspar1']
minimum, xstal = ptt.mineral_cosaturation(
Model = "MELTSv1.0.2", bulk= bulk,
P_bar = np.linspace(250, 5000, 64),
T_initial_C= 900.0,
find_min=True, H2O_Sat=True, T_maxdrop_C=25.0,
dt_C=1.0, fO2_buffer = "NNO", T_cut_C = 5,
phases = ['quartz1', 'plagioclase1', 'alkali-feldspar1'])
To evaluate these results we can do a few different things.
First we can examine the two outputs. Here minimum is a dictionary containing two variables, one (a DataFrame called Output) contains the saturation temperature of the three phases of interest and the temperature offsets between the different mineral pairs.
[6]:
minimum['Output'].head()
[6]:
| P_bar | quartz1 | plagioclase1 | alkali-feldspar1 | quartz1 - plagioclase1 | quartz1 - alkali-feldspar1 | plagioclase1 - alkali-feldspar1 | 3 Phase Saturation | |
|---|---|---|---|---|---|---|---|---|
| 0 | 250.000000 | NaN | 903.7 | NaN | NaN | NaN | NaN | NaN |
| 1 | 325.396825 | NaN | 888.3 | NaN | NaN | NaN | NaN | NaN |
| 2 | 400.793651 | NaN | 874.9 | NaN | NaN | NaN | NaN | NaN |
| 3 | 476.190476 | NaN | 863.0 | NaN | NaN | NaN | NaN | NaN |
| 4 | 551.587302 | NaN | 852.3 | 830.3 | NaN | NaN | 22.0 | NaN |
The second variable in minimum is called CurveMin and represents a dictionary containing a series of key-value pairs, including the pressure estimated for mineral co-saturation (P_min) for each mineral-mineral pair and for 3 phase saturation (pressure here is given in bars).
[12]:
minimum['CurveMin']['quartz1 - plagioclase1']['P_min']
[12]:
2078.846614022493
The second output (xstal) contains the full results of the crystallization calculations, as you would receive if you ran these calculations through the isobaric_crystallisation function:
[8]:
xstal['Run 0']['All'].head()
[8]:
| T_C | P_bar | mass_g | H_J | S_J/K | V_cm^3 | rho_kg/m^3 | log10(fO2) | dVdP_cm^3/bar | H2O_Fl | ... | 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 | 904.7 | 250.0 | 100.015640 | -1.329559e+06 | 373.782312 | 361.478627 | 276.684796 | -11.805683 | -1.325382 | 100.0 | ... | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN |
| 1 | 903.7 | 250.0 | 100.015651 | -1.329767e+06 | 373.605310 | 361.252958 | 276.857666 | -11.823668 | -1.324726 | 100.0 | ... | -4057.577301 | 0.664751 | 1221.686553 | 0.000029 | 0.000004 | NaN | 4.255283e-10 | 9.477135e-13 | 1.065071e-12 | 264.964408 |
| 2 | 902.7 | 250.0 | 100.015662 | -1.329977e+06 | 373.426732 | 361.030413 | 277.028355 | -11.841684 | -1.324086 | 100.0 | ... | -8245.402384 | 1.350171 | 1221.597218 | 0.000059 | 0.000007 | NaN | 8.552248e-10 | 1.955436e-12 | 2.164643e-12 | 264.958944 |
| 3 | 901.7 | 250.0 | 100.015672 | -1.330187e+06 | 373.248579 | 360.806568 | 277.200254 | -11.859730 | -1.323440 | 100.0 | ... | -12395.610216 | 2.028751 | 1221.507895 | 0.000089 | 0.000011 | NaN | 1.271205e-09 | 2.984838e-12 | 3.254631e-12 | 264.953646 |
| 4 | 900.7 | 250.0 | 100.015683 | -1.330396e+06 | 373.070849 | 360.581424 | 277.373365 | -11.877807 | -1.322789 | 100.0 | ... | -16508.162110 | 2.700490 | 1221.418605 | 0.000118 | 0.000014 | NaN | 1.673374e-09 | 4.036203e-12 | 4.334996e-12 | 264.948509 |
5 rows × 87 columns
Now that we have looked out the outputs we can use a series of in-built plotting functions to visualise the results. For example, we can view the saturation surfaces of the three phases against pressure. For both of the plotting functions we need to pass the minimum variable as the Results:
[13]:
ptt.plot_surfaces(Results = minimum, P_bar = P_bar, phases = phases)
[13]:
(<Figure size 500x400 with 1 Axes>,
<Axes: xlabel='P (bars)', ylabel='T ($\\degree$C)'>)
Alternatively we can examine the residual temperatures and how these vary with pressure (along with the spline fit).
[14]:
ptt.residualT_plot(Results = minimum, P_bar = P_bar, phases = phases)
Alternatively, we can extract the pressure of magma storage directly from the ‘Results’ output variable. Doing so reveals that ‘Results’ contains a variety of information. First, it contains arrays for the saturation temperature of the three phases (at each pressure). Second, the liquidus temperature and the melt H\(_2\)O content at the liquidus is also included. Third, the temperature residuals between the saturation points of the different phases are also stored in this output. Finally, if ‘find_min = True’ has been selected in the initial function, the ‘Results’ variable will also contain the minimum pressure determined by the spline-fit to the data.
Examine the quartz saturation temperature (pressure is given by the variable P_bar):
[15]:
minimum['Output']['quartz1'].values
[15]:
array([ nan, nan, nan, nan, nan, nan, 810. , 804. , 797.7,
793.1, 788. , 784.3, 780.1, 776.3, 773.8, 770.6, 767.6, 765.9,
763.4, 761.1, 758.9, 757.9, 756.1, 753.9, 752.7, 751.7, 749.8,
749.3, 749.3, 749.4, 749.5, 749.6, 749.7, 749.8, 750. , 750.1,
750.3, 750.5, 750.7, 750.8, 751. , 751.2, 751.4, 751.6, 751.7,
751.9, 752.1, 752.3, 752.4, 752.6, 752.7, 752.9, 753. , 753.1,
753.3, 753.4, 753.5, 753.6, 753.7, 753.8, 753.9, 754. , 754.1,
754.2])
Examine the temperature residual between quartz and plagioclase saturation:
[16]:
minimum['Output']['quartz1 - plagioclase1'].values
[16]:
array([nan, nan, nan, nan, nan, nan, 24., 22., 21., 19., 18., 16., 15.,
14., 12., 11., 10., 8., 7., 6., 5., 3., 2., 0., 0., 1.,
1., 1., 3., 4., 6., 7., 9., 10., 12., 14., 15., 16., 18.,
19., 20., 22., 23., nan, nan, nan, nan, nan, nan, nan, nan, nan,
nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan])
[13]:
print('Pressure of minimum residual: ' + str(round(minimum['CurveMin']['3 Phase Saturation']['P_min'],2)) + ' bars')
print('Minimum temperature residual: ' + str(round(minimum['CurveMin']['3 Phase Saturation']['Res_min'], 2)) + ' $^o$C')
How does parallel processing affect calculation speeds?
Here we run the same calculations but investigate how parallel processing might influence the speed of the calculations by running the calculation first with only 1 process in parallel, then increasing this number to 2, 4, and 8.
[17]:
import time
[18]:
P_bar = np.linspace(250, 5000, 64)
phases = ['quartz1', 'plagioclase1', 'alkali-feldspar1']
c = np.array([1, 2, 3, 4, 6, 8])
t = np.zeros(len(c))
for i in range(len(c)):
s = time.time()
Results = ptt.mineral_cosaturation(Model = "MELTSv1.0.2",
cores = c[i],
bulk= bulk,
H2O_Sat=True,
P_bar = P_bar,
T_initial_C=900,
T_maxdrop_C=25,
dt_C=1,
fO2_buffer="NNO",
timeout = 200,
phases = phases)
t[i] = time.time() - s
[19]:
t
[19]:
array([37.09252 , 22.28453302, 20.35465813, 17.82389593, 18.2386713 ,
18.45382714])
[20]:
t[3]/t[0]
[20]:
0.4805253439761024
[16]:
t_null = 32.6 # time taken to run these calculations in the main process (i.e., no child processes or parallization)
[18]:
t
[18]:
array([33.07583618, 23.79052997, 21.82468081, 20.55444407, 19.59500694,
18.49505687])
[21]:
# load in time constraints for when this code was run on different machines
# lenovo P620
t_P620_null = 68.05
c_P620 = np.array([1,2,3,4,6,8,10,12,16,20,24,28,32])
t_P620 = np.array([67.58162308, 46.56659913, 33.45296741, 31.14283967, 31.25809503,
33.38872027, 41.64300919, 38.08958292, 44.94679379, 48.21261454,
52.22937322, 52.43954849])
[23]:
f, a = plt.subplots(3,1, figsize = (4, 8), sharex = True)
f.tight_layout()
a[0].plot(c, t, 'ob', mec = 'k', label = "multiprocessing")
a[0].plot(1, t_null, 'dr', mec = 'k', label = "Main Process Only")
# a[0].plot(c, t_P620, 'or', mec = 'k', label = "2022 Lenovo P620 - i9 12th gen")
a[0].set_ylabel('Time (secs)')
a[0].legend()
a[1].plot(c, 100*t/t_null, 'ob', mec = 'k', label = "2022 MacBook Pro - M2")
a[1].plot(1, 100*t_null/t_null, 'dr', mec = 'k', label = "2022 MacBook - Main Process Only")
# a[1].plot(c, 100*t_P620/t_P620[0], 'or', mec = 'k', label = "2022 Lenovo P620 - i9 12th gen")
a[1].set_ylabel('Run time relative to 1 process \n on each machine')
a[2].plot(c, 100*t/t[0], 'ob', mec = 'k', label = "2022 MacBook Pro M2")
# a[2].plot(c, 100*t_P620/t[0], 'or', mec = 'k', label = "2022 Lenovo P620 - i9 12th gen")
a[2].set_ylabel('Run time relative to 1 process \n on MacBook Pro')
a[2].set_xlabel('No. of Parallel Processes')
[23]:
Text(0.5, 58.7222222222222, 'No. of Parallel Processes')
[44]:
f, a = plt.subplots(1,3, figsize = (10,3.5))
a[0].set_title('Saturation Temperatures')
a[1].set_title('Temperature residuals')
a[2].set_title('Calculation time')
a[0].plot(minimum['Output']['P_bar'], minimum['Output']['quartz1'], '-k', label = 'Quartz')
a[0].plot(minimum['Output']['P_bar'], minimum['Output']['plagioclase1'], '-r', label = 'Plagioclase')
a[0].plot(minimum['Output']['P_bar'], minimum['Output']['alkali-feldspar1'], '-b', label = 'Alkali Feldspar')
a[0].legend()
a[0].set_xlabel('Pressure (bar)')
a[0].set_ylabel('Temperature ($^o$C)')
a[1].plot(minimum['Output']['P_bar'], minimum['Output']['3 Phase Saturation'], 'ok', mfc = 'b')
a[1].plot(minimum['CurveMin']['3 Phase Saturation']['P_new'], minimum['CurveMin']['3 Phase Saturation']['y_new'], '--k')
a[1].set_xlabel('Pressure (bar)')
a[1].set_ylabel('Temperature Residual ($^o$C)')
a[2].plot(c, t, 'ob', mec = 'k', label = "multiprocessing")
a[2].plot(1, t_null, 'dr', mec = 'k', label = "Main Process Only")
a[2].set_ylabel('Time (secs)')
a[2].set_xlabel('Number of processes')
a[2].legend()
f.tight_layout()
plt.savefig('SpeedComparisons.png', format = 'png', dpi = 600)
The Kernel crashed while executing code in the the current cell or a previous cell. Please review the code in the cell(s) to identify a possible cause of the failure. Click <a href='https://aka.ms/vscodeJupyterKernelCrash'>here</a> for more info. View Jupyter <a href='command:jupyter.viewOutput'>log</a> for further details.
[ ]: