This page was generated from docs/Examples/CO2_in_Melt_Inclusion_Vapour_Bubbles/Example8a_PropagatingCO2Uncertainties.ipynb. Interactive online version: .

Python Notebook Download

Progapagating uncertainty in melt inclusion vapour bubble reconstructions

• This notebook shows how to use Monte Carlo methods to propagate uncertainty in bubble density, bubble volume and melt density to calculate the uncertainty in the equivalent amount of CO2 in the bubble - e.g. how many ppm you add back into the glass

• In this instance, we get uncertainty in bubble density from repeated Raman measurements, uncertainty in melt density from the code DensityX, and estimate the uncertainty in bubble volume from optical measurements as ~30-50%.

## Install DiadFit if you havent already
#!pip install DiadFit --upgrade

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
# The code to do the MC is in DiadFit, make sure you cite!
import DiadFit as pf
pf.__version__

'1.0.0'

Scenario 1 - You have measured x and y using an optical microscope, but you didnt measure Z, so you are assuming it is the average of x and y.

df_noZ=pd.read_excel('Melt_inclusions_all_dimensions.xlsx', sheet_name='noZ')
df_noZ.head()

	SampleID	Melt_x_um	Melt_y_um	VB_x_um	VB_y_um	CO2_dens_gcm3	err_CO2_dens_gcm3	melt_dens_kgm3	err_melt_dens_kgm3	Unnamed: 9	constant	Vol_melt	Vol_VB	Vol%	CO2_in_melt
0	Test_melt_1	50	33	10	15.0	0.030	0.010	2601	130.05	NaN	4.18879	35853.426159	981.747704	2.738226	3.158276
1	Test_melt_2	40	42	5	5.0	0.040	0.020	2603	130.15	NaN	4.18879	36065.483663	65.449847	0.181475	0.278871
2	Test_melt_3	30	50	6	5.0	0.050	0.010	2700	135.00	NaN	4.18879	31415.926536	86.393798	0.275000	0.509259
3	Test_melt_4	20	30	3	2.5	0.036	0.020	2606	130.30	NaN	4.18879	7853.981634	10.799225	0.137500	0.189946
4	Test_melt_5	100	21	11	12.0	0.042	0.015	2704	135.20	NaN	4.18879	66523.224440	794.822941	1.194805	1.855836

• We are assuming the error on the measurements is 1um, but for Z, lets assume the error is half the difference between x and y.

vol_noZ_av, vol_noZ_all, fig=pf.propagate_CO2_in_bubble(sample_ID=df_noZ['SampleID'],  N_dup=1000, vol_perc_bub=None,
CO2_bub_dens_gcm3=df_noZ['CO2_dens_gcm3'], melt_dens_kgm3=df_noZ['melt_dens_kgm3'],
MI_x=df_noZ['Melt_x_um'], MI_y=df_noZ['Melt_y_um'],VB_x=df_noZ['VB_x_um'], VB_y=df_noZ['VB_y_um'],
error_MI_x=1, error_MI_y=1,error_MI_z=0.5*np.abs(df_noZ['Melt_y_um']-df_noZ['Melt_x_um']),error_VB_x=1, error_VB_y=1,
 error_VB_z=np.abs(df_noZ['VB_y_um']-df_noZ['VB_x_um']),
error_type_dimension='Abs', error_dist_dimension='normal',
error_CO2_bub_dens_gcm3=df_noZ['err_CO2_dens_gcm3'], error_type_CO2_bub_dens_gcm3='Abs', error_dist_CO2_bub_dens_gcm3='normal',
error_melt_dens_kgm3=df_noZ['err_melt_dens_kgm3'], error_type_melt_dens_kgm3='Abs', error_dist_melt_dens_kgm3='normal', neg_values=True)

vol_noZ_av.head()

working on sample number 0

	Filename	CO2_eq_in_melt_noMC	mean_MC_CO2_equiv_melt_ppm	med_MC_CO2_equiv_melt_ppm	std_MC_IQR_CO2_equiv_melt_ppm	16th_quantile_melt_ppm	84th_quantile_melt_ppm
0	Test_melt_1	315.827639	342.152023	305.403692	181.101886	151.279370	513.483142
1	Test_melt_2	27.887058	28.974437	26.647298	16.954758	12.480985	46.390500
2	Test_melt_3	50.925926	31.805252	48.978685	23.642208	28.854819	76.139236
3	Test_melt_4	18.994628	19.983244	16.084909	15.292613	5.066951	35.652177
4	Test_melt_5	185.583647	162.166698	163.250032	149.272392	72.941783	371.486566

Example 1 - Run for a single samle

• This function lets you investigate for 1 MI how different values propagate to uncertainty in calculated CO2 contents

• Here, we show a scenario where you dont use the x-y-z dimensions, you just input the bubble vol% and the error on it (as reported in many supporting datasets)

#  Which of your samples do you want to look at? here we look at the first one
sample_i=0
## How many duplicates do you want in the simulation?
N_dup=1000

######################### VAPOUR BUBBLE VOLUMES ###########################
# What is the measured vol% of your bubble
vol_perc_bub=5
# What is your estimated error on this value? Is this an absolute error, a percentage error - do you think the error is distributed normally or uniformly etc
error_vol_perc_bub=40
error_type_vol_perc_bub='Perc'
error_dist_vol_perc_bub='normal'


######################### VAPOUR BUBBLE DENSITIES ###########################
# What is your measured bubble density in g/cm3 by Raman or microthermometry?
CO2_bub_dens_gcm3=0.02
# What is the error on that value? Here we use the absolute error outputted by DiadFit. This is an absolute error and it should follow a normal distribution as its a 1 sigma value
error_CO2_bub_dens_gcm3=0.01
error_type_CO2_bub_dens_gcm3='Abs'
error_dist_CO2_bub_dens_gcm3='normal'

######################### SILICATE MELT DENSITIES ###########################
# Melt density (e.g. from DensityX)
melt_dens_kgm3=2700
# Estimated error on that - here we use 200 kg/m3, its absolute, and distributed normally
error_melt_dens_kgm3=80
error_type_melt_dens_kgm3='Abs'
error_dist_melt_dens_kgm3='normal'

## Now lets calculate the distribution
test_dist=pf.propagate_CO2_in_bubble_ind(sample_i=sample_i,  N_dup=N_dup,
vol_perc_bub=vol_perc_bub,
error_vol_perc_bub=error_vol_perc_bub, error_type_vol_perc_bub=error_type_vol_perc_bub,
error_dist_vol_perc_bub=error_dist_vol_perc_bub,
CO2_bub_dens_gcm3=CO2_bub_dens_gcm3, error_CO2_bub_dens_gcm3=error_CO2_bub_dens_gcm3,
error_type_CO2_bub_dens_gcm3=error_type_CO2_bub_dens_gcm3,
error_dist_CO2_bub_dens_gcm3=error_dist_CO2_bub_dens_gcm3,
melt_dens_kgm3=melt_dens_kgm3, error_melt_dens_kgm3=error_melt_dens_kgm3,
error_type_melt_dens_kgm3=error_type_melt_dens_kgm3,
error_dist_melt_dens_kgm3=error_dist_melt_dens_kgm3, len_loop=1, neg_values=True)

using error on the bubble volume percent, not the entered dimensions, as error_vol_perc_bub was not None

# Now lets look at the simulation for this one sapmle. The column 'CO2_eq_melt_ppm_noMC' is the same answer as above. The column CO2_eq_melt_ppm_MC gives 1 MonteCarlo result.
test_dist

	CO2_eq_melt_ppm_MC	CO2_eq_melt_ppm_noMC	vol_perc_bub_with_noise	CO2_bub_dens_gcm3_with_noise	melt_dens_kgm3_with_noise	vol_perc_bub	CO2_bub_dens_gcm3	error_type_vol_perc_bub	error_dist_Vol	error_CO2_bub_dens_gcm3	error_type_CO2_bub_dens_gcm3	error_dist_CO2_bub_dens_gcm3
0	184.027060	370.37037	4.831051	0.010431	2738.272230	5	0.02	Perc	normal	0.01	Abs	normal
1	331.485215	370.37037	7.000886	0.012853	2714.556235	5	0.02	Perc	normal	0.01	Abs	normal
2	558.692180	370.37037	3.619383	0.039636	2567.751507	5	0.02	Perc	normal	0.01	Abs	normal
3	530.544696	370.37037	4.255368	0.034486	2766.049050	5	0.02	Perc	normal	0.01	Abs	normal
4	909.139425	370.37037	8.238086	0.028807	2610.302722	5	0.02	Perc	normal	0.01	Abs	normal
...	...	...	...	...	...	...	...	...	...	...	...	...
995	166.182144	370.37037	2.466667	0.018581	2757.999512	5	0.02	Perc	normal	0.01	Abs	normal
996	263.174731	370.37037	3.855987	0.018293	2680.219760	5	0.02	Perc	normal	0.01	Abs	normal
997	344.687962	370.37037	8.222711	0.011141	2657.776998	5	0.02	Perc	normal	0.01	Abs	normal
998	158.276474	370.37037	2.400432	0.017365	2633.646470	5	0.02	Perc	normal	0.01	Abs	normal
999	853.307478	370.37037	5.146933	0.044449	2681.073321	5	0.02	Perc	normal	0.01	Abs	normal

1000 rows × 12 columns

# Lets do the basic math just to persuade ourselves this is working
vol_perc_bub=5
melt_dens_gcm3=2.7
bub_dens_gcm3=0.02
CO2_eq_bubble=10**4*(bub_dens_gcm3*vol_perc_bub)/melt_dens_gcm3
CO2_eq_bubble

370.3703703703703

## Lets plot some of these outputs to check they look how we expect
fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(10,8))
ax1.hist(test_dist['vol_perc_bub_with_noise'], bins=50);
ax1.set_xlabel('MC bubble volume (vol%)')
ax1.set_ylabel('# of simulations')

ax2.hist(test_dist['melt_dens_kgm3_with_noise'], bins=50);
ax2.set_xlabel('MC Melt Density (kg/cm3)')

ax3.hist(test_dist['CO2_bub_dens_gcm3_with_noise'], bins=50);
ax3.set_xlabel('MC Bubble Density (g/cm3)')
ax3.set_ylabel('# of simulations')

ax4.hist(test_dist['CO2_eq_melt_ppm_MC'], bins=50, color='red');
ax4.set_xlabel('CO2 equivalent in the melt held in the bubble (ppm)')
# Plotting preferred value - no Monte Carlo (as no negs can skew)
ax4.plot([test_dist['CO2_eq_melt_ppm_noMC'].iloc[0],
          test_dist['CO2_eq_melt_ppm_noMC'].iloc[0]], [0, N_dup/15], '-k', label='Preferred value');

# Plot standard deviation of simulation
ax4.plot([test_dist['CO2_eq_melt_ppm_noMC'].iloc[0]+np.std(test_dist['CO2_eq_melt_ppm_MC']),
          test_dist['CO2_eq_melt_ppm_noMC'].iloc[0]+np.std(test_dist['CO2_eq_melt_ppm_MC'])], [0, N_dup/15], ':k', label='Preferred+value+1s MC');
ax4.plot([test_dist['CO2_eq_melt_ppm_noMC'].iloc[0]-np.std(test_dist['CO2_eq_melt_ppm_MC']),
          test_dist['CO2_eq_melt_ppm_noMC'].iloc[0]-np.std(test_dist['CO2_eq_melt_ppm_MC'])], [0, N_dup/15], ':k', label='Preferred-value+1s MC');
ax4.legend()

<matplotlib.legend.Legend at 0x1404dfdbad0>

Now lets load in some real data

• This data has the values for errors as columns for each melt inclusion (a row)

• The function outputs 2 things, av is just the mean, and std dev of the simulations. test_dist is all the simulations so you can plot them up!

# Download here - https://github.com/PennyWieser/DiadFit/blob/main/docs/Examples/CO2_in_Melt_Inclusion_Vapour_Bubbles/Wieser_2021_Kilauea.xlsx
df=pd.read_excel('Wieser_2021_Kilauea.xlsx', sheet_name='Sheet2')
df.head()

	sample	PEC_amount	Na2O_PEC	Al2O3_PEC	P2O5_PEC	CaO_PEC	K2O_PEC	TiO2_PEC	SiO2_PEC	MgO_PEC	FeOt_PEC	MnO_PEC	H2O_Liq_meas	CO2_Liq_meas	Melt_dens	Melt_dens_err	Vol_%	CO2_dens_gcm3	CO2_dens_gcm3_std_Dev
0	LL8_613b	16.60	2.469	12.900	0.229	10.702	0.349	2.394	49.770	9.328	11.336	0.142	0.240274	28.979265	2.728376	0.081851	2.723415	0.021977	0.000662
1	LL8_615	30.59	2.131	10.928	0.229	9.435	0.334	1.815	49.676	13.612	11.336	0.151	0.238325	37.495154	2.705747	0.081172	4.254140	0.022912	0.006003
2	LL8_617_a	30.02	2.092	11.323	0.291	9.757	0.324	2.009	49.028	13.349	11.338	0.135	0.235590	41.302324	2.716938	0.081508	3.996225	0.029050	0.000521
3	LL8_623_b	26.68	2.103	11.895	0.299	10.022	0.449	2.448	48.418	12.514	11.331	0.177	0.216237	23.837194	2.728426	0.081853	5.110120	0.052225	0.009579
4	LL8_626	27.81	2.181	11.612	0.209	9.789	0.330	2.027	49.215	12.783	11.331	0.172	0.228748	5.660589	2.714561	0.081437	5.645210	0.029996	0.008963

av, test_dist, fig=pf.propagate_CO2_in_bubble(sample_ID=df['sample'],
N_dup=1000, vol_perc_bub=df['Vol_%'],
CO2_bub_dens_gcm3=df['CO2_dens_gcm3'], melt_dens_kgm3=df['Melt_dens']*1000,
error_vol_perc_bub=50, error_type_vol_perc_bub='Perc', error_dist_vol_perc_bub='normal',
error_CO2_bub_dens_gcm3=df['CO2_dens_gcm3_std_Dev'],
error_type_CO2_bub_dens_gcm3='Abs', error_dist_CO2_bub_dens_gcm3='normal',
error_melt_dens_kgm3=df['Melt_dens_err']*1000, error_type_melt_dens_kgm3='Abs', error_dist_melt_dens_kgm3='normal',
 neg_values=True)
av.head()

working on sample number 0
didnt get inside else loop
using error on the bubble volume percent, not the entered dimensions, as error_vol_perc_bub was not None
didnt get inside else loop
using error on the bubble volume percent, not the entered dimensions, as error_vol_perc_bub was not None
didnt get inside else loop
using error on the bubble volume percent, not the entered dimensions, as error_vol_perc_bub was not None
didnt get inside else loop
using error on the bubble volume percent, not the entered dimensions, as error_vol_perc_bub was not None
didnt get inside else loop
using error on the bubble volume percent, not the entered dimensions, as error_vol_perc_bub was not None
didnt get inside else loop
using error on the bubble volume percent, not the entered dimensions, as error_vol_perc_bub was not None
didnt get inside else loop
using error on the bubble volume percent, not the entered dimensions, as error_vol_perc_bub was not None
didnt get inside else loop
using error on the bubble volume percent, not the entered dimensions, as error_vol_perc_bub was not None
didnt get inside else loop
using error on the bubble volume percent, not the entered dimensions, as error_vol_perc_bub was not None
didnt get inside else loop
using error on the bubble volume percent, not the entered dimensions, as error_vol_perc_bub was not None
didnt get inside else loop
using error on the bubble volume percent, not the entered dimensions, as error_vol_perc_bub was not None
didnt get inside else loop
using error on the bubble volume percent, not the entered dimensions, as error_vol_perc_bub was not None
didnt get inside else loop
using error on the bubble volume percent, not the entered dimensions, as error_vol_perc_bub was not None
didnt get inside else loop
using error on the bubble volume percent, not the entered dimensions, as error_vol_perc_bub was not None
didnt get inside else loop
using error on the bubble volume percent, not the entered dimensions, as error_vol_perc_bub was not None
didnt get inside else loop
using error on the bubble volume percent, not the entered dimensions, as error_vol_perc_bub was not None
didnt get inside else loop
using error on the bubble volume percent, not the entered dimensions, as error_vol_perc_bub was not None
didnt get inside else loop
using error on the bubble volume percent, not the entered dimensions, as error_vol_perc_bub was not None
didnt get inside else loop
using error on the bubble volume percent, not the entered dimensions, as error_vol_perc_bub was not None
didnt get inside else loop
using error on the bubble volume percent, not the entered dimensions, as error_vol_perc_bub was not None
working on sample number 20
didnt get inside else loop
using error on the bubble volume percent, not the entered dimensions, as error_vol_perc_bub was not None
didnt get inside else loop
using error on the bubble volume percent, not the entered dimensions, as error_vol_perc_bub was not None
didnt get inside else loop
using error on the bubble volume percent, not the entered dimensions, as error_vol_perc_bub was not None
didnt get inside else loop
using error on the bubble volume percent, not the entered dimensions, as error_vol_perc_bub was not None
didnt get inside else loop
using error on the bubble volume percent, not the entered dimensions, as error_vol_perc_bub was not None
didnt get inside else loop
using error on the bubble volume percent, not the entered dimensions, as error_vol_perc_bub was not None
didnt get inside else loop
using error on the bubble volume percent, not the entered dimensions, as error_vol_perc_bub was not None
didnt get inside else loop
using error on the bubble volume percent, not the entered dimensions, as error_vol_perc_bub was not None
didnt get inside else loop
using error on the bubble volume percent, not the entered dimensions, as error_vol_perc_bub was not None
didnt get inside else loop
using error on the bubble volume percent, not the entered dimensions, as error_vol_perc_bub was not None
didnt get inside else loop
using error on the bubble volume percent, not the entered dimensions, as error_vol_perc_bub was not None
didnt get inside else loop
using error on the bubble volume percent, not the entered dimensions, as error_vol_perc_bub was not None
didnt get inside else loop
using error on the bubble volume percent, not the entered dimensions, as error_vol_perc_bub was not None
didnt get inside else loop
using error on the bubble volume percent, not the entered dimensions, as error_vol_perc_bub was not None
didnt get inside else loop
using error on the bubble volume percent, not the entered dimensions, as error_vol_perc_bub was not None
didnt get inside else loop
using error on the bubble volume percent, not the entered dimensions, as error_vol_perc_bub was not None
didnt get inside else loop
using error on the bubble volume percent, not the entered dimensions, as error_vol_perc_bub was not None
didnt get inside else loop
using error on the bubble volume percent, not the entered dimensions, as error_vol_perc_bub was not None
didnt get inside else loop
using error on the bubble volume percent, not the entered dimensions, as error_vol_perc_bub was not None
didnt get inside else loop
using error on the bubble volume percent, not the entered dimensions, as error_vol_perc_bub was not None
working on sample number 40
didnt get inside else loop
using error on the bubble volume percent, not the entered dimensions, as error_vol_perc_bub was not None
didnt get inside else loop
using error on the bubble volume percent, not the entered dimensions, as error_vol_perc_bub was not None
didnt get inside else loop
using error on the bubble volume percent, not the entered dimensions, as error_vol_perc_bub was not None
didnt get inside else loop
using error on the bubble volume percent, not the entered dimensions, as error_vol_perc_bub was not None
didnt get inside else loop
using error on the bubble volume percent, not the entered dimensions, as error_vol_perc_bub was not None
didnt get inside else loop
using error on the bubble volume percent, not the entered dimensions, as error_vol_perc_bub was not None
didnt get inside else loop
using error on the bubble volume percent, not the entered dimensions, as error_vol_perc_bub was not None
didnt get inside else loop
using error on the bubble volume percent, not the entered dimensions, as error_vol_perc_bub was not None
didnt get inside else loop
using error on the bubble volume percent, not the entered dimensions, as error_vol_perc_bub was not None
didnt get inside else loop
using error on the bubble volume percent, not the entered dimensions, as error_vol_perc_bub was not None
didnt get inside else loop
using error on the bubble volume percent, not the entered dimensions, as error_vol_perc_bub was not None
didnt get inside else loop
using error on the bubble volume percent, not the entered dimensions, as error_vol_perc_bub was not None
didnt get inside else loop
using error on the bubble volume percent, not the entered dimensions, as error_vol_perc_bub was not None
didnt get inside else loop
using error on the bubble volume percent, not the entered dimensions, as error_vol_perc_bub was not None
didnt get inside else loop
using error on the bubble volume percent, not the entered dimensions, as error_vol_perc_bub was not None
didnt get inside else loop
using error on the bubble volume percent, not the entered dimensions, as error_vol_perc_bub was not None

	Filename	CO2_eq_in_melt_noMC	mean_MC_CO2_equiv_melt_ppm	med_MC_CO2_equiv_melt_ppm	std_MC_IQR_CO2_equiv_melt_ppm	16th_quantile_melt_ppm	84th_quantile_melt_ppm
0	LL8_613b	219.365581	217.215534	218.985842	114.200710	101.939183	330.340602
1	LL8_615	360.243929	358.298977	350.438083	195.562328	161.084693	552.209349
2	LL8_617_a	427.279031	427.858356	423.736977	217.283396	210.811546	645.378337
3	LL8_623_b	978.127362	980.211686	969.344386	529.055224	437.304821	1495.415269
4	LL8_626	623.806538	611.552947	553.109461	361.732326	256.610321	980.074974