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

Python Notebook Download

Olivine-Liquid thermometry

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

[1]:
#!pip install Thermobar
[2]:
# Loading various python things
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import Thermobar as pt
[3]:
pt.__version__
[3]:
'1.0.3dev'
[4]:
# Setting plotting parameters
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

Step 1 - load data

[5]:
out=pt.import_excel('Liquid_only_Thermometry.xlsx', sheet_name="Ol-Liq")
my_input=out['my_input']
myLiquids1=out['Liqs']
myOls1=out['Ols']
display(myOls1.head())
display(myLiquids1.head())
SiO2_Ol TiO2_Ol Al2O3_Ol FeOt_Ol MnO_Ol MgO_Ol CaO_Ol Na2O_Ol K2O_Ol Cr2O3_Ol NiO_Ol Sample_ID_Ol
0 40.5 0.02 0.08 12.40 0.17 47.4 0.30 0.0 0 0.03 0.0 0
1 41.3 0.03 0.11 9.59 0.14 50.2 0.31 0.0 0 0.00 0.0 1
2 39.7 0.05 0.11 15.60 0.18 44.5 0.31 0.0 0 0.03 0.0 2
3 40.5 0.05 0.10 13.20 0.18 46.8 0.29 0.0 0 0.02 0.0 3
4 40.5 0.00 0.10 9.41 0.10 49.3 0.31 0.0 0 0.00 0.0 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 NiO_Liq CoO_Liq CO2_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.0 0.0 0.0 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 0.0 0.0 0.0 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 0.0 0.0 0.0 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 0.0 0.0 0.0 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 0.0 0.0 0.0 4

Example 1 - Olivine-Liquid temperatures

  • It has been shown many times that olivines are not in equilirium with their co-erupted carrier melts. Thus, it is difficult to know what a meaningful Ol-Liq match is. Because of this, by default we return the value of Kd calculated for olivine-liquid as well as the calculated temperature in the “calculate_ol_liq_temp” function.

  • The function uses the Fe3Fet_Liq column in the user-entered spreadsheet, and calculates Fe/Mg in the liquid using only Fe2+ by default.

  • You can also specify a Fe3Fet_Liq ratio in the function itself (see later examples) to overwrite this value.

  • If the user doesn’t specify a column for Fe3Fet in the input, this value is set to zero.

1a - Putirka eq 22

[6]:
eq22_PHinput=pt.calculate_ol_liq_temp(liq_comps=myLiquids1, ol_comps=myOls1,  equationT="T_Put2008_eq22",
                                      P=my_input['P_kbar'])
eq22_PHinput
[6]:
T_K_calc Kd (Fe-Mg) Meas
0 1289.947705 0.314264
1 1229.813416 0.175383
2 1285.857491 0.269582
3 1198.240159 0.173799
4 1259.174284 0.152172
5 1244.931462 0.159554
6 1241.300118 0.215915
7 1199.601831 0.404699
8 1162.693468 0.314115
9 1339.466029 0.429135

1b - Specify eq_tests=True

  • FIf you specify eq_tests=True, the function also returns calculated Kd values for olivine-Liquid pairs so you can assess if olivine-liquids are in equilibrium.

  • A number of different equilibrium tests are returned, using Toplis, Matzen and Roeder and Emslie (preferred value=0.3 for Roeder and Emslie, 1970, =0.34 for Matzen, function of melt comp, temp and press for Toplis).

  • As before, here the Fe3Fet_Liq ratio in the spreadsheet is being used for calculations.

[7]:
eq22_PHinput=pt.calculate_ol_liq_temp(liq_comps=myLiquids1, ol_comps=myOls1,
                                      equationT="T_Put2008_eq22", P=my_input['P_kbar'], eq_tests=True)
eq22_PHinput
[7]:
T_K_calc Kd Meas Kd calc (Toplis) ΔKd, Toplis (M-P) ΔKd, Roeder (M-P) ΔKd, Matzen (M-P) SiO2_Liq TiO2_Liq Al2O3_Liq FeOt_Liq ... K2O_Ol Cr2O3_Ol NiO_Ol Sample_ID_Ol P2O5_Ol DMg_Meas CNML CSiO2L NF Den_Beat93
0 1289.947705 0.314264 0.325040 -0.010776 0.014264 -0.025736 57.023602 0.623106 16.332899 4.36174 ... 0 0 0.0 0 0 9.472619 0.170395 0.557471 -0.761843 9.819071
1 1229.813416 0.175383 0.308684 -0.133301 -0.124617 -0.164617 57.658600 0.654150 17.194799 3.90621 ... 0 0 0.0 1 0 14.186789 0.137941 0.570023 -0.816654 10.303635
2 1285.857491 0.269582 0.315891 -0.046309 -0.030418 -0.070418 60.731201 0.862054 17.144199 4.07781 ... 0 0 0.0 2 0 15.531216 0.122184 0.579288 -0.792785 10.250501
3 1198.240159 0.173799 0.296719 -0.122920 -0.126201 -0.166201 61.532799 0.440860 16.508801 3.32990 ... 0 0 0.0 3 0 23.753206 0.097626 0.605753 -0.767076 10.715090
4 1259.174284 0.152172 0.323766 -0.171594 -0.147828 -0.187828 52.969101 0.803412 17.563000 5.93217 ... 0 0 0.0 4 0 10.781292 0.186877 0.523030 -0.842157 10.215321
5 1244.931462 0.159554 0.324986 -0.165432 -0.140446 -0.180446 54.050201 0.857348 17.333300 5.60072 ... 0 0 0.0 5 0 12.276130 0.168070 0.533835 -0.833465 10.295063
6 1241.300118 0.215915 0.333878 -0.117963 -0.084085 -0.124085 55.656300 0.897984 17.117800 5.31307 ... 0 0 0.0 6 0 14.083460 0.145845 0.547483 -0.820914 10.324021
7 1199.601831 0.404699 0.391378 0.013321 0.104699 0.064699 49.054699 0.488832 14.665000 5.83168 ... 0 0 0.0 7 0 7.372298 0.220053 0.554209 -0.789502 9.845140
8 1162.693468 0.314115 0.361834 -0.047718 0.014115 -0.025885 50.625099 0.334074 16.875000 5.01968 ... 0 0 0.0 8 0 10.509878 0.180293 0.551305 -0.873502 10.229233
9 1339.466029 0.429135 0.299720 0.129416 0.129135 0.089135 51.403301 0.663880 18.019600 5.98440 ... 0 0 0.0 9 0 7.070601 0.228386 0.498172 -0.840523 9.673734

10 rows × 74 columns

  • If you want to access just the temperature from this panda dataframe, you do name[‘column heading’]

  • To get the output in Celcius, subtract -273.15

[7]:
Teq22_PHinput=eq22_PHinput['T_K_calc']-273.15 #converting temp to C
Teq22_PHinput
[7]:
0    1016.797705
1     956.663416
2    1012.707491
3     925.090159
4     986.024284
5     971.781462
6     968.150118
7     926.451831
8     889.543468
9    1066.316029
Name: T_K_calc, dtype: float64
  • Can also filter outputs, to only look at temps for pairs with delta K<=0.03 using the pandas loc function.

  • This basically says “give me the rows of eq22_PHinput where the column “ΔKd, Roeder” is less than or equal to 0.03 (e.g., in equilibrium)

  • Only 2 Ol-Liq pairs are in equilibrium here!

[8]:
T_in_eq=eq22_PHinput.loc[eq22_PHinput['ΔKd, Roeder']<=0.03]
T_in_eq
[8]:
T_K_calc Kd Meas Kd calc (Toplis) ΔKd, Toplis ΔKd, Roeder ΔKd, Matzen SiO2_Liq TiO2_Liq Al2O3_Liq FeOt_Liq ... K2O_Ol Cr2O3_Ol NiO_Ol Sample_ID_Ol P2O5_Ol DMg_Meas CNML CSiO2L NF Den_Beat93
0 1289.947705 0.314264 0.325040 0.010776 0.014264 0.025736 57.023602 0.623106 16.332899 4.36174 ... 0 0 0.0 0 0 9.472619 0.170395 0.557471 -0.761843 9.819071
8 1162.693468 0.314115 0.361834 0.047718 0.014115 0.025885 50.625099 0.334074 16.875000 5.01968 ... 0 0 0.0 8 0 10.509878 0.180293 0.551305 -0.873502 10.229233

2 rows × 74 columns

1c - You can also overwrite the spreadsheet Fe3Fet_Liq in the function itself

  • Here, we peform calculations using 30% Fe3+

[9]:
eq22_PHinput_30Fe3=pt.calculate_ol_liq_temp(liq_comps=myLiquids1, ol_comps=myOls1,
                                      equationT="T_Put2008_eq22", P=my_input['P_kbar'],
                                        Fe3Fet_Liq=0.3)
eq22_PHinput_30Fe3
[9]:
T_K_calc Kd (Fe-Mg) Meas
0 1289.947705 0.359158
1 1229.813416 0.200438
2 1285.857491 0.308094
3 1198.240159 0.198628
4 1259.174284 0.173911
5 1244.931462 0.182347
6 1241.300118 0.246761
7 1199.601831 0.462513
8 1162.693468 0.358989
9 1339.466029 0.490441

Example 2 - Considering all liquid-olivine matches

  • In reality, you have probably measured a load of glasses, and a load of olivines, but you dont know which ones match with which.

  • So, we have a matching algorithm, that returns all possible olivine-liquid pairs, which you can then filter.

[13]:
eq22_PHinput_30Fe3_Matching=pt.calculate_ol_liq_temp_matching(liq_comps=myLiquids1, ol_comps=myOls1,
                                      equationT="T_Put2008_eq22", P=5,
                                        Fe3Fet_Liq=0.3, eq_tests=True)
eq22_PHinput_30Fe3_Matching
Considering N=10 Ol & N=10 Liqs, which is a total of N=100 Liq-Ol pairs, be patient if this is >>1 million!
[13]:
T_K_calc Kd Meas Kd calc (Toplis) ΔKd, Toplis ΔKd, Roeder ΔKd, Matzen SiO2_Liq TiO2_Liq Al2O3_Liq FeOt_Liq ... Al_Ol_cat_frac Na_Ol_cat_frac K_Ol_cat_frac Mn_Ol_cat_frac Ti_Ol_cat_frac DMg_Meas CNML CSiO2L NF Den_Beat93
0 1306.094822 0.314264 0.331693 0.017429 0.014264 0.025736 57.023602 0.623106 16.332899 4.36174 ... 0.000772 0.0 0.0 0.001179 0.000123 9.490384 0.170076 0.556427 -0.760264 9.813745
1 1246.358585 0.240169 0.320346 0.080177 0.059831 0.099831 57.658600 0.654150 17.194799 3.90621 ... 0.000772 0.0 0.0 0.001179 0.000123 13.711088 0.137693 0.568999 -0.815019 10.226597
2 1278.487099 0.201173 0.308304 0.107131 0.098827 0.138827 60.731201 0.862054 17.144199 4.07781 ... 0.000772 0.0 0.0 0.001179 0.000123 16.264149 0.121869 0.577797 -0.790522 10.330149
3 1196.823888 0.161199 0.295190 0.133991 0.138801 0.178801 61.532799 0.440860 16.508801 3.32990 ... 0.000772 0.0 0.0 0.001179 0.000123 24.065168 0.097450 0.604665 -0.765547 10.732471
4 1265.970449 0.208562 0.331112 0.122550 0.091438 0.131438 52.969101 0.803412 17.563000 5.93217 ... 0.000772 0.0 0.0 0.001179 0.000123 10.410636 0.186517 0.522022 -0.840343 10.135819
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
95 1258.696461 0.280399 0.341138 0.060738 0.019601 0.059601 54.050201 0.857348 17.333300 5.60072 ... 0.000589 0.0 0.0 0.001411 0.000188 11.366691 0.167758 0.532846 -0.831744 10.131992
96 1243.573070 0.241747 0.337524 0.095777 0.058253 0.098253 55.656300 0.897984 17.117800 5.31307 ... 0.000589 0.0 0.0 0.001411 0.000188 13.960289 0.145510 0.546226 -0.818818 10.295163
97 1195.069349 0.351613 0.385106 0.033493 0.051613 0.011613 49.054699 0.488832 14.665000 5.83168 ... 0.000589 0.0 0.0 0.001411 0.000188 7.596985 0.220036 0.554167 -0.789435 9.904811
98 1160.902635 0.302323 0.360085 0.057761 0.002323 0.037677 50.625099 0.334074 16.875000 5.01968 ... 0.000589 0.0 0.0 0.001411 0.000188 10.651575 0.180240 0.551142 -0.873211 10.254547
99 1339.482412 0.429135 0.299722 0.129414 0.129135 0.089135 51.403301 0.663880 18.019600 5.98440 ... 0.000589 0.0 0.0 0.001411 0.000188 7.071114 0.228369 0.498136 -0.840454 9.673521

100 rows × 76 columns

Example 3 - Calculating equilibrium olivine Fo contents based on the liquid composition

3a - Roeder and Emslie, 1970

  • Here, using the Kd model of Roeder and Emslie, 1970, Kd=0.03+-0.03, for 20% Fe3

  • The function returns the Mg# of the liquid using Fe2+ based on the Fe3Fet ratio from the input spreadsheet, or specified by the user in the function (as in this example)

  • The Fet assumes that all Fe in the olivine and liquid is Fe2+

  • The Eq Fo contents are calculated using the Mg# with just Fe2, based on the Kd model shown in brackets

[9]:
Eq_Ol_Roeder=pt.calculate_eq_ol_content(liq_comps=myLiquids1, Kd_model="Roeder1970")
Eq_Ol_Roeder
[9]:
Mg#_Liq_Fe2 Mg#_Liq_Fet Eq Fo (Roeder, Kd=0.3) Eq Fo (Roeder, Kd=0.33) Eq Fo (Roeder, Kd=0.27)
0 0.681666 0.631416 0.877117 0.866470 0.888030
1 0.620707 0.566947 0.845080 0.832188 0.858378
2 0.578196 0.523041 0.820442 0.805970 0.835443
3 0.523445 0.467721 0.785468 0.768971 0.802688
4 0.586967 0.532031 0.825694 0.811549 0.840342
5 0.565465 0.510056 0.812653 0.797709 0.828169
6 0.528731 0.473003 0.789019 0.772716 0.806025
7 0.620033 0.566244 0.844705 0.831788 0.858030
8 0.583864 0.528846 0.823846 0.809585 0.838619
9 0.665729 0.614386 0.869086 0.857856 0.880615
[8]:
Eq_Ol_Roeder=pt.calculate_eq_ol_content(liq_comps=myLiquids1, Kd_model="Roeder1970", Fe3Fet_Liq=0.2)
Eq_Ol_Roeder
[8]:
Mg#_Liq_Fe2 Mg#_Liq_Fet Eq Fo (Roeder, Kd=0.3) Eq Fo (Roeder, Kd=0.33) Eq Fo (Roeder, Kd=0.27)
0 0.681666 0.631416 0.877117 0.866470 0.888030
1 0.620707 0.566947 0.845080 0.832188 0.858378
2 0.578196 0.523041 0.820442 0.805970 0.835443
3 0.523445 0.467721 0.785468 0.768971 0.802688
4 0.586967 0.532031 0.825694 0.811549 0.840342
5 0.565465 0.510056 0.812653 0.797709 0.828169
6 0.528731 0.473003 0.789019 0.772716 0.806025
7 0.620033 0.566244 0.844705 0.831788 0.858030
8 0.583864 0.528846 0.823846 0.809585 0.838619
9 0.665729 0.614386 0.869086 0.857856 0.880615
  • If you also specify the olivine compositions, it will add the measured Fo content as a column for comparison

[11]:
Eq_Ol_Roeder=pt.calculate_eq_ol_content(liq_comps=myLiquids1, ol_comps=myOls1,
                                     Kd_model="Roeder1970", Fe3Fet_Liq=0.2)
Eq_Ol_Roeder
[11]:
Mg#_Liq_Fe2 Mg#_Liq_Fet Eq Fo (Roeder, Kd=0.3) Eq Fo (Roeder, Kd=0.33) Eq Fo (Roeder, Kd=0.27) Fo_meas
0 0.681666 0.631416 0.877117 0.866470 0.888030 0.872023
1 0.620707 0.566947 0.845080 0.832188 0.858378 0.903203
2 0.578196 0.523041 0.820442 0.805970 0.835443 0.835656
3 0.523445 0.467721 0.785468 0.768971 0.802688 0.863386
4 0.586967 0.532031 0.825694 0.811549 0.840342 0.903278
5 0.565465 0.510056 0.812653 0.797709 0.828169 0.890781
6 0.528731 0.473003 0.789019 0.772716 0.806025 0.838609
7 0.620033 0.566244 0.844705 0.831788 0.858030 0.801278
8 0.583864 0.528846 0.823846 0.809585 0.838619 0.817074
9 0.665729 0.614386 0.869086 0.857856 0.880615 0.822724
  • We can also specify to use the Kd model of Matzen (2011), where Kd=0.34 +-0.012

[12]:
Eq_Ol_Matzen=pt.calculate_eq_ol_content(liq_comps=myLiquids1,
                                        Kd_model="Matzen2011", Fe3Fet_Liq=0.2)
Eq_Ol_Matzen
[12]:
Mg#_Liq_Fe2 Mg#_Liq_Fet Eq Fo (Matzen, Kd=0.34) Eq Fo (Matzen, Kd=0.352) Eq Fo (Matzen, Kd=0.328)
0 0.681666 0.631416 0.862978 0.858825 0.867172
1 0.620707 0.566947 0.827977 0.822980 0.833035
2 0.578196 0.523041 0.801259 0.795678 0.806919
3 0.523445 0.467721 0.763625 0.757307 0.770049
4 0.586967 0.532031 0.806941 0.801479 0.812477
5 0.565465 0.510056 0.792848 0.787094 0.798688
6 0.528731 0.473003 0.767431 0.761182 0.773782
7 0.620033 0.566244 0.827569 0.822563 0.832637
8 0.583864 0.528846 0.804941 0.799437 0.810521
9 0.665729 0.614386 0.854176 0.849803 0.858595

3b - Toplis, 2005

  • This Kd expression is a bit more complicated, because the Toplis model is dependent on melt composition, as well as pressure, temperature, water content, and the olivine forsterite content

  • There are a number of options, firstly, you can specify a fixed olivine Fo content, P and T from input spreadsheet to perform calculations at

[13]:
Eq_Ol_Toplis_FixedFo=pt.calculate_eq_ol_content(liq_comps=myLiquids1, Kd_model="Toplis2005", Fe3Fet_Liq=0.2,
                     P=my_input['P_kbar'], T=my_input['Temperature_C']+273.15, ol_fo=0.8)
Eq_Ol_Toplis_FixedFo
[13]:
Mg#_Liq_Fe2 Mg#_Liq_Fet Kd (Toplis, input Fo) Eq Fo (Toplis, input Fo)
0 0.681666 0.631416 0.343939 0.861610
1 0.620707 0.566947 0.354888 0.821787
2 0.578196 0.523041 0.339995 0.801262
3 0.523445 0.467721 0.343036 0.762017
4 0.586967 0.532031 0.373208 0.792006
5 0.565465 0.510056 0.372049 0.777663
6 0.528731 0.473003 0.365932 0.754055
7 0.620033 0.566244 0.451357 0.783331
8 0.583864 0.528846 0.434877 0.763388
9 0.665729 0.614386 0.336976 0.855286
  • Can do the same, but input olivine compositions, so the function uses the Fo content of each row to perform calculations (rather than a fixed Fo content).

[14]:
Eq_Ol_Toplis_FixedFo2=pt.calculate_eq_ol_content(liq_comps=myLiquids1, ol_comps=myOls1,
        Kd_model="Toplis2005", Fe3Fet_Liq=0.2, P=my_input['P_kbar'], T=my_input['Temperature_C'])
Eq_Ol_Toplis_FixedFo2
[14]:
Mg#_Liq_Fe2 Mg#_Liq_Fet Kd (Toplis, input Fo) Eq Fo (Toplis, input Fo) Fo_meas
0 0.681666 0.631416 0.262240 0.890897 0.872023
1 0.620707 0.566947 0.270643 0.858089 0.903203
2 0.578196 0.523041 0.276728 0.832031 0.835656
3 0.523445 0.467721 0.274170 0.800249 0.863386
4 0.586967 0.532031 0.290580 0.830239 0.903278
5 0.565465 0.510056 0.292062 0.816702 0.890781
6 0.528731 0.473003 0.297262 0.790541 0.838609
7 0.620033 0.566244 0.375747 0.812833 0.801278
8 0.583864 0.528846 0.358296 0.796579 0.817074
9 0.665729 0.614386 0.276607 0.878049 0.822724
  • However, this is clearly a bit of a backwards arguement if olivines and liquids aren’t in equilibrium, because that Fo content isn’t relevant, and you are having to use a Fo content to calculate a Fo content…

  • Instead, if you don’t specify Fo or an olivine content, the function will iterate Kd and olivine Fo content, starting from an olivine Fo content of 0.95 to reach an olivine content in equilibrium with your liquid for the toplis model

[15]:
Eq_Ol_Toplis_IterFo=pt.calculate_eq_ol_content(liq_comps=myLiquids1, Kd_model="Toplis2005", Fe3Fet_Liq=0.2,
                                                P=10, T=my_input['Temperature_C']+273.15)
Eq_Ol_Toplis_IterFo
[15]:
Mg#_Liq_Fe2 Mg#_Liq_Fet Kd (Toplis, Iter) Eq Fo (Toplis, Iter)
0 0.681666 0.631416 0.340829 0.862690
1 0.620707 0.566947 0.358919 0.820127
2 0.578196 0.523041 0.345482 0.798700
3 0.523445 0.467721 0.356696 0.754863
4 0.586967 0.532031 0.381278 0.788460
5 0.565465 0.510056 0.383301 0.772469
6 0.528731 0.473003 0.382311 0.745844
7 0.620033 0.566244 0.463462 0.778805
8 0.583864 0.528846 0.451827 0.756413
9 0.665729 0.614386 0.332066 0.857093
  • You might not know temperature, which is needed for Toplis.

  • But, you have olivine-only thermometers at your fingertips, which can be input into the Eq ol content function.

  • Here we use one of the adapted olivine-liquid thermometers which use calulated DMg from the liquid rather than measured DMg.

[16]:
T_Calc_22=pt.calculate_liq_only_temp(liq_comps=myLiquids1, equationT="T_Put2008_eq22_BeattDMg",
                                     P=my_input['P_kbar'])

Eq_Ol_Toplis_IterFo_calcT=pt.calculate_eq_ol_content(liq_comps=myLiquids1, Kd_model="Toplis2005",
                          Fe3Fet_Liq=0.2, P=my_input['P_kbar'], T=T_Calc_22)

Eq_Ol_Toplis_IterFo_calcT
[16]:
Mg#_Liq_Fe2 Mg#_Liq_Fet Kd (Toplis, Iter) Eq Fo (Toplis, Iter)
0 0.681666 0.631416 0.336288 0.864271
1 0.620707 0.566947 0.332721 0.831037
2 0.578196 0.523041 0.324101 0.808775
3 0.523445 0.467721 0.317484 0.775769
4 0.586967 0.532031 0.354413 0.800390
5 0.565465 0.510056 0.353012 0.786612
6 0.528731 0.473003 0.349961 0.762237
7 0.620033 0.566244 0.404831 0.801226
8 0.583864 0.528846 0.385667 0.784390
9 0.665729 0.614386 0.319173 0.861875

3c - Putirka (2016)

  • Equations 8a-8c are for when the proportion of Fe3 is known

  • Equations 9a-c are for when Fe is entered as Fet

  • Equation 8b and 9b require pressure, so if a presure isn’t entered, you dont get these outputs

  • Similarly, equation 9b requires an estimate of fo2

[17]:
## Here, we dont specify fo2 or P, so only get 8a, 8c, 9a +- errors back
Eq_Ol_Put_noP_fo2=pt.calculate_eq_ol_content(liq_comps=myLiquids1, Kd_model="Putirka2016",
                  Fe3Fet_Liq=0.2)
Eq_Ol_Put_noP_fo2.head()
g:\my drive\postdoc\pymme\mybarometers\thermobar_outer\src\Thermobar\mineral_equilibrium.py:193: UserWarning: Putirka (2016) Kd models equation 8b and 9b are P-dependent you need to enter a P in kbar to get these outputs
  w.warn(
[17]:
Mg#_Liq_Fe2 Mg#_Liq_Fet Eq Fo (Putirka 8a Fe2, Kd=0.33) Eq Fo (Putirka 8a Fe2, Kd=0.33-0.044) Eq Fo (Putirka 8a Fe2, Kd=0.33+0.044) Calc Kd (Putirka 8c, Fe2) Eq Fo (Putirka 8c Fe2) Eq Fo (Putirka 9a Fet, Kd=0.29) Eq Fo (Putirka 9a Fet, Kd=0.29-0.051) Eq Fo (Putirka 9a Fet, Kd=0.29+0.051)
0 0.681666 0.631416 0.866470 0.882176 0.851313 0.346039 0.860883 0.855223 0.877567 0.877567
1 0.620707 0.566947 0.832188 0.851234 0.813975 0.346004 0.825469 0.818658 0.845626 0.845626
2 0.578196 0.523041 0.805970 0.827375 0.785645 0.348050 0.797507 0.790858 0.821056 0.821056
3 0.523445 0.467721 0.768971 0.793411 0.745992 0.349700 0.758510 0.751864 0.786171 0.786171
4 0.586967 0.532031 0.811549 0.832466 0.791657 0.339147 0.807332 0.796761 0.826295 0.826295
[18]:
## Here, we specify fo2 andd P, so only get 8b and 9b outputs as well.
Eq_Ol_Put_P_fo2=pt.calculate_eq_ol_content(liq_comps=myLiquids1, Kd_model="Putirka2016",
                                           Fe3Fet_Liq=0.2, logfo2=-10, P=5)
Eq_Ol_Put_P_fo2.head()
[18]:
Mg#_Liq_Fe2 Mg#_Liq_Fet Eq Fo (Putirka 8a Fe2, Kd=0.33) Eq Fo (Putirka 8a Fe2, Kd=0.33-0.044) Eq Fo (Putirka 8a Fe2, Kd=0.33+0.044) Calc Kd (Putirka 8c, Fe2) Eq Fo (Putirka 8c Fe2) Eq Fo (Putirka 9a Fet, Kd=0.29) Eq Fo (Putirka 9a Fet, Kd=0.29-0.051) Eq Fo (Putirka 9a Fet, Kd=0.29+0.051) Calc Kd (Putirka 8b, Fe2) Eq Fo (Putirka 8b Fe2) Calc Kd (Putirka 9b, Fet) Eq Fo (Putirka 9b Fet)
0 0.681666 0.631416 0.866470 0.882176 0.851313 0.346039 0.860883 0.855223 0.877567 0.877567 0.349228 0.859781 0.308935 0.873919
1 0.620707 0.566947 0.832188 0.851234 0.813975 0.346004 0.825469 0.818658 0.845626 0.845626 0.349509 0.824013 0.307039 0.842019
2 0.578196 0.523041 0.805970 0.827375 0.785645 0.348050 0.797507 0.790858 0.821056 0.821056 0.353321 0.795068 0.305758 0.817624
3 0.523445 0.467721 0.768971 0.793411 0.745992 0.349700 0.758510 0.751864 0.786171 0.786171 0.355555 0.755456 0.308271 0.780850
4 0.586967 0.532031 0.811549 0.832466 0.791657 0.339147 0.807332 0.796761 0.826295 0.826295 0.339543 0.807150 0.299996 0.825696

3d - All Models

  • We can also specify Kd_model=”All” to get results from all of these models

  • Here we are using the temperature from the adapted ol-liq thermometer applied to just liquid compositions

[19]:
Eq_Ol_IterFo_calcT=pt.calculate_eq_ol_content(liq_comps=myLiquids1, Kd_model="All", Fe3Fet_Liq=0.2,
                                                P=my_input['P_kbar'], T=T_Calc_22)
Eq_Ol_IterFo_calcT
[19]:
Mg#_Liq_Fe2 Mg#_Liq_Fet Eq Fo (Roeder, Kd=0.3) Eq Fo (Roeder, Kd=0.33) Eq Fo (Roeder, Kd=0.27) Eq Fo (Matzen, Kd=0.34) Eq Fo (Matzen, Kd=0.352) Eq Fo (Matzen, Kd=0.328) Kd (Toplis, Iter) Eq Fo (Toplis, Iter) Eq Fo (Putirka 8a Fe2, Kd=0.33) Eq Fo (Putirka 8a Fe2, Kd=0.33-0.044) Eq Fo (Putirka 8a Fe2, Kd=0.33+0.044) Calc Kd (Putirka 8c, Fe2) Eq Fo (Putirka 8c Fe2) Eq Fo (Putirka 9a Fet, Kd=0.29) Eq Fo (Putirka 9a Fet, Kd=0.29-0.051) Eq Fo (Putirka 9a Fet, Kd=0.29+0.051) Calc Kd (Putirka 8b, Fe2) Eq Fo (Putirka 8b Fe2)
0 0.681666 0.631416 0.877117 0.866470 0.888030 0.862978 0.858825 0.867172 0.336288 0.864271 0.866470 0.882176 0.851313 0.346039 0.860883 0.855223 0.877567 0.877567 0.346828 0.860610
1 0.620707 0.566947 0.845080 0.832188 0.858378 0.827977 0.822980 0.833035 0.332721 0.831037 0.832188 0.851234 0.813975 0.346004 0.825469 0.818658 0.845626 0.845626 0.347909 0.824677
2 0.578196 0.523041 0.820442 0.805970 0.835443 0.801259 0.795678 0.806919 0.324101 0.808775 0.805970 0.827375 0.785645 0.348050 0.797507 0.790858 0.821056 0.821056 0.353321 0.795068
3 0.523445 0.467721 0.785468 0.768971 0.802688 0.763625 0.757307 0.770049 0.317484 0.775769 0.768971 0.793411 0.745992 0.349700 0.758510 0.751864 0.786171 0.786171 0.355555 0.755456
4 0.586967 0.532031 0.825694 0.811549 0.840342 0.806941 0.801479 0.812477 0.354413 0.800390 0.811549 0.832466 0.791657 0.339147 0.807332 0.796761 0.826295 0.826295 0.339543 0.807150
5 0.565465 0.510056 0.812653 0.797709 0.828169 0.792848 0.787094 0.798688 0.353012 0.786612 0.797709 0.819821 0.776758 0.339433 0.793122 0.782127 0.813288 0.813288 0.340403 0.792654
6 0.528731 0.473003 0.789019 0.772716 0.806025 0.767431 0.761182 0.773782 0.349961 0.762237 0.772716 0.796865 0.749988 0.340052 0.767403 0.755798 0.789713 0.789713 0.341896 0.766436
7 0.620033 0.566244 0.844705 0.831788 0.858030 0.827569 0.822563 0.832637 0.404831 0.801226 0.831788 0.850871 0.813541 0.337757 0.828512 0.818233 0.845252 0.845252 0.336036 0.829236
8 0.583864 0.528846 0.823846 0.809585 0.838619 0.804941 0.799437 0.810521 0.385667 0.784390 0.809585 0.830675 0.789540 0.339023 0.805392 0.794683 0.824452 0.824452 0.338229 0.805759
9 0.665729 0.614386 0.869086 0.857856 0.880615 0.854176 0.849803 0.858595 0.319173 0.861875 0.857856 0.874428 0.841899 0.338564 0.854703 0.846012 0.869560 0.869560 0.338111 0.854869

Example 4 - Rhodes Diagrams

  • We can also calculate equilibrium lines to draw on a Rhodes diagram to assess olivine-liquid equilibrium

  • We don’t have the Toplis model built in, because it depends on liquid composition, so you can’t use it to draw lines for multiple liquid compositions.

  • By default, the user just specifies the max and min liquid Mgno they want to calculate lines for

[20]:
Rhodes=pt.calculate_ol_rhodes_diagram_lines(Min_Mgno=0.5, Max_Mgno=0.7)
Rhodes.head()
[20]:
Mg#_Liq Eq_Ol_Fo_Roeder (Kd=0.3) Eq_Ol_Fo_Roeder (Kd=0.27) Eq_Ol_Fo_Roeder (Kd=0.33) Eq_Ol_Fo_Matzen (Kd=0.34) Eq_Ol_Fo_Matzen (Kd=0.328) Eq_Ol_Fo_Matzen (Kd=0.352)
0 0.500000 0.769231 0.787402 0.751880 0.746269 0.753012 0.739645
1 0.502020 0.770662 0.788751 0.753384 0.747796 0.754512 0.741198
2 0.504040 0.772087 0.790095 0.754883 0.749317 0.756006 0.742745
3 0.506061 0.773506 0.791432 0.756375 0.750832 0.757493 0.744286
4 0.508081 0.774919 0.792763 0.757861 0.752341 0.758975 0.745822
  • We also need to calculate olivine forsterite contents for our measured olivine compositions:

[21]:
Ol_Fo_Calc=pt.calculate_ol_fo(ol_comps=myOls1)
  • And the Mg#s of the matrix glass using the ratio inputted in the spreadsheet (as here, 0.2), or you could also enter Fe3Fet_Liq=value here

[22]:
Liq_Mgno_calc=pt.calculate_liq_mgno(liq_comps=myLiquids1)
  • Here, we use these calculated lines, Mg#s and Ol conetnts to draw min and max equilibrium fields for Roeder and Emslie on the LH axis, and Matzen on the RH axis

[23]:

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12,5)) # Plotting for Roeder and Emslie ax1.set_title('Roeder and Emslie, 1970') # Plotting equilibrium lines ax1.plot(Rhodes['Mg#_Liq'], Rhodes['Eq_Ol_Fo_Roeder (Kd=0.27)'], ':k') ax1.plot(Rhodes['Mg#_Liq'], Rhodes['Eq_Ol_Fo_Roeder (Kd=0.33)'], ':k') ax1.plot(Rhodes['Mg#_Liq'], Rhodes['Eq_Ol_Fo_Roeder (Kd=0.3)'], '-k') # Plotting data ax1.plot(Liq_Mgno_calc, Ol_Fo_Calc, 'ok', mfc='red') ax2.set_title('Matzen (2011)') # Plotting equilibrium lines ax2.plot(Rhodes['Mg#_Liq'], Rhodes['Eq_Ol_Fo_Matzen (Kd=0.328)'], ':k') ax2.plot(Rhodes['Mg#_Liq'], Rhodes['Eq_Ol_Fo_Matzen (Kd=0.352)'], ':k') ax2.plot(Rhodes['Mg#_Liq'], Rhodes['Eq_Ol_Fo_Matzen (Kd=0.34)'], '-k') # Plotting data ax2.plot(Liq_Mgno_calc, Ol_Fo_Calc, 'ok', mfc='red') ax1.set_ylabel('Ol Fo content') ax1.set_xlabel('Glass Mg#') ax2.set_xlabel('Glass Mg#') # Adjust axis limits here ax1.set_ylim([0.74, 0.92]) ax2.set_ylim([0.74, 0.92]) ax1.set_xlim([0.5, 0.7]) ax2.set_xlim([0.5, 0.7])
[23]:
(0.5, 0.7)
../../_images/Examples_Liquid_Ol_Liq_Themometry_Olivine_Liquid_thermometry_49_1.png
  • We can also use the fill_between function to show the equilibrium field or you could show the min and max to provide the full range using the fill between function

[24]:
# Calculate Ol Fo contents to plot
Ol_Fo_Calc=pt.calculate_ol_fo(ol_comps=myOls1)
Liq_Mgno_calc=pt.calculate_liq_mgno(liq_comps=myLiquids1, Fe3Fet_Liq=0.2)
## Here we plot these results
fig, (ax1) = plt.subplots(1, 1, figsize=(6,5))

# This section of code fills between the two lines, where xfill_pap is the sorted x co-ordinate,
# and the y1fill and y2fill are the two y coordinates. You have to have a shared x axis (here, same
# generated set of Mg#s)
xfill_pap = np.sort(Rhodes['Mg#_Liq'])
y1fill_pap = Rhodes['Eq_Ol_Fo_Roeder (Kd=0.27)']
y2fill_pap = Rhodes['Eq_Ol_Fo_Matzen (Kd=0.352)']
ax1.fill_between(xfill_pap, y1fill_pap, y2fill_pap, where=y1fill_pap < y2fill_pap,
                 interpolate=True, color='grey',  alpha=0.3)
ax1.fill_between(xfill_pap, y1fill_pap, y2fill_pap, where=y1fill_pap > y2fill_pap,
                 interpolate=True, color='grey', linewidth=0.5, alpha=0.3)


# # Plotting equilibrium lines
ax1.plot(Rhodes['Mg#_Liq'], Rhodes['Eq_Ol_Fo_Roeder (Kd=0.27)'], '-k', linewidth=1)

ax1.plot(Rhodes['Mg#_Liq'], Rhodes['Eq_Ol_Fo_Matzen (Kd=0.352)'], '-k', linewidth=1)
# Plotting data
ax1.plot(Liq_Mgno_calc, Ol_Fo_Calc, 'ok', mfc='red')
ax1.set_ylabel('Ol Fo content')
ax1.set_xlabel('Glass Mg#')
ax1.set_ylim([0.74, 0.92]);
ax1.set_xlim([0.5, 0.7]);
../../_images/Examples_Liquid_Ol_Liq_Themometry_Olivine_Liquid_thermometry_51_0.png

Example 5 - Calculating Fe3Fet_Liq using a buffer position to assess equilibrium

  • Say we don’t know the Fe3FeT_Liq we want to perform equilibrium calculations at, but we have reason to believe our system is buffered at NNO+1

  • As liquid-olivine temperatures aren’t sensitive to fo2, first we calculate temperature (as buffers are highly temperature-sensitive)

[25]:
eq22_PHinput=pt.calculate_ol_liq_temp(liq_comps=myLiquids1, ol_comps=myOls1,
                                      equationT="T_Put2008_eq22", P=my_input['P_kbar'])
  • Then we use this function to convert NNO+1 to a Fe3/FeT ratio (the function returns a dataframe, and the variable we care about is Fe3Fet_Liq)

[26]:
myLiquids_Fe3_Kress_no_norm_Fe3=pt.convert_fo2_to_fe_partition(liq_comps=myLiquids1,
                                T_K=eq22_PHinput['T_K_calc'], P_kbar=3, fo2="NNO", fo2_offset=1,
                                model="Kress1991", renorm=False).Fe3Fet_Liq
myLiquids_Fe3_Kress_no_norm_Fe3
[26]:
0    0.286549
1    0.294522
2    0.290725
3    0.304642
4    0.288451
5    0.293145
6    0.296381
7    0.276202
8    0.290563
9    0.278905
Name: Fe3Fet_Liq, dtype: float64
  • We then input this into the calculate_ol_liq_temp function as Fe3Fet_Liq=….

[27]:
eq22_PHinput_EqTests=pt.calculate_ol_liq_temp(liq_comps=myLiquids1, ol_comps=myOls1,
                                      equationT="T_Put2008_eq22", P=my_input['P_kbar'], eq_tests=True,
                                    Fe3Fet_Liq=myLiquids_Fe3_Kress_no_norm_Fe3)
eq22_PHinput_EqTests
[27]:
T_K_calc Kd Meas Kd calc (Toplis) ΔKd, Toplis ΔKd, Roeder ΔKd, Matzen SiO2_Liq TiO2_Liq Al2O3_Liq FeOt_Liq ... Cr2O3_Ol NiO_Ol Sample_ID_Ol Fo_meas P2O5_Ol DMg_Meas CNML CSiO2L NF Den_Beat93
0 1289.947705 0.352387 0.325040 0.027347 0.052387 0.012387 57.023602 0.623106 16.332899 4.36174 ... 0 0.0 0 0.872023 0 9.472619 0.170395 0.557471 -0.761843 9.819071
1 1229.813416 0.198881 0.308684 0.109802 0.101119 0.141119 57.658600 0.654150 17.194799 3.90621 ... 0 0.0 1 0.903203 0 14.186789 0.137941 0.570023 -0.816654 10.303635
2 1285.857491 0.304065 0.315891 0.011826 0.004065 0.035935 60.731201 0.862054 17.144199 4.07781 ... 0 0.0 2 0.835656 0 15.531216 0.122184 0.579288 -0.792785 10.250501
3 1198.240159 0.199954 0.296719 0.096765 0.100046 0.140046 61.532799 0.440860 16.508801 3.32990 ... 0 0.0 3 0.863386 0 23.753206 0.097626 0.605753 -0.767076 10.715090
4 1259.174284 0.171088 0.323766 0.152678 0.128912 0.168912 52.969101 0.803412 17.563000 5.93217 ... 0 0.0 4 0.903278 0 10.781292 0.186877 0.523030 -0.842157 10.215321
5 1244.931462 0.180579 0.324986 0.144407 0.119421 0.159421 54.050201 0.857348 17.333300 5.60072 ... 0 0.0 5 0.890781 0 12.276130 0.168070 0.533835 -0.833465 10.295063
6 1241.300118 0.245491 0.333878 0.088387 0.054509 0.094509 55.656300 0.897984 17.117800 5.31307 ... 0 0.0 6 0.838609 0 14.083460 0.145845 0.547483 -0.820914 10.324021
7 1199.601831 0.447306 0.391378 0.055928 0.147306 0.107306 49.054699 0.488832 14.665000 5.83168 ... 0 0.0 7 0.801278 0 7.372298 0.220053 0.554209 -0.789502 9.845140
8 1162.693468 0.354214 0.361834 0.007620 0.054214 0.014214 50.625099 0.334074 16.875000 5.01968 ... 0 0.0 8 0.817074 0 10.509878 0.180293 0.551305 -0.873502 10.229233
9 1339.466029 0.476093 0.299720 0.176373 0.176093 0.136093 51.403301 0.663880 18.019600 5.98440 ... 0 0.0 9 0.822724 0 7.070601 0.228386 0.498172 -0.840523 9.673734

10 rows × 75 columns