This page was generated from
docs/Examples/DecompressionTests/StHelens_decompress_simple.ipynb.
Interactive online version:
.
Isothermal decompression calculations
In addition to crystallization and melting calculations, thermodynamic models are often used to investigate the processes of magma decompression.
In this notebook we focus on the fluid composition during magma decompression by investigating different initial H2O contents.
Before any calculations can be run users need to download the alphaMELTS for Python files. Please see the installation guide on ReadTheDocs.
You will need to also download this excel spreadsheet: https://github.com/gleesonm1/PetThermoTools/blob/master/docs/Examples/DecompressionTests/StHelens.xlsx
[1]:
# import core python packages that we'll use for plotting and data manipulation.
import sys
import os
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
# 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'C:\Users\penny\Box\Berkeley_new\MELTS_Installation\alphamelts_py\alphamelts-py-2.3.1-win64')
import petthermotools as ptt
ptt.__version__
alphaMELTS for Python files successfully located.
[1]:
'0.3.21dev'
To suppress outputs from the MELTS calculations run the cell below twice.
[3]:
import platform
if platform.system() == "Darwin":
import sys
import os
sys.stdout = open(os.devnull, 'w')
sys.stderr = open(os.devnull, 'w')
In this example we’ll load melt inclusion compositions from the White Pumice of Mt St Helens published the study of Blundy and Cashman (2005).
[4]:
# Melt inclusion and groundmass glass data from Mt St Helens (Blundy and Cashman, 2005)
StHelens = pd.read_excel('StHelens.xlsx', sheet_name='1980_WhitePumice')
StHelens = StHelens.fillna(0.0)
StHelens.head()
[4]:
| sample | Type | SiO2 | TiO2 | Al2O3 | FeOt | MnO | MgO | CaO | Na2O | K2O | H2O | Total | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | KCHB-4Aa | A1 | 68.05 | 0.33 | 12.96 | 2.72 | 0.12 | 0.98 | 1.95 | 5.60 | 1.93 | 4.63 | 99.26 |
| 1 | KCHB-7Ab | A1 | 69.01 | 0.37 | 14.07 | 2.19 | 0.02 | 0.36 | 2.25 | 5.54 | 1.85 | 3.19 | 99.09 |
| 2 | KCHB-8Aa | A1 | 70.19 | 0.35 | 14.70 | 2.52 | 0.03 | 0.34 | 2.24 | 6.51 | 1.72 | 1.83 | 100.43 |
| 3 | KCHB-I IR | A1 | 65.81 | 0.51 | 15.03 | 2.13 | 0.00 | 0.37 | 2.35 | 5.65 | 2.44 | 4.92 | 99.2 |
| 4 | KCHB-I IB' | A1 | 69.71 | 0.21 | 14.15 | 2.48 | 0.00 | 0.45 | 2.23 | 6.39 | 2.03 | 2.66 | 100.29 |
For these simple decompression calculations we can simply pick any of the inclusions from the Blundy and Cashman (2005) study. Here we select a plagioclase-hosted melt incluion:
[5]:
# select a representative composition of the Mt St Helens white pumice. Here we've chosen a more mafic (~65 wt% SiO2) water rich inclusion.
WhitePumice = StHelens.loc[23]
WhitePumice
[5]:
sample 006-IOAO
Type PI
SiO2 64.91
TiO2 0.41
Al2O3 13.31
FeOt 2.78
MnO 0.02
MgO 0.71
CaO 1.95
Na2O 5.49
K2O 1.99
H2O 6.38
Total 98.27
Name: 23, dtype: object
The inclusion has a measured H2O content of 6.38 wt%, but it is worth investigating how the degassing behaviour - and volatile composition - might change with different assumed initial CO2 contents. In the following cell we set up a decompression calculation using the measured H2O content and 4 different initial CO2 contents (in wt%).
Fluid_sat = True means if there are excess volatiles at the start of the calculation, they are removed. If you want to start with an exsolved fluid phase, put this parameter to False
Right now this calculation is holding temp constant, in the more realistic case where you hold entropy constant, just swap out ptt.isothermal_decompression for ptt.isentropic_decompression
[6]:
# Lets run at 100, 500, 1000 ppm CO2
# This is closed system because we aren't fractionating fluids
Helens_decompress_closed = ptt.isothermal_decompression(Model="MELTSv1.2.0",
bulk=WhitePumice, find_liquidus=True, fO2_buffer="FMQ", fO2_offset=2,
P_start_bar=3000, P_end_bar=50, dp_bar=20, fluid_sat = True,
CO2_init = np.array([0.01, 0.05, 0.1]), Frac_fluid=False,
Frac_solid=False, label = 'CO2_init')
# This is open system for volatiles because we are fractionating fluids but not solids
Helens_decompress_open = ptt.isothermal_decompression(Model="MELTSv1.2.0",
bulk=WhitePumice, find_liquidus=True, fO2_buffer="FMQ", fO2_offset=2,
P_start_bar=3000, P_end_bar=50, dp_bar=20, fluid_sat = True,
CO2_init = np.array([0.01, 0.05, 0.1]), Frac_fluid=True,
Frac_solid=False, label = 'CO2_init')
## If you wanted isentropic instead
Helens_decompress_open_isentropic = ptt.isentropic_decompression(Model="MELTSv1.2.0",
bulk=WhitePumice, find_liquidus=True, fO2_buffer="FMQ", fO2_offset=2,
P_start_bar=3000, P_end_bar=50, dp_bar=20, fluid_sat = True,
CO2_init = np.array([0.01, 0.05, 0.1]), Frac_fluid=True,
Frac_solid=False, label = 'CO2_init')
We can inspect the outputs and see that each run is now labeled by the initial CO2 content of that model.
[7]:
Helens_decompress_closed.keys()
[7]:
dict_keys(['CO2 = 0.01 wt%', 'CO2 = 0.05 wt%', 'CO2 = 0.1 wt%'])
Examining the results from one model reveals the various Dictionaries and Pandas DataFrames that are created during the thermodynamic calculations.
If you’re following the YouTube example please note that there have been a few name changes from v0.2.55 onwards. Most importantly the keys ‘Mass’, ‘Volume’, and ‘rho’ in the output Dictionaries have been replaced by the keys ‘mass_g’, ‘volume_cm3’, and ‘rho_kg/m3’. In addition, nearly all outputs are now associated with their respective unit.
[8]:
Helens_decompress_closed['CO2 = 0.01 wt%'].keys()
[8]:
dict_keys(['Conditions', 'fluid1', 'fluid1_prop', 'liquid1', 'liquid1_prop', 'plagioclase1', 'plagioclase1_prop', 'spinel1', 'spinel1_prop', 'All', 'PhaseList', 'mass_g', 'volume_cm3', 'rho_kg/m3', 'Input'])
Specifically focusing on the fluid, we can use the cell below to examine the composition of the fluid phase for a single model.
[9]:
Helens_decompress_closed['CO2 = 0.05 wt%']['fluid1'].tail()
[9]:
| H2O_Fl | CO2_Fl | X_H2O_mol_Fl | X_CO2_mol_Fl | |
|---|---|---|---|---|
| 144 | 99.188994 | 0.811006 | 0.996666 | 0.003334 |
| 145 | 99.203554 | 0.796446 | 0.996726 | 0.003274 |
| 146 | 99.217925 | 0.782075 | 0.996786 | 0.003214 |
| 147 | 99.232278 | 0.767722 | 0.996845 | 0.003155 |
| 148 | 99.246917 | 0.753083 | 0.996905 | 0.003095 |
Lets look at fluid composition vs. pressure and H2O and CO2 in the melt remaining
We can further examine the results by plotting up the fluid composition against pressure. In the cell below we can loop through the main results Dictionary to access each model in turn.
[10]:
# set up figure layout
f, (ax1, ax2, ax3) = plt.subplots(1,3, figsize = (12, 5))
colors=['red', 'blue', 'k']
# loop through main Dictionary
for key, color in zip(Helens_decompress_closed, colors) :
# extract individial calculations
res_closed = Helens_decompress_closed[key]
res_open = Helens_decompress_open[key]
# plot fluid chemistry vs pressure
ax1.plot(res_closed['All']['X_H2O_mol_Fl'],
res_closed['All']['P_bar'],'-', color=color, label=key+' closed')
ax1.plot(res_open['All']['X_H2O_mol_Fl'],
res_open['All']['P_bar'],':',color=color,label=key+' open')
# plot melt CO2 vs pressure
ax2.plot(res_closed['All']['H2O_Liq'],
res_closed['All']['CO2_Liq'],'-',color=color,label=key)
ax2.plot(res_open['All']['H2O_Liq'],
res_open['All']['CO2_Liq'],':',color=color,label=key)
# plot melt CO2 vs pressure
ax3.plot(res_closed['All']['SiO2_Liq'],
res_closed['All']['P_bar'],'-',color=color,label=key)
ax3.plot(res_open['All']['SiO2_Liq'],
res_open['All']['P_bar'],':',color=color,label=key)
ax1.legend(fontsize = 8)
ax1.set_xlabel('X$_{H2O}$')
ax1.set_ylabel('Pressure (bar)')
ax2.set_xlabel('H$_2$O Liq')
ax2.set_ylabel('CO$_2$ (Liq)')
ax3.set_ylabel('Pressure (bar)')
ax3.set_xlabel('SiO$_2$ (Liq)')
ax1.legend()
f.tight_layout()
Why is SiO2 changing?
[11]:
# Lets look at what phases are forming
Helens_decompress_closed['CO2 = 0.05 wt%']['mass_g']
[11]:
| fluid1 | liquid1 | plagioclase1 | spinel1 | |
|---|---|---|---|---|
| 0 | 0.000224 | 100.095381 | 0.000000 | 0.000000 |
| 1 | 0.003674 | 100.091945 | 0.000000 | 0.000000 |
| 2 | 0.007216 | 100.088418 | 0.000000 | 0.000000 |
| 3 | 0.010855 | 100.084794 | 0.000000 | 0.000000 |
| 4 | 0.014598 | 100.081066 | 0.000000 | 0.000000 |
| ... | ... | ... | ... | ... |
| 144 | 5.634193 | 76.960340 | 16.568217 | 0.962348 |
| 145 | 5.737642 | 74.846723 | 18.495453 | 1.047435 |
| 146 | 5.843503 | 72.635386 | 20.514815 | 1.135787 |
| 147 | 5.953161 | 70.299081 | 22.651231 | 1.228362 |
| 148 | 6.069271 | 67.782636 | 24.955273 | 1.327157 |
149 rows × 4 columns
Lets plot what the phases are doing
[12]:
res_closed['All'].loc[:,res_closed['All'].columns.str.contains('mass')]
[12]:
| mass_g | mass_g_Fl | mass_g_Liq | mass_g_Plag | mass_g_Sp | |
|---|---|---|---|---|---|
| 0 | 99.958509 | 0.000226 | 99.958282 | NaN | NaN |
| 1 | 99.958523 | 0.003509 | 99.955014 | NaN | NaN |
| 2 | 99.958538 | 0.006872 | 99.951666 | NaN | NaN |
| 3 | 99.958553 | 0.010321 | 99.948232 | NaN | NaN |
| 4 | 99.958568 | 0.013860 | 99.944708 | NaN | NaN |
| ... | ... | ... | ... | ... | ... |
| 144 | 99.987757 | 5.543447 | 77.015964 | 16.474114 | 0.954231 |
| 145 | 99.989920 | 5.646954 | 74.899383 | 18.403999 | 1.039584 |
| 146 | 99.992165 | 5.752877 | 72.684512 | 20.426546 | 1.128229 |
| 147 | 99.994517 | 5.862599 | 70.344089 | 22.566702 | 1.221126 |
| 148 | 99.997029 | 5.978775 | 67.822881 | 24.875096 | 1.320277 |
149 rows × 5 columns
[13]:
# set up figure layout
f, (ax1, ax2, ax3) = plt.subplots(1,3, figsize = (12, 5))
# loop through main Dictionary
for key, color in zip(Helens_decompress_closed, colors) :
# extract individial calculations
res_closed = Helens_decompress_closed[key]
res_open = Helens_decompress_open[key]
# plot fluid chemistry vs pressure
ax1.plot(100*res_closed['All']['mass_g_Fl']/res_closed['All']['mass_g'],
res_closed['All']['P_bar'],'-', color=color, label=key + ' closed')
ax1.plot(100*res_open['All']['mass_g_Fl']/res_open['All']['mass_g'],
res_open['All']['P_bar'],':', color=color, label=key + ' open')
# plot melt CO2 vs pressure
ax2.plot(200*res_closed['All']['mass_g_Plag']/res_closed['All']['mass_g'],
res_closed['All']['P_bar'],'-', color=color, label=key)
ax2.plot(200*res_open['All']['mass_g_Plag']/res_open['All']['mass_g'],
res_open['All']['P_bar'],':', color=color, label=key)
# plot melt CO2 vs pressure
ax3.plot(100*res_closed['All']['mass_g_Sp']/res_closed['All']['mass_g'],
res_closed['All']['P_bar'],'-',label=key)
ax3.plot(100*res_open['All']['mass_g_Sp']/res_open['All']['mass_g'],
res_open['All']['P_bar'],':',label=key)
# ax1.legend(fontsize = 8)
ax1.set_xlabel('Fluid mass %')
ax1.set_ylabel('Pressure (bar)')
ax2.set_xlabel('Plag mass %')
ax2.set_ylabel('Pressure (bar)')
ax3.set_xlabel('Spinel mass %')
ax3.set_ylabel('Pressure (bar)')
ax1.legend()
f.tight_layout()
We can do the same for phase volume
[14]:
# set up figure layout
f, (ax1, ax2, ax3) = plt.subplots(1,3, figsize = (12, 5))
# loop through main Dictionary
for key, color in zip(Helens_decompress_closed, colors) :
# extract individial calculations
res_closed = Helens_decompress_closed[key]
res_open = Helens_decompress_open[key]
# plot fluid chemistry vs pressure
ax1.plot(100*res_closed['All']['V_cm^3_Fl']/res_closed['All']['V_cm^3'],
res_closed['All']['P_bar'],'-', color=color, label=key + ' closed')
ax1.plot(100*res_open['All']['V_cm^3_Fl']/res_open['All']['V_cm^3'],
res_open['All']['P_bar'],':', color=color, label=key + ' open')
# plot melt CO2 vs pressure
ax2.plot(200*res_closed['All']['V_cm^3_Plag']/res_closed['All']['V_cm^3'],
res_closed['All']['P_bar'],'-', color=color, label=key)
ax2.plot(200*res_open['All']['V_cm^3_Plag']/res_open['All']['V_cm^3'],
res_open['All']['P_bar'],':', color=color, label=key)
# plot melt CO2 vs pressure
ax3.plot(100*res_closed['All']['V_cm^3_Sp']/res_closed['All']['V_cm^3'],
res_closed['All']['P_bar'],'-',label=key)
ax3.plot(100*res_open['All']['V_cm^3_Sp']/res_open['All']['V_cm^3'],
res_open['All']['P_bar'],':',label=key)
ax1.set_xlabel('Fluid vol %')
ax1.set_ylabel('Pressure (bar)')
ax2.set_xlabel('Plag vol %')
ax2.set_ylabel('Pressure (bar)')
ax3.set_xlabel('Spinel vol %')
ax3.set_ylabel('Pressure (bar)')
ax1.legend()
f.tight_layout()
Lets think about what happens to liquid viscosity as we ascend
We are going to use the python package Thermobar to do this
[15]:
#! pip install Thermobar
[16]:
import Thermobar as pt
pt.__version__
[16]:
'1.0.67'
[17]:
# Example of how to calculate viscosity
Viscosity_500ppmCO2=pt.calculate_viscosity_giordano_2008(liq_comps=Helens_decompress_closed['CO2 = 0.05 wt%']['liquid1'],
T_K=Helens_decompress_closed['CO2 = 0.05 wt%']['All']['T_C']+273.15,
H2O_Liq=Helens_decompress_closed['CO2 = 0.05 wt%']['liquid1']['H2O_Liq'])
Viscosity_500ppmCO2.head()
[17]:
| n_melt | logn_melt | T_K | A | B | C | SiO2_Liq | TiO2_Liq | Al2O3_Liq | Cr2O3_Liq | ... | b13 | c1 | c2 | c3 | c4 | c5 | c6 | c11 | T_K | logn_melt | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 268.281673 | 2.428591 | 1274.19842 | -4.55 | 8513.383031 | 54.269780 | 66.172254 | 0.417973 | 13.568829 | 0.0 | ... | 805.151471 | 167.816457 | 120.474977 | 26.597082 | 20.035407 | -76.160248 | -302.833877 | 98.339981 | 1274.19842 | 2.428591 |
| 1 | 268.643049 | 2.429176 | 1274.19842 | -4.55 | 8513.866243 | 54.302730 | 66.174526 | 0.417987 | 13.569295 | 0.0 | ... | 805.266902 | 167.828486 | 120.483613 | 26.598989 | 20.036843 | -76.165707 | -302.806617 | 98.327123 | 1274.19842 | 2.429176 |
| 2 | 269.019899 | 2.429784 | 1274.19842 | -4.55 | 8514.369415 | 54.337043 | 66.176858 | 0.418002 | 13.569773 | 0.0 | ... | 805.387100 | 167.841011 | 120.492604 | 26.600974 | 20.038338 | -76.171391 | -302.778225 | 98.313731 | 1274.19842 | 2.429784 |
| 3 | 269.413245 | 2.430419 | 1274.19842 | -4.55 | 8514.893821 | 54.372806 | 66.179254 | 0.418017 | 13.570265 | 0.0 | ... | 805.512369 | 167.854064 | 120.501975 | 26.603042 | 20.039897 | -76.177315 | -302.748629 | 98.299772 | 1274.19842 | 2.430419 |
| 4 | 269.824199 | 2.431081 | 1274.19842 | -4.55 | 8515.440843 | 54.410113 | 66.181719 | 0.418033 | 13.570770 | 0.0 | ... | 805.643038 | 167.867678 | 120.511748 | 26.605200 | 20.041522 | -76.183493 | -302.717751 | 98.285210 | 1274.19842 | 2.431081 |
5 rows × 41 columns
[18]:
# set up figure layout
fig, (ax1, ax2, ax3) = plt.subplots(1,3, figsize = (12, 5))
# loop through main Dictionary
for key, color in zip(Helens_decompress_closed, colors):
# extract individial calculations
res = Helens_decompress_closed[key].copy()
# Lets calculate viscosity now
Viscosity=pt.calculate_viscosity_giordano_2008(liq_comps=res['liquid1'],
T_K=res['All']['T_C']+273.15,
H2O_Liq=res['liquid1']['H2O_Liq'])
# plot fluid chemistry vs pressure
ax1.plot(Viscosity['n_melt'],
res['All']['P_bar'],'-', color=color, label=key)
ax2.plot(res['liquid1']['SiO2_Liq'], Viscosity['n_melt'], '-', color=color)
ax3.plot(res['liquid1']['H2O_Liq'], Viscosity['n_melt'], '-', color=color)
ax1.set_xscale('log')
ax1.set_xlabel('Melt Viscosity (PaS)')
ax1.set_ylabel('Pressure (bar)')
ax2.set_yscale('log')
ax2.set_ylabel('Melt Viscosity (PaS)')
ax2.set_xlabel('Melt SiO$_2$ (wt%)')
ax3.set_yscale('log')
ax3.set_ylabel('Melt Viscosity (PaS)')
ax3.set_xlabel('Melt H$_2$O (wt%)')
fig.tight_layout()
Lets plot the Plag compositions on a ternary diagram to compare to measured microlite chemistry
[19]:
fig, tax = pt.plot_fspar_classification(figsize=(6, 6), major_grid=True, ticks=False, labels=True)
for key in Helens_decompress_closed:
# extract individial calculations
plag = Helens_decompress_closed[key]['plagioclase1']
# calculate plag compositions
tern_points=pt.tern_points_fspar(fspar_comps=Helens_decompress_closed[key]['plagioclase1'])
tax.scatter(
tern_points,
edgecolor="k",
marker="o",
s=50
)
[ ]: