This page was generated from docs/Examples/EOS_calculations/Example5a_Introducing_EOS_Calcs.ipynb. Interactive online version: Binder badge.

Python Notebook Download

Calculations using CO2 EOS

  • DiadFit includes the Span and Wagner (1996) and Sterner and Pitzer (1994) EOS.

  • This allows lots of useful calculations to be performed which rely on these equation of states

  • Spreadsheets we use can be downloaded here:

https://github.com/PennyWieser/DiadFit/blob/main/docs/Examples/EOS_calculations/Batch_processing_TC.xlsx

https://github.com/PennyWieser/DiadFit/blob/main/docs/Examples/EOS_calculations/FI_densities.xlsx

  • If you want to use the Span and Wanger (1996) EOS you will have to download coolprop separatly. Either pip install CoolProp or if you use anaconda, see these!

https://anaconda.org/conda-forge/coolprop

[1]:

#!pip install --upgrade DiadFit

[2]:

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import DiadFit as pf
# Make sure this reads at least 0.0.56
pf.__version__

[2]:

'0.0.91'

Set some plotting parameters

[3]:

plt.rcParams["font.family"] = 'arial'
plt.rcParams["font.size"] =12
plt.rcParams["mathtext.default"] = "regular"
plt.rcParams["mathtext.fontset"] = "dejavusans"
plt.rcParams['patch.linewidth'] = 1
plt.rcParams['axes.linewidth'] = 1
plt.rcParams["xtick.direction"] = "in"
plt.rcParams["ytick.direction"] = "in"
plt.rcParams["ytick.direction"] = "in"
plt.rcParams["xtick.major.size"] = 6 # Sets length of ticks
plt.rcParams["ytick.major.size"] = 4 # Sets length of ticks
plt.rcParams["ytick.labelsize"] = 12 # Sets size of numbers on tick marks
plt.rcParams["xtick.labelsize"] = 12 # Sets size of numbers on tick marks
plt.rcParams["axes.titlesize"] = 14 # Overall title
plt.rcParams["axes.labelsize"] = 14 # Axes labels

There are 3 main classes of calculations you can do using an EOS

  • There are 3 variables related by the EOS, T, density and pressure. If you know 2 of these, you can calculate the third:

  1. Calculate Pressure for a given T and CO2 density, using the function ‘calculate_P_for_rho_T’, can use EOS=’SP94’ or ‘SW96’

  2. Calculate CO2 density for a given P-T. using the function ‘calculate_rho_for_P_T’ with EOS=’SP94’ or ‘SW96’

  3. Calculate Temperature for a specified CO2 density and P using the function ‘calculate_T_for_rho_P’ with EOS=’SP94’ or ‘SW96’

You also have to use equation of states for microthermometry

  1. calculate CO2 density given a homogenization temp. This requires knowledge of the position of the L-V phase boundary. This is only currently supported in Span and Wanger (1996), using the function ‘calculate_CO2_density_homog_T’

Calc 1 - Calculating Pressures from densities and entrapment temps

  • The fundamental tenet of fluid inclusion barometry is that your inclusion has a fixed mass and volume. This means if you can get its density today in the lab using microthermometry or Raman Spectroscopy, and you can choose an entrapment temperature, you can calculate a pressure

  • We can extrapolate along an isopleth using Span and Wanger (1996), or Sterner and Pitzer (1994)

  • This uses the function ‘calculate_P_for_rho_T’

  • Lets perform the calculation for a single fluid inclusion using the Span and Wanger (1996) EOS. This requires you to have installed CoolProp!

[3]:

P_SW96=pf.calculate_P_for_rho_T(
CO2_dens_gcm3=0.5, T_K=1200, EOS='SW96')
P_SW96

[3]:
P_kbar P_MPa T_K CO2_dens_gcm3
0 1.721807 172.180688 1200 0.5

  • Lets use Sterner and Pitzer (1994) instead - if you cant install CoolProp, this is the EOS for you!

[4]:

P_SP94=pf.calculate_P_for_rho_T(
CO2_dens_gcm3=0.5, T_K=1200, EOS='SP94')
P_SP94

[4]:
P_kbar P_MPa T_K CO2_dens_gcm3
0 1.756787 175.678735 1200 0.5

  • Now lets calculate Depth using SP94, using the crustal density model of Ryan-Lerner

[5]:


D=pf.convert_pressure_to_depth(P_kbar=P_SP94['P_kbar'], model='ryan_lerner')
D

[5]:

0    7.139669
dtype: float64

  • Lets perform the calculation for an entire spreadsheet of fluid inclusions

[6]:

# Get from here: https://github.com/PennyWieser/DiadFit/blob/main/docs/Examples/EOS_calculations/FI_densities.xlsx
df=pd.read_excel('FI_densities.xlsx')
P_kbar_SW96=pf.calculate_P_for_rho_T(
T_K=df['Temp in C']+273.15,
CO2_dens_gcm3=df['Density_g_cm3'],
EOS='SW96', Sample_ID=df['Sample'])
P_kbar_SW96.head()

[6]:
P_kbar P_MPa T_K CO2_dens_gcm3 Sample_ID
0 0.276826 27.682633 1373.15 0.1 FI1
1 1.559737 155.973743 1473.15 0.4 FI2
2 3.789037 378.903651 1423.15 0.7 FI3
3 10.626112 1062.611174 1473.15 1.1 FI4
4 26.062249 2606.224906 1573.15 1.5 FI5

  • Lets convert these pressures into depths using the Ryan-Lerner density profile

  • Note, Nans are returned for pressures outside the calibration range of this model

[7]:

RL_D=pf.convert_pressure_to_depth(P_kbar=P_kbar_SW96['P_kbar'],
model='ryan_lerner')
RL_D

[7]:

0     1.256957
1     6.419763
2    14.157451
3          NaN
4          NaN
dtype: float64

  • Lets make synthetic data across a range of temperatures and densities (between 0-2000 K, 0.1 and 1.1 g/cm3) to visualize how pressure and temperature are related

[8]:

# Chose an array of temperatures in celcius
T_C=np.linspace(34, 2000, 100)
# Convert to Kelvin
T_K=T_C+273.15
# Choose an array of densities (g/cm3)
rho_lin=np.linspace(1.1, 0.1, 11)

# Now write a loop to fill out a pressure array
Pressure_SP94=np.empty([len(rho_lin), len(T_K)], float)
Pressure_SW96=np.empty([len(rho_lin), len(T_K)], float)
for i in range(0, len(rho_lin)):
    P_SW96=pf.calculate_P_for_rho_T(CO2_dens_gcm3=rho_lin[i], T_K=T_K, EOS='SW96').P_kbar
    Pressure_SW96[i, :]=P_SW96
    P_SP94=pf.calculate_P_for_rho_T(CO2_dens_gcm3=rho_lin[i], T_K=T_K, EOS='SP94').P_kbar
    Pressure_SP94[i, :]=P_SP94

  • Now lets plot it up!

[9]:

fig, (ax1) = plt.subplots(1, 1, figsize=(4,3), sharex=True, sharey=True)
import matplotlib.colors as colors
#ax1.plot(T_C, P_kbar, '-r')
import matplotlib.cm as cm
for i in range(0, len(rho_lin)):
    color=cm.plasma(rho_lin[i])
    ax1.plot(T_C, Pressure_SW96[i, :], '-',
    color=color,label=str(np.round(rho_lin[i], 1)) + ' g/cm$^3$')


norm = colors.Normalize(vmin=rho_lin[0], vmax=rho_lin[-1])
sm = cm.ScalarMappable(cmap=cm.plasma, norm=norm)
sm.set_array([])
cbar = fig.colorbar(sm, ax=ax1)
cbar.set_label('Density (g/cm$^3$)')



ax1.set_ylabel('Pressure (kbar)')
#legend = ax1.legend(loc='upper left', ncol=2, frameon=True, fontsize=7)
ax1.set_xlim([0, 2000])
ax1.set_ylim([0, 16])

ax1.set_xlabel('T (°C)')
ax1.grid(color = 'k', linestyle = '-', linewidth = 1, alpha = 0.1)

fig.savefig('EOS.png', dpi=200, bbox_inches='tight')
../../_images/Examples_EOS_calculations_Example5a_Introducing_EOS_Calcs_20_0.png

  • That is for one EOS, but we might want to compare SW96 and SP94 - we do that here

[10]:

fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(14,5), sharex=True, sharey=True)
#ax1.plot(T_C, P_kbar, '-r')
import matplotlib.cm as cm
for i in range(0, len(rho_lin)):
    color=cm.plasma(i/len(rho_lin))
    ax1.plot(T_C, Pressure_SP94[i, :], '-.',
    color=color,label=str(np.round(rho_lin[i], 1)) + ' g/cm$^3$')
    ax2.plot(T_C, Pressure_SW96[i, :], '-',
    color=color,label=str(np.round(rho_lin[i], 1)) + ' g/cm$^3$')

    ax3.plot(T_C, Pressure_SW96[i, :], '-',
    color=color,label=str(np.round(rho_lin[i], 1)) + ' g/cm$^3$')
    ax3.plot(T_C, Pressure_SP94[i, :], '-.',
    color=color, label=str(np.round(rho_lin[i], 1)) + ' g/cm$^3$', alpha=0.5)


ax1.set_ylabel('Pressure (kbar)')
legend = ax1.legend(loc='upper left', ncol=2, frameon=True, fontsize=7)
ax1.set_xlim([0, 2000])
ax1.set_ylim([0, 16])
ax2.set_xlabel('T (°C)')
ax3.set_xlabel('T (°C)')
ax1.set_xlabel('T (°C)')
ax1.grid(color = 'k', linestyle = '-', linewidth = 1, alpha = 0.1)
ax2.grid(color = 'k', linestyle = '-', linewidth = 1, alpha = 0.1)
ax3.grid(color = 'k', linestyle = '-', linewidth = 1, alpha = 0.1)

ax2.yaxis.set_tick_params(which='both', labelbottom=True)
ax3.yaxis.set_tick_params(which='both', labelbottom=True)
ax1.set_title('Sterner and Pitzer (1994)')
ax2.set_title('Span and Wanger (1996)')
ax3.set_title('Comparing SP1994 and SW1996')

[10]:

Text(0.5, 1.0, 'Comparing SP1994 and SW1996')
../../_images/Examples_EOS_calculations_Example5a_Introducing_EOS_Calcs_22_1.png

Calc 2 - Calculating density for given Pressure and temperature

  • This might be useful for experimental petrologists, or anyone using a gas cell to calibrate their Raman (where they have measured P and T inside the cell).

  • Lets do the calculation firstly for one pressure and temperature we are interested in

  • You could also load in a spreadsheet as before

[11]:

SW96_rho=pf.calculate_rho_for_P_T(P_kbar=1,
T_K=1400, EOS='SW96')
SW96_rho

[11]:

0    0.300608
dtype: float64

  • Or for lots and lots of different pressures, to see how density changes with pressure at fixed temperatures

[12]:

P_kbar_linspace=np.linspace(0.01, 10, 50)
SW96_Parray_1400=pf.calculate_rho_for_P_T(P_kbar=P_kbar_linspace,
T_K=1400, EOS='SW96')
SP94_Parray_1400=pf.calculate_rho_for_P_T(P_kbar=P_kbar_linspace,
T_K=1400, EOS='SP94')

SW96_Parray_1200=pf.calculate_rho_for_P_T(P_kbar=P_kbar_linspace,
T_K=1200, EOS='SW96')
SP94_Parray_1200=pf.calculate_rho_for_P_T(P_kbar=P_kbar_linspace,
T_K=1200, EOS='SP94')

[13]:

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12,5), sharex=True)

ax1.plot(P_kbar_linspace, SW96_Parray_1400, '-r', label='SW96 1400K')
ax1.plot(P_kbar_linspace, SP94_Parray_1400, ':r', label='SP94 1400K')


ax1.plot(P_kbar_linspace, SW96_Parray_1200, '-k', label='SW96 1200K')
ax1.plot(P_kbar_linspace, SP94_Parray_1200, ':', color='k', label='SP94 1200K')
ax1.legend()

ax2.plot(P_kbar_linspace, SW96_Parray_1400-SP94_Parray_1400, '--r', label='1400 K')
ax2.plot(P_kbar_linspace, SW96_Parray_1200-SP94_Parray_1200, '--k', label='1200 K')
ax2.legend()


ax1.set_xlabel('Pressure (kbar)')
ax1.set_ylabel('Density (g/cm3)')

ax2.set_xlabel('Pressure (kbar)')
ax2.set_ylabel('Density SW96-SP94 (g/cm3)')

[13]:

Text(0, 0.5, 'Density SW96-SP94 (g/cm3)')
../../_images/Examples_EOS_calculations_Example5a_Introducing_EOS_Calcs_28_1.png

  • We can also use this to compare EOS in another way

[14]:

fig, (( ax3, ax4)) = plt.subplots(1, 2, figsize=(12,5))
ax5 = fig.add_axes([0.05, 0.95, 0.9, 0.05])



#ax1.plot(T_C, P_kbar, '-r')
import matplotlib.cm as cm
for i in range(0, len(rho_lin)):
    color=cm.plasma(i/len(rho_lin))
    if i==0:
        ax3.plot(T_C, Pressure_SW96[i, :], '-',
    color=color,label='SW96')
        ax3.plot(T_C, Pressure_SP94[i, :], '-.',
            color=color, alpha=0.5, label='SP94')


    ax3.plot(T_C, Pressure_SW96[i, :], '-',
    color=color)
    ax3.plot(T_C, Pressure_SP94[i, :], '-.',
    color=color, alpha=0.5)

    ax4.plot(T_C,  100*(Pressure_SP94[i, :]-Pressure_SW96[i, :])/(0.5*(Pressure_SP94[i, :]+Pressure_SW96[i, :])),  label=str(np.round(rho_lin[i], 1)) + ' g/cm$^3$',color=color)
    ax5.plot(0*T_C,  0*100*(Pressure_SP94[i, :]-Pressure_SW96[i, :])/(0.5*(Pressure_SP94[i, :]+Pressure_SW96[i, :])),  label=str(np.round(rho_lin[i], 1)) + ' g/cm$^3$', color=color)

legend = ax5.legend(loc='center', ncol=6, frameon=False, fontsize=12)
ax5.axis('off')
ax4.plot([33, 2000], [0, 0], '-k')



ax3.legend(loc='upper center')

ax3.set_ylabel('Pressure (kbar)')
ax4.set_ylabel('% offset (SP94-SW96)')
#legend = ax1.legend(loc='upper left', ncol=2, frameon=True, fontsize=9)
ax4.set_xlim([33, 2000])
ax3.set_xlim([33, 2000])
ax2.set_ylim([0, 16])

ax3.set_xlabel('T (°C)')

ax4.set_xlabel('T (°C)')


ax3.grid(color = 'k', linestyle = '-', linewidth = 1, alpha = 0.1)

ax2.yaxis.set_tick_params(which='both', labelbottom=True)
ax3.yaxis.set_tick_params(which='both', labelbottom=True)



ax1.set_title('Sterner and Pitzer (1994)')
ax2.set_title('Span and Wanger (1996)')

plt.subplots_adjust(top=0.89, bottom=0.1, left=0.1, right=0.9, hspace=0.4, wspace=0.2)
ax3.annotate("a)", xy=(0.02, 0.95), xycoords="axes fraction", fontsize=14)
ax4.annotate("b)", xy=(0.02, 0.95), xycoords="axes fraction", fontsize=14)


fig.savefig('EOS_diffs.png', dpi=200, bbox_inches='tight')
../../_images/Examples_EOS_calculations_Example5a_Introducing_EOS_Calcs_30_0.png

Calc 3 - Calculating temperature for a known density and pressure

  • This function behaves the same as the previous 2 discussed, here we show a single calculation, it could also be applied to a spreadsheet

[15]:

SW96_T=pf.calculate_T_for_rho_P(CO2_dens_gcm3=0.3, P_kbar=1,
 EOS='SW96')
SW96_T

[15]:

0    1403.252835
dtype: float64

[16]:

SP94_T=pf.calculate_T_for_rho_P(CO2_dens_gcm3=0.3, P_kbar=1,
 EOS='SP94')
SP94_T

[16]:

0    1389.233164
dtype: float64

Calc 4: Calculating densities in 2 phase bubbles at room temp

  • We have parameterized the L-V curve using the Span and Wanger (1996) EOS.

  • This means you can enter a temperature, and find out what the density of the CO2 gas and liquid would be for an inclusion following the L-V curve on the phase diagram

  • This can be very useful, because if you are conducting Raman analyses at room temperature for bubbles with densities higher than the critical density, you will only analyse the interior gas (with complex interactions with the liquid phase, see DeVitre et al. 2023, Volcanica).

[18]:

# At 18C, this tells us the coexisting liquid is 0.79 g/cm3, and gas is 0.179 g/cm3
df_18_SW96=pf.calculate_CO2_density_homog_T(T_h_C=18, EOS='SW96')
df_18_SW96

[18]:
Bulk_gcm3 Liq_gcm3 Gas_gcm3 T_h_C homog_to
0 NaN 0.793752 0.179577 18 None 
[20]:

# Lets do the same calc but at 22C
df_22_SW96=pf.calculate_CO2_density_homog_T(T_h_C=22, EOS='SW96')
df_22_SW96

[20]:
Bulk_gcm3 Liq_gcm3 Gas_gcm3 T_h_C homog_to
0 NaN 0.750772 0.211079 22 None 
[21]:

# And at 20C
df_20_SW96=pf.calculate_CO2_density_homog_T(T_h_C=20, EOS='SW96')
df_20_SW96

[21]:
Bulk_gcm3 Liq_gcm3 Gas_gcm3 T_h_C homog_to
0 NaN 0.773392 0.194204 20 None 
[22]:

 # Lets consider 'room Temps' between 18 and 23C using Span and Wagner 1996
Temps=np.linspace(17, 28)
df_roomT_SW96=pf.calculate_CO2_density_homog_T(T_h_C=Temps, EOS='SW96')
df_roomT_SW96.head()

[22]:
Bulk_gcm3 Liq_gcm3 Gas_gcm3 T_h_C homog_to
0 NaN 0.803262 0.172932 17.000000 None
1 NaN 0.801161 0.174389 17.224490 None
2 NaN 0.799041 0.175866 17.448980 None
3 NaN 0.796902 0.177363 17.673469 None
4 NaN 0.794741 0.178881 17.897959 None 
[23]:

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10,5))
ax1.plot(Temps, df_roomT_SW96['Liq_gcm3'], '-r', label='SW96')
ax2.plot(Temps, df_roomT_SW96['Gas_gcm3'], '-b', label='SW96')

ax1.set_xlabel('Temp (C)')
ax1.set_ylabel('Liq density (g/cm3)')
ax2.set_xlabel('Temp (C)')
ax2.set_ylabel('Gas density (g/cm3)')

[23]:

Text(0, 0.5, 'Gas density (g/cm3)')
../../_images/Examples_EOS_calculations_Example5a_Introducing_EOS_Calcs_39_1.png

Calc 5 - Calculating homogenization temperatures from microthermometry

  • Microthermometry works by cooling CO2- rich fluid inclusions down, and then heating up until the liquid and gas homogenizes to a single phase. This gives you the density of CO2. You can then calculate the pressure by extrapolating

  • Again, we use the Span and Wagner (1996) EOS. We note that some existing spreadsheets use the L-V curve from Span and Wanger, but extrapolate using Sterner and Pitzer (1994). We do not wish to mix and match EOS, so use SW96 for the L-V curve and extrapolating along the isochore

  • Lets perform calcs for one FI, e.g. while we are on the microscope and are very excited to get the density!

  • Say we see the homogenization at -18C, and it homogenizes to a liquid, we know the density is 1.022 g/cm3

[24]:

Density1=pf.calculate_CO2_density_homog_T(T_h_C=-18, Sample_ID='FI2',
                            homog_to='L', EOS='SW96')
Density1

[24]:
Bulk_gcm3 Liq_gcm3 Gas_gcm3 T_h_C homog_to Sample_ID
0 1.022329 1.022329 0.055156 -18 L FI2

Lets load an entire day of microthermometry results, and calculate density for all at once

[25]:

# Get from here: https://github.com/PennyWieser/DiadFit/blob/main/docs/Examples/EOS_calculations/Batch_processing_TC.xlsx
homog_t=pd.read_excel('Batch_processing_TC.xlsx')
homog_t.head()

Whole_day_dens=pf.calculate_CO2_density_homog_T(T_h_C=homog_t['T_c_Homog'],
            Sample_ID=homog_t['Sample'], homog_to=homog_t['homog_to'], EOS='SW96')
Whole_day_dens.head()

[25]:
Bulk_gcm3 Liq_gcm3 Gas_gcm3 T_h_C homog_to Sample_ID
0 0.593263 0.593263 0.345003 30 L FI1
1 0.710471 0.710471 0.242710 25 L FI2
2 0.773392 0.773392 0.194204 20 L FI3
3 0.160735 0.821187 0.160735 15 V FI4
4 0.861096 0.861096 0.135157 10 L FI5

Propagating errors in microthermometry

  • Lets say we measure a homogenization temp of -30 C, but there is a 0.3 C error on our stage temp, or we werent exactly sure when the bubble disapeared

  • N dup says how many monte carlo simulations we are going to do. Lets say 1000.

[26]:

# First calculate density
Av_outputs_SW96, All_outputs_SW96=pf.propagate_microthermometry_uncertainty(T_h_C=-30,
Sample_ID='test', error_T_h_C=0.3, N_dup=1000,
EOS='SW96', homog_to='L',
error_dist_T_h_C='normal', error_type_T_h_C='Abs')

Lets inspect the outputs

[28]:

# This is the average for all 1000 simulations
Av_outputs_SW96.head()

[28]:
Sample_ID Density_Gas_noMC Density_Liq_noMC Mean_density_Gas_gcm3 Std_density_Gas_gcm3 Std_density_Gas_gcm3_from_percentiles Std_density_Liq_gcm3_from_percentiles Mean_density_Liq_gcm3 Std_density_Liq_gcm3 Input_temp error_T_h_C
0 test 0.0371 1.075785 0.037109 0.000371 0.000364 0.001219 1.075759 0.001241 -30 0.3 
[29]:

# This gives you the result of each simulation!
All_outputs_SW96.head()

[29]:
Bulk_gcm3 Liq_gcm3 Gas_gcm3 T_h_C homog_to Sample_ID
0 1.075155 1.075155 0.037289 -29.850986 L test
1 1.075960 1.075960 0.037048 -30.041479 L test
2 1.074963 1.074963 0.037346 -29.805693 L test
3 1.073850 1.073850 0.037681 -29.543091 L test
4 1.076082 1.076082 0.037012 -30.070246 L test

Lets plot up these outputs

[30]:

fig, (ax0, ax1, ax2) = plt.subplots(1, 3, figsize=(14,5))
ax0.hist(All_outputs_SW96['T_h_C'], ec='k');
ax1.hist(All_outputs_SW96['Liq_gcm3'], ec='k');
ax2.hist(All_outputs_SW96['Gas_gcm3'], ec='k');
ax0.set_xlabel('Input T_C')
ax1.set_xlabel('Liq (g/cm3)')
ax2.set_xlabel('gas (g/cm3)')

[30]:

Text(0.5, 0, 'gas (g/cm3)')
../../_images/Examples_EOS_calculations_Example5a_Introducing_EOS_Calcs_51_1.png