This page was generated from docs/Examples/Integration_with_VESIcal/Combining_VESIcal_Thermobar_SatPs.ipynb. Interactive online version: .

Combining Thermobar with VESIcal for Sat P at known T

We show how to use various thermometers in Thermobar to calculate temperature, which can then be fed into saturation pressure calculations
We also show how to combine thermometry (which is often sensitive to water content) with VESIcal calculate dissolved volatile calculations for insights into different H\(_2\)O contents at a range of crustal depths
You can find the excel spreadsheet here: https://github.com/PennyWieser/Thermobar/blob/main/docs/Examples/Integration_with_VESIcal/Ol_hosted_melt_inclusions.xlsx

!!!!! Please Note !!!!!!!

VESIcal requires Thermoengine to run. If you do not have Thermoengine installed locally on your laptop, you will either need to instal it, or more simply to start with, run calculations using the ENKI server (http://enki-portal.org/). For more information, watch the following YouTube videos:

https://www.youtube.com/watch?v=BuwUhu9YdG4&t=2s (Enki server)
https://www.youtube.com/watch?v=FRpsDbouuec&t=763s (Worked example)

First, load the python things you might need

[ ]:

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

Now, install VESIcal if you aren’t running on the server

[ ]:

#!pip install VESIcal

Now install Thermobar if you dont already have it

[ ]:

#!pip install Thermobar

Now import both Thermobar and VESIcal

[ ]:

# You need to uncomment this if you are running on the enki server
import VESIcal as v
print('VESIcal version')
print(v.__version__)

[ ]:

import Thermobar as pt
# This prints the versions
print('Thermobar version')
print(pt.__version__)

Example 1 - Simple integration of thermometry and saturation pressures

Here, we calculating saturation pressures for olivine-hosted melt inclusions from Kilauea using temperatures from a specific thermometer

Step 1: import the data

Here, we import the data using Thermobars input structure, e.g., oxides for a liquid are followed by _Liq, and oxides for olivine are followed by _Ol

[ ]:

out=pt.import_excel('Ol_hosted_melt_inclusions.xlsx', sheet_name="Ol-Liq")

# This subdivdes outputs into a dataframe for all inputs (my_input), ols, and liqs
my_input=out['my_input']
myOls=out['Ols']
myLiquids1=out['Liqs']

## Lets check the outputs have loaded right
display(myOls.head())
display(myLiquids1.head())

Step 2: lets calculate the temperature using the Helz and Thornber (1987) thermometer which only uses liquid compositions.

[ ]:

T_HT87=pt.calculate_liq_only_temp(liq_comps=myLiquids1, equationT="T_Helz1987_MgO")
T_HT87

Step 3: lets calculate the saturation pressure using VESIcal at this temperature.

The function “convert_To_VESIcal” converts a liquid dataframe into the form needed for VESIcal, and adds a new column called “Temperature” in celcius based on an input you have given it in Kelvin from Thermobar
It strips away the _Liq from each oxide used by Thermobar, and partitions into FeO and Fe2O3 depending user-entered Fe3Fet_Liq, as some solubility models are sensitive to redox.

[ ]:

# This function
df_HT87=pt.convert_to_vesical(liq_comps=myLiquids1,
                              T1=T_HT87, unit='Kelvin', Fe3Fet_Liq=0.15)
df_HT87.head()

Step 4: Now we do the final step to convert this dataframe into an object that we can enter into the calculation structures of VESIcal

[ ]:

VESIcal_out=v.BatchFile(filename=None, dataframe=df_HT87, label='Sample_ID')
VESIcal_out.data

Step 5: Final step - Calculate saturation pressure!

[ ]:

# First lets use the default model MagmaSat
SatP_HT87=VESIcal_out.calculate_saturation_pressure(temperature="Temp")
SatP_HT87.head()

Example 2 - Using a P-sensitive thermometer

The first example was relatively simple, because the thermometer we choose wasn’t sensitive to pressure
Many more recent thermometers also have a pressure term. For example, the Sugawara (2000) eq 3 olivine thermometer requires a pressure
Thus, we show a way to iterate towards a solution
First, we use the saturation pressure calculated above using the Helz and Thornber thermometer as a first guess. We then use the Sugawara thermometer as a P-sensitive thermometer to calculate a new temperature based on this pressure.

Step 1 - Calc a new temp using a P-sensitive thermometer

[ ]:

# Here, we get temperature from the sat Ps calculated above, remember to divide by 1000 as Thermobar wants kbar not bar.
T_Sug=pt.calculate_liq_only_temp(liq_comps=myLiquids1, equationT="T_Sug2000_eq3_ol",
                                 P=SatP_HT87['SaturationP_bars_VESIcal'].values/1000)

Step 2 - Lets add a new column called Temp_Sug and store our temp calculated (converted to values and to float)

[ ]:

VESIcal_out.data["Temp_Sug"] =T_Sug.values.astype(float)-273.15 # Remmeber to store as celcius
VESIcal_out.data

Step 3: Now we calculate saturation pressure using this new temperature

[ ]:

# Remember to subtract 273.15 if using direct outputs from thermobar
# (which use Kelvin, VESical wants celcius)
SatP_Sug1=VESIcal_out.calculate_saturation_pressure(temperature="Temp_Sug")

Example 3 - Using an olivine-Liquid thermometer which is sensitive to temperature and pressure

Here, we use the solubility model of IaconoMarziano

Step 1: Calculate temperature using a best guess of your pressure (e.g from knowledge of your system, 10 kbar here)

[ ]:

T_Put_10kbar=pt.calculate_ol_liq_temp(ol_comps=myOls, liq_comps=myLiquids1,
                                equationT="T_Put2008_eq22", P=10)

Step 3: Append a new temperature onto the dataframe

[ ]:

# Now we append this as a third temperature. We need T_K_calc
# because this thermometer returns a dataframe
# Remmeber to store as celcius
VESIcal_out.data["Temp_Put_10kbar"] =T_Put_10kbar['T_K_calc'].values.astype(float)-273.15
VESIcal_out.data

Step 4: Calculate sat P using this temperature

[ ]:

# Now we calculate saturation pressure using this temperature
SatP_Put10kbar=VESIcal_out.calculate_saturation_pressure(temperature="Temp_Put_10kbar",
                                model="IaconoMarziano")

[ ]:

# We can see our guess of 10 kbar was way too high
plt.hist(SatP_Put10kbar['SaturationP_bars_VESIcal']/1000)
plt.plot([10, 10], [0, 10], '-r')

Step 5: We could then repeat this process to basically “iterate” down to the right pressure

[ ]:

# Here we calculate a new temprature using these newly calculated pressures
T_Put_PSat1=pt.calculate_ol_liq_temp(ol_comps=myOls, liq_comps=myLiquids1,
equationT="T_Put2008_eq22", P=SatP_Put10kbar['SaturationP_bars_VESIcal'].values/1000)

# Add another new column with this new temperature
# Remmeber to store as celcius
VESIcal_out.data["Temp_Put_PSat1"] =T_Put_PSat1['T_K_calc'].values.astype(float)-273.15

# Calculate saturation pressure again using this new temperature.
SatP_Put_PSat1=VESIcal_out.calculate_saturation_pressure(temperature="Temp_Put_PSat1",
                model="IaconoMarziano")

Step 6 - Lets compare the pressures calculated vs. the first step where we assumed 10 kbar.

[ ]:

plt.plot(VESIcal_out.data['Temp_Put_PSat1'] - VESIcal_out.data['Temp_Put_10kbar'],
    SatP_Put_PSat1['SaturationP_bars_VESIcal']-SatP_Put10kbar['SaturationP_bars_VESIcal'], 'ok')
plt.xlabel('Difference in T \n (P=10 kbar - P calc)')
plt.ylabel('Difference in Pressure \n for different T (bars)')

Example 4 - Lets generally compare thermometer models and how much difference they make

[ ]:

T_Put22_DMg=pt.calculate_liq_only_temp(liq_comps=myLiquids1,
            equationT="T_Put2008_eq22_BeattDMg", P=2)

# Remmeber to store as celcius
VESIcal_out.data["T_Put22_DMg"] =T_Put22_DMg.values.astype(float)-273.15

SatP_T_Put22_DMg=VESIcal_out.calculate_saturation_pressure(temperature="T_Put22_DMg")

[ ]:

plt.plot(SatP_HT87['H2O'],
    SatP_T_Put22_DMg['SaturationP_bars_VESIcal']-SatP_HT87['SaturationP_bars_VESIcal'], 'ok', mfc='red')
plt.xlabel('H$_2$O content',  fontsize=13)
plt.ylabel('Difference in Pressure: \n HT87 \n vs. P2008 eq22 (bars)', fontsize=13)

Summary

This notebook shows that VESIcal combined with Thermobar offers you lots of options. In reality, you will find that basaltic solubility models really aren’t that sensitive to temperature, so in general there is no need to iterate tempeature and pressure (but you can confirm this yourself)
But, in water-rich arcs, choosing a temperature from say Helz and Thornber vs. a thermometer which has H2O included could make more of a difference.