Liquid-only thermometry¶

This script shows how to use the various options for liquid-only thermometry
This tutorial is also a great starting point when using other thermobarometers, as the choice of inputs is the same as for other phases
You can download the Excel spreadsheet you need here: https://github.com/PennyWieser/Thermobar/blob/main/docs/Examples/Liquid_Ol_Liq_Themometry/Liquid_only_Thermometry.xlsx

You need to install Thermobar once on your machine, if you haven't done this yet, uncomment the line below (remove the #)¶

In [1]:

#!pip install Thermobar

In [2]:

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import Thermobar as pt

Load in data¶

The function "import_excel" allows you to specify the name of the excel file and the sheet name.
Data should have the headings SiO2_Liq, TiO2_Liq etc. etc. The order of headings doesn't matter.
As in the 5 min intro example, you could also have no _Liq in the headings, and load them in the functoin using suffix="_Liq"
You can also have any other columns, e.g., estimate of pressure from any other proxy (melt inclusions, geophysics), and anything else you might want to plot (e.g., latitude, longitude)

In [3]:

out=pt.import_excel('Liquid_only_Thermometry.xlsx', sheet_name="Liquid_only")
my_input=out['my_input']
myLiquids1=out['Liqs']

Inspect inputs to check you didn't have any funny column headings.¶

It is always useful to inspect the outputs from import_excel, sometimes your column headings may have funny characters due to use of spaces, subscripts etc. in journal pdf tables. Check that all the columns you entered have numbers. If, say your SiO2_Liq heading had funny characters, this column will be filled with zeros when you inspect it.

In [4]:

myLiquids1.head()

Out[4]:

	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	Sample_ID_Liq
0	57.023602	0.623106	16.332899	4.36174	0.103851	4.19180	6.94858	3.59702	0.896895	0.000000	0.226584	5.59	0.2	0
1	57.658600	0.654150	17.194799	3.90621	0.084105	2.86892	5.91538	3.85948	1.018600	0.000000	0.214935	6.55	0.2	1
2	60.731201	0.862054	17.144199	4.07781	0.077488	2.50867	5.22075	4.45556	1.414160	0.000000	0.319638	3.14	0.2	2
3	61.532799	0.440860	16.508801	3.32990	0.037520	1.64150	4.34294	4.40860	1.407000	0.000000	0.215740	6.20	0.2	3
4	52.969101	0.803412	17.563000	5.93217	0.149472	3.78351	7.65110	3.80219	0.551178	0.037368	0.196182	6.58	0.2	4

If at any point you need help with what functions inputs to use, you can type help(pt.function_name)¶

In [5]:

# Options for different thermometers laid out in the "equationT" =
help(pt.calculate_liq_only_temp)

Help on function calculate_liq_only_temp in module Thermobar.liquid_thermometers:

calculate_liq_only_temp(*, liq_comps, equationT, P=None, H2O_Liq=None, print=False)
     Liquid-only thermometery. Returns a temperature in Kelvin.
    
    Parameters
     -------
    
     liq_comps: pandas.DataFrame
         liquid compositions with column headings SiO2_Liq, MgO_Liq etc.
    
     equationT: str
         If has _sat at the end, represents the saturation surface of that mineral.
    
         Equations from Putirka et al. (2016).
             | T_Put2016_eq3_amp_sat (saturation surface of amphibole)
    
         Equations from Putirka (2008) and older studies:
    
             | T_Put2008_eq13
             | T_Put2008_eq14
             | T_Put2008_eq15
             | T_Put2008_eq16
             | T_Put2008_eq34_cpx_sat
             | T_Put2008_eq28b_opx_sat
             | T_Put1999_cpx_sat
             * Following 3 thermometers are adaptations of olivine-liquid thermometers with  DMg calculated using Beattie 1993,
             This means you can use them without knowing an olivine composition. ocan be applied when you haven't measured an olivine composiiton.
             | T_Put2008_eq19_BeattDMg
             | T_Put2008_eq21_BeattDMg
             | T_Put2008_eq22_BeattDMg
    
         Equations from Sugawara (2000):
    
             | T_Sug2000_eq1
             | T_Sug2000_eq3_ol
             | T_Sug2000_eq3_opx
             | T_Sug2000_eq3_cpx
             | T_Sug2000_eq3_pig
             | T_Sug2000_eq6a
             | T_Sug2000_eq6b
    
         Equations from Helz and Thornber (1987):
    
             | T_Helz1987_MgO
             | T_Helz1987_CaO
    
         Equation from Molina et al. (2015)
    
             | T_Molina2015_amp_sat
    
         Equation from Montrieth 1995
            | T_Montierth1995_MgO
    
         Equation from Beattie (1993)
            | T_Beatt1993_opx
    
     P: float, int, pandas.Series, str  ("Solve")
         Pressure in kbar
         Only needed for P-sensitive thermometers.
         If enter P="Solve", returns a partial function
         Else, enter an integer, float, or panda series
    
     H2O_Liq: optional.
         If None, uses H2O_Liq column from input.
         If int, float, pandas.Series, uses this instead of H2O_Liq Column
    
    
     Returns
     -------
     pandas series
        Temperature in K

Example 1 - P and H2O-independent thermometers¶

The thermometer of Helz and Thornber 1987 is not dependent on pressure or H2O.
Thus, the function only requires 2 inputs, the name of pandas dataframe with liquid compositions (here, myLiquids1) and the name of the equationT you have choosen
All functions output in Kelvin and Kbar, so to get celcius, we subtract 273.15

In [6]:

T_Helz1987=pt.calculate_liq_only_temp(liq_comps=myLiquids1,  equationT="T_Helz1987_MgO")-273.15
# In jupyter lab/notebooks, outputs don't automatically display, 
# but it prints the last line in each cell
T_Helz1987

Out[6]:

0    1098.255182
1    1071.665294
2    1064.424269
3    1046.994150
4    1090.048550
5    1079.746499
6    1067.772325
7    1099.844089
8    1077.533086
9    1121.514497
Name: MgO_Liq, dtype: float64

Example 2 - Pressure-dependent thermometer¶

a) If you select an equation which is Pressure-dependent, and don't specify a pressure, the code returns an error

b) You can either select P="Solve" which returns a partial function. This means you can evaluate it at any pressure you want easily

c) Or you can specify a fixed pressure, e.g., P=5 (which runs all calculations at 5 kbar)

d) Or you can specify pressure based on a column in your original spreadsheet containing pressure

Example 2a - Don't specify a pressure¶

returns an error, telling you to input a pressure

In [7]:

Teq15_2kbar_6wt=pt.calculate_liq_only_temp(liq_comps=myLiquids1,  
                equationT="T_Put2008_eq15")-273.15 # Convert to Celcius
Teq15_2kbar_6wt

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Input In [7], in <cell line: 1>()
----> 1 Teq15_2kbar_6wt=pt.calculate_liq_only_temp(liq_comps=myLiquids1,  
      2                 equationT="T_Put2008_eq15")-273.15 # Convert to Celcius
      3 Teq15_2kbar_6wt

File g:\my drive\postdoc\pymme\mybarometers\thermobar_outer\src\Thermobar\liquid_thermometers.py:526, in calculate_liq_only_temp(liq_comps, equationT, P, H2O_Liq, print)
    524 if sig.parameters['P'].default is not None:
    525     if P is None:
--> 526         raise ValueError(f'{equationT} requires you to enter P, or specify P="Solve"')
    528 #else:
    529     # print('gothere2')
    530     # if P is not None:
   (...)
    533     #     P=None
    534     #     print('got here')
    538 kwargs = {name: anhyd_cat_frac[name] for name, p in sig.parameters.items() if p.kind == inspect.Parameter.KEYWORD_ONLY}

ValueError: T_Put2008_eq15 requires you to enter P, or specify P="Solve"

Example 2b - specify P="Solve" to return a partial function¶

In [8]:

Teq15_partial=pt.calculate_liq_only_temp(liq_comps=myLiquids1,  
            equationT="T_Put2008_eq15", P="Solve")# Output is in Kelvin

In [9]:

# You can then evaluate this partial at any pressure you want (here, 10 kbar).
# This can save computation time
Teq15_10kbar=Teq15_partial(10)-273.15
Teq15_10kbar

Out[9]:

0    1082.230315
1    1031.095340
2    1065.719538
3     991.644728
4    1049.453769
5    1033.581813
6    1021.013571
7     946.501846
8     919.539200
9    1108.248454
dtype: float64

Example 2c - evaluate at constant pressure¶

In [10]:

# Here P=5 kbar
Teq15_5kbar=pt.calculate_liq_only_temp(liq_comps=myLiquids1,  
            equationT="T_Put2008_eq15", P=5)-273.15 # Output is in Kelvin, convert to C
Teq15_5kbar 

Out[10]:

0    1062.650315
1    1011.515340
2    1046.139538
3     972.064728
4    1029.873769
5    1014.001813
6    1001.433571
7     926.921846
8     899.959200
9    1088.668454
dtype: float64

Example 2d - Evaluate at values given by P_kbar column in inputted excel spreadsheet¶

your column can be called anything, you just need state P=my_input['column name']. If you stored pressure in GPa say, you could do P=10*my_input['P_GPa']

In [11]:

Teq15_input=pt.calculate_liq_only_temp(liq_comps=myLiquids1,  
            equationT="T_Put2008_eq15", P=my_input['P_kbar'])-273.15 
# Output is in Kelvin, convert to c.
Teq15_input

Out[11]:

0    1050.902315
1    1003.683340
2    1046.139538
3     972.064728
4    1029.873769
5    1014.001813
6    1001.433571
7     926.921846
8     899.959200
9    1088.668454
dtype: float64

Example 3 - Specifying water content¶

The dataframe for liquids in this example has a H2O column it read from the excel spreadsheet. By default, H2O-dependent thermometers will use the values entered in this column
If you don't have this column, it'll be filled with zeros, so calculations will be done at H2O=0
However, you can also override this by specifying H2O_Liq=value in the function.

In [12]:

Teq15_2H2O=pt.calculate_liq_only_temp(liq_comps=myLiquids1,  
        equationT="T_Put2008_eq15", P=5, H2O_Liq=2)-273.15 # Output is in Kelvin
Teq15_2H2O

Out[12]:

0    1108.710017
1    1069.891842
2    1060.765739
3    1025.950726
4    1088.635168
5    1075.714113
6    1058.655371
7    1083.961043
8    1058.281402
9    1129.211252
dtype: float64

We can also investigate how sensitive our calculations are to water, say we think maybe water ranges from 2-4 wt%, we can do the same calculation at 4 wt%¶

In [13]:

Teq15_4H2O=pt.calculate_liq_only_temp(liq_comps=myLiquids1,  
                                      equationT="T_Put2008_eq15", P=5, H2O_Liq=4)-273.15 # Output is in Kelvin
Teq15_4H2O

Out[13]:

0    1083.050017
1    1044.231842
2    1035.105739
3    1000.290726
4    1062.975168
5    1050.054113
6    1032.995371
7    1058.301043
8    1032.621402
9    1103.551252
dtype: float64

We can now calculate the difference between these by subtracting one from the other. Here, as H2O_Liq is just a constant on the function, the difference is constant for all liquids¶

In [14]:

Teq15_4H2O-Teq15_2H2O

Out[14]:

0   -25.66
1   -25.66
2   -25.66
3   -25.66
4   -25.66
5   -25.66
6   -25.66
7   -25.66
8   -25.66
9   -25.66
dtype: float64

Exporting data¶

There are a number of ways you can export data to excel
The first example appends these new calculations onto a copy of the "my_input" dataframe. This means you can export a spreadsheet that looks exactly like your input spreadsheet, but with the values for the new calculations

In [15]:

# First make a copy of the dataframe so you don't overwrite the original
my_input_c=my_input.copy()
# To make a new column, you specify my_input_c['new column name']= value
my_input_c['T_Helz1987']=T_Helz1987
# You can add as many new columns as you wish
my_input_c['Teq15_4H2O']=Teq15_4H2O
# You can give your columns in the new dataframe any name you wish...
my_input_c['Temp eq 15 4 wt% H2O']=Teq15_4H2O
# Now, you ue the "to_excel" panda function to save this new appended dataframe to excel. The name of the dataframe you've been appended to goes on the let, and the name of the
# file you want to make goes bewteen the ' ' signs
my_input_c.to_excel('Thermometry_out1.xlsx')

In [ ]: