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 and install the alphaMELTS for Python files. This can be easily achieved with one simple function, please see the installation guide available on ReadTheDocs.

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

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import petthermotools as ptt
from tqdm.notebook import tqdm, trange
# import pickle

# 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'path_to_alphaMELTS_files')

alphaMELTS for Python files successfully located.

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

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 (Wieser et al. 2025).

Data = pd.read_excel('MELTS_MAGEMIN_Liquidus.xlsx')
Data = Data[~np.isinf(Data["H2O_Liq"])]
Data = Data.reset_index(drop=True)

The entire experimental dataset contains over 2000 samples. Therefore, to reduce the computational time in this example 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'])), 200, replace=False)
Test = Data.copy()
Test = Test.loc[Ch]

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

new_comp = ptt.comp_fix(comp=Test)
new_comp

	SiO2_Liq	TiO2_Liq	Al2O3_Liq	FeOt_Liq	MnO_Liq	MgO_Liq	CaO_Liq	Na2O_Liq	K2O_Liq	Cr2O3_Liq	P2O5_Liq	H2O_Liq	Fe3Fet_Liq	CO2_Liq
0	64.066798	0.874263	13.831041	6.218075	0.147348	1.640472	4.518664	5.157171	1.227898	0.000000	0.255403	2.062868	0.15	0.0
1	63.000000	0.500000	16.400000	3.240000	0.100000	1.460000	4.180000	4.890000	1.010000	0.000000	0.000000	5.220000	0.15	0.0
2	49.871506	1.746048	14.087436	12.420753	0.218256	6.726254	11.646936	2.271847	0.198415	0.019841	0.158732	0.633976	0.15	0.0
3	51.976977	1.561151	17.509859	8.666432	0.146426	4.394229	6.773000	3.322450	0.992945	0.000000	0.253837	4.402694	0.15	0.0
4	46.958829	0.583149	16.113323	7.560473	0.112537	10.486448	9.688456	2.772515	0.358074	0.122768	0.214844	5.028583	0.15	0.0
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
195	65.164083	0.242702	14.659290	2.293986	0.147794	0.526606	2.088591	3.429156	3.027627	0.000000	0.220769	8.199396	0.15	0.0
196	66.200000	0.580000	14.900000	3.190000	0.060000	0.890000	3.110000	5.050000	1.380000	0.000000	0.000000	4.640000	0.15	0.0
197	69.566753	0.186007	12.648501	0.985839	0.102304	0.437117	1.320652	3.124924	4.352572	0.000000	0.167407	7.107924	0.15	0.0
198	60.503949	0.149985	22.157784	1.359864	0.319968	0.099990	0.829917	8.629137	5.949405	0.000000	0.000000	0.000000	0.15	0.0
199	67.227273	0.345455	12.945455	2.554545	0.063636	0.390909	2.872727	2.100000	2.300000	0.000000	0.109091	9.090909	0.15	0.0

200 rows × 14 columns

T = ptt.estimate_t(comp = new_comp, P_bar = Test['P_kbar']*1000)

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 = ptt.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)
                        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 reveals the information that is recorded by this function. 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_C	liquidus_phase	fluid_saturated	SiO2_Liq	TiO2_Liq	Al2O3_Liq	Cr2O3_Liq	FeOt_Liq	MnO_Liq	MgO_Liq	CaO_Liq	Na2O_Liq	K2O_Liq	P2O5_Liq	H2O_Liq	CO2_Liq	Fe3Fet_Liq
0	1046.589636	plagioclase1	No	64.067405	0.874305	13.830595	0.000000	6.217798	0.147355	1.640551	4.518543	5.157135	1.227931	0.255415	2.062968	0.0	0.149921
1	963.322929	plagioclase1	No	63.002055	0.500139	16.396922	0.000000	3.240599	0.100028	1.460405	4.178467	4.889776	1.010159	0.000000	5.221450	0.0	0.149921
2	1196.603712	olivine1	No	49.872499	1.746130	14.088097	0.019842	12.419525	0.218056	6.725265	11.647463	2.271954	0.198424	0.158739	0.634005	0.0	0.149929
3	1059.622056	olivine1	Yes	53.089794	0.327077	17.999068	0.000000	8.737000	0.130302	4.170554	6.959927	3.415276	1.020686	0.260929	3.889386	0.0	0.152851
4	1229.343347	olivine1	No	46.959433	0.583179	16.114143	0.122774	7.559720	0.112519	10.484852	9.688938	2.772656	0.358092	0.214855	5.028839	0.0	0.149930

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 = ptt.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)
                        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_MELTSv102.head()

	T_Liq_C	liquidus_phase	fluid_saturated	SiO2_Liq	TiO2_Liq	Al2O3_Liq	Cr2O3_Liq	FeOt_Liq	MnO_Liq	MgO_Liq	CaO_Liq	Na2O_Liq	K2O_Liq	P2O5_Liq	H2O_Liq	CO2_Liq	Fe3Fet_Liq
0	1040.789636	orthopyroxene1	No	64.067685	0.874305	13.831670	0.000000	6.216893	0.147355	1.639473	4.518827	5.157438	1.227962	0.255416	2.062976	0.0	0.149938
1	958.422929	orthopyroxene1	No	63.000489	0.500015	16.400524	0.000000	3.239280	0.100004	1.459187	4.180103	4.890175	1.010036	0.000000	5.220188	0.0	0.149941
2	1202.103712	clinopyroxene1	No	49.871918	1.746284	14.089198	0.019845	12.420523	0.218296	6.724750	11.645662	2.272224	0.198450	0.158760	0.634090	0.0	0.149913
3	1095.522056	olivine1	Yes	52.516888	1.577381	17.691899	0.000000	8.755330	0.147945	4.439146	6.843411	3.356992	1.003268	0.256476	3.411264	0.0	0.149928
4	1306.743347	olivine1	No	46.959380	0.583174	16.114021	0.122773	7.559750	0.112540	10.485116	9.688868	2.772635	0.358089	0.214854	5.028801	0.0	0.149928

Results_MELTSv120 = ptt.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)
                        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.head()

	T_Liq_C	liquidus_phase	fluid_saturated	SiO2_Liq	TiO2_Liq	Al2O3_Liq	Cr2O3_Liq	FeOt_Liq	MnO_Liq	MgO_Liq	CaO_Liq	Na2O_Liq	K2O_Liq	P2O5_Liq	H2O_Liq	CO2_Liq	Fe3Fet_Liq
0	1041.189636	plagioclase1	No	64.067515	0.874309	13.830485	0.000000	6.217820	0.147355	1.640557	4.518456	5.157161	1.227950	0.255416	2.062976	0.0	0.149921
1	971.222929	plagioclase1	No	63.002159	0.500098	16.397399	0.000000	3.240332	0.100020	1.460285	4.178394	4.890115	1.010180	0.000000	5.221019	0.0	0.149921
2	1195.303712	clinopyroxene1	No	49.871992	1.746273	14.089123	0.019845	12.420486	0.218296	6.724819	11.645645	2.272221	0.198451	0.158760	0.634091	0.0	0.149911
3	1094.722056	spinel1	Yes	52.568804	1.578587	17.709042	0.000000	8.763007	0.148094	4.444079	6.850119	3.360281	1.004251	0.256727	3.317010	0.0	0.149877
4	1389.243347	spinel1	No	46.959313	0.583153	16.113383	0.122683	7.559774	0.112539	10.486498	9.688555	2.772543	0.358078	0.214847	5.028635	0.0	0.149917

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_C'] > 0)]-273.15,
       Results_pMELTS['T_Liq_C'][(Test['H2O_Liq'] > 0) & (Results_pMELTS['T_Liq_C'] > 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(800, 1460, 'pMELTS')


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


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

a[2].set_ylim([700,1500])

(700.0, 1500.0)