This page was generated from docs/Examples/LiquidusTests/LiquidusTest.ipynb. Interactive online version: .

Python Notebook Download

Using MELTS to find the liquidus temperature

• Finding the liquidus temperature of a magma is one of the first steps of most crystallisation and/or decompression calculations.

• Usually this command is used as part of the isobaric_crystallisation and isothermal_decompression functions (amongst others).

• However, it’s still useful to consider how the liquidus function works independently, and examine how well MELTS performs as a Thermometer in igneous systems.

Before any calculations can be run users need to download the alphaMELTS for MATLAB files (https://magmasource.caltech.edu/gitlist/MELTS_Matlab.git/) and store them locally on their computer. These files then need to be added to the Python path by using the \(sys.path.append()\) command below.

Data used in the calculations below can be downloaded from here: https://github.com/gleesonm1/pyMELTScalc/blob/master/docs/Examples/LiquidusTests/LEPR_Wet_Stitched_Sept2022.xlsx

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import pyMELTScalc as M

import sys
sys.path.append(r'MELTS')

To test the ability of MELTS to act as a liquid-only thermometer we need some experiments to test them against. Here, we use a database of experimental data compiled by P. Wieser and isolate all the melt compositions in the \(Liq\) sheet.

Data = pd.read_excel('LEPR_Wet_stitched_Sept2022.xlsx', sheet_name = 'Liq')

The entire experimental dataset contains over 2000 samples. Therefore, to reduce the computational time in this example spreadsheet I isolate 200 experiments (randomly selected) to be used in the following calculations. A new DataFrame is constructed using only these experiments.

# randomly select 200 rows from the original DataFrame
Ch = np.random.choice(range(len(Data['SiO2_Liq'])), 2000, replace=False)
Test = Data.copy()
Test = Test.loc[Ch]

# reset the index in the new DataFrame
Test = Test.reset_index(drop = True)

# used to suppress MELTS outputs in MacOS systems (run twice)
import os
sys.stdout = open(os.devnull, 'w')
sys.stderr = open(os.devnull, 'w')

Now that we have this ‘Test’ DataFrame we can use the findLiq_multi function to identify the liquidus temperature of each sample. Ideally, we’d have information about either the Fe redox state or an estimate of the oxygen fugacity of each melt. Unfortunately, we don’t have this information for every sample so instead we will define a constant Fe\(^{3+}\)/Fe\(_{tot}\) ratio for all samples:

Results_pMELTS = M.findLiq_multi(Model = "pMELTS", # which MELTS model to use - "pMELTS", "MELTSv1.0.2", "MELTSv1.1.0", or "MELTSv1.2.0"
                        bulk = Test, # initial composition(s) wither a dictionary (if one composition) or a DataFrame (multiple compositions)
                        T_initial_C = Test['T_K'].values-200, # Initial guess for the liquidus temperature (optional)
                        P_bar = Test['P_kbar'].values*1000, # Pressure (in bars) of the calculation
                        Fe3Fet_Liq = 0.15) # initial Fe3Fet_Liq ratio of the sample(s)

Examining the results demonstrates the information that is recorded by this code. Notably, the findLiq_multi code will return the composition of the melt at the liquidus (in wt% hydrous normalized), whether or not the melt is fluid saturated at the liquidus, the liquidus phase and (critically) the liquidus temperature (in degrees Celsius).

Results_pMELTS.head()

	T_Liq	liquidus_phase	fluid_saturated	SiO2_Liq	TiO2_Liq	Al2O3_Liq	Fe2O3	Cr2O3	FeO	MnO_Liq	MgO_Liq	CaO_Liq	Na2O_Liq	K2O_Liq	P2O5_Liq	H2O_Liq	CO2_Liq	FeOt_Liq	Fe3Fet_Liq
0	1157.35	olivine1	No	52.997792	0.670484	16.555032	1.279656	0.0	6.528835	0.108461	6.438616	10.106556	3.766541	0.069020	0.000000	1.479008	0.0	7.680261	0.14992
1	980.75	clinopyroxene1	No	61.938726	0.539009	15.981137	0.531496	0.0	2.711704	0.113476	2.108754	4.822710	4.321526	1.219862	0.151301	5.560301	0.0	3.189940	0.14992
2	885.75	biotite1	No	67.007141	0.343422	13.375135	0.386015	0.0	1.969458	0.088343	0.509983	2.072056	4.489019	2.866449	0.000000	6.892978	0.0	2.316791	0.14992
3	902.75	plagioclase1	No	70.409287	0.292452	13.346339	0.282837	0.0	1.443042	0.155660	0.176578	0.859103	4.320344	4.310279	0.000000	4.404078	0.0	1.697537	0.14992
4	1338.55	olivine1	No	46.724667	0.724176	14.197934	1.457974	0.0	7.438616	0.163195	14.259132	10.781046	2.019534	0.142795	0.000000	2.090931	0.0	8.750491	0.14992

In this initial calculation, the liquidus temperature is calculated using the pMELTS model, which was designed for use on mantle-like bulk compositions at 1 - 3 GPa (Ghiorso et al. 2002). For the majorty of samples in this database it may be more appropriate to use the rhyoliteMELTSv1.0.2 (Gualda et al. 2012) or rhyoliteMELTSv1.2.0 (Ghiorso and Gualda, 2015) models. We can rerun the same calculations but change the Model selected:

Results_MELTSv102 = M.findLiq_multi(Model = "MELTSv1.0.2", # which MELTS model to use - "pMELTS", "MELTSv1.0.2", "MELTSv1.1.0", or "MELTSv1.2.0"
                        bulk = Test, # initial composition(s) wither a dictionary (if one composition) or a DataFrame (multiple compositions)
                        T_initial_C = Test['T_K'].values-200, # Initial guess for the liquidus temperature (optional)
                        P_bar = Test['P_kbar'].values*1000, # Pressure (in bars) of the calculation
                        Fe3Fet_Liq = 0.15) # initial Fe3Fet_Liq ratio of the sample(s)

Results_MELTSv120 = M.findLiq_multi(Model = "MELTSv1.2.0", # which MELTS model to use - "pMELTS", "MELTSv1.0.2", "MELTSv1.1.0", or "MELTSv1.2.0"
                        bulk = Test, # initial composition(s) wither a dictionary (if one composition) or a DataFrame (multiple compositions)
                        T_initial_C = Test['T_K'].values-200, # Initial guess for the liquidus temperature (optional)
                        P_bar = Test['P_kbar'].values*1000, # Pressure (in bars) of the calculation
                        Fe3Fet_Liq = 0.15) # initial Fe3Fet_Liq ratio of the sample(s)

Now that we have the results of this analysis, we can compare the results of these calculations to the true experimental temperature. Experiments that have no reported water typically show a very poor match to the experimental temperature, and are thus excluded from the following comparison. Additionally, we exlude any calculations that did not return a result (Results[‘T_C_liq’] = 0).

Overall, MELTS does an acceptable job in reproducing the experimental temperatures, although it is notable that the liquidus temperature is typically overpredicted by around 50 \(^{o}\)C.

f, a = plt.subplots(1,3, figsize = (12,4), sharex = True, sharey = True)

a[0].plot(Test['T_K'][(Test['H2O_Liq'] > 0) & (Results_pMELTS['T_Liq'] > 0)]-273.15,
       Results_pMELTS['T_Liq'][(Test['H2O_Liq'] > 0) & (Results_pMELTS['T_Liq'] > 0)],
        'ok')
a[0].plot([800,1400],[800,1400], 'r-')
a[0].set_ylabel('Estimated Temperature ($^{o}$C)')
a[0].set_xlabel('Experimental Temperature ($^{o}$C)')
a[0].text(750, 1460, 'pMELTS')


a[1].plot(Test['T_K'][(Test['H2O_Liq'] > 0) & (Results_MELTSv102['T_Liq'] > 0)]-273.15,
       Results_MELTSv102['T_Liq'][(Test['H2O_Liq'] > 0) & (Results_MELTSv102['T_Liq'] > 0)],
        'ok')
a[1].plot([800,1400],[800,1400], 'r-')
a[1].set_xlabel('Experimental Temperature ($^{o}$C)')
a[1].text(750, 1460, 'rhyoliteMELTSv1.0.2')


a[2].plot(Test['T_K'][(Test['H2O_Liq'] > 0) & (Results_MELTSv120['T_Liq'] > 0)]-273.15,
       Results_MELTSv120['T_Liq'][(Test['H2O_Liq'] > 0) & (Results_MELTSv120['T_Liq'] > 0)],
        'ok')
a[2].plot([800,1400],[800,1400], 'r-')
a[2].set_xlabel('Experimental Temperature ($^{o}$C)')
a[2].text(750, 1460, 'rhyoliteMELTSv1.2.0')

Text(750, 1460, 'rhyoliteMELTSv1.2.0')

findfont: Font family ['cmsy10'] not found. Falling back to DejaVu Sans.
findfont: Font family ['cmr10'] not found. Falling back to DejaVu Sans.
findfont: Font family ['cmtt10'] not found. Falling back to DejaVu Sans.
findfont: Font family ['cmmi10'] not found. Falling back to DejaVu Sans.
findfont: Font family ['cmb10'] not found. Falling back to DejaVu Sans.
findfont: Font family ['cmss10'] not found. Falling back to DejaVu Sans.
findfont: Font family ['cmex10'] not found. Falling back to DejaVu Sans.