solubility.m 5.71 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
% solubility calculated from Iacono-Marzian et al., 2012
% valied only for T in 900-1500 °C and p in 0 - 10 000 bar (1 GPa)

function [H2O_dissolved,CO2_dissolved] = solubility(P_bar,T_K,H2O_fl_mole,X_H2O_wt,hydrous) % X_H2O here in wt-fraction, not molar!!!

% molar masses, g/mol
M_H = 1.0079;
M_C = 12.0107;
M_O = 15.9994;
M_Na = 22.9897;
M_Mg = 24.305;
M_Al = 26.9815;
M_Si = 28.0855;
M_K = 39.0983;
M_Ca = 40.078;
M_Ti = 47.867;
M_Fe = 55.845;
M_CO2 = M_C + 2*M_O;
M_H2O = 2*M_H + M_O;
M_carb = M_C + 3*M_O; %molar mass of carbonate CO_2^3-

% partial pressures P_CO2 from previous degassing, how to integrate that?
%P_bar = 10000; %p_lid;
%T_K = 1200+273.15; %T_lid;
%H2O_fl_mole = 0.5; % pre-defined mole fraction of H2O (over H2O+CO2) in gas bubble, calc partial pressures P_CO2/P_H2O from that

%X_H2O = %0.01; %molar fraction of water in the melt

%hydrous = 0; % unhydrous okay for water contents below ~4 wt-%; in online calculator only anhydrous version used

%if (H2O_fl_mole<0.5) % when to be anhydrous? for now if more CO2 moles than H2O moles
%  hydrous = 0;
%else
%  hydrous = 1;
%endif

37
38
39
40
41
42
43
44
45
46
47
% basaltic example composition, aparently in molar fractions, Saal et al 2012, D-6-1 sample
X_K2O = 0.00023;
X_Na2O = 0.0194; 
X_CaO = 0.1201;
X_MgO = 0.1134;
X_FeO = 0.08;
X_Al2O3 = 0.164;
X_SiO2 = 0.4957;
X_TiO2 = 0.0084;


48
% basaltic example composition, all in weight fractions
49
50
51
52
53
54
55
56
%X_K2O = 0.0772;
%X_Na2O = 0.0195; 
%X_CaO = 0.114;
%X_MgO = 0.0575;
%X_FeO = 0.0782;
%X_Al2O3 = 0.1557;
%X_SiO2 = 0.4989;
%X_TiO2 = 0.0089;
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72

%% Dixon 1997, MORB Thol. Basalt, wt-fractions:
%X_K2O = 0.07/100;
%X_Na2O = 2.13/100; 
%X_CaO = 11.7/100;
%X_MgO = 10.2/100;
%X_FeO = 8.86/100;
%X_Al2O3 = 16.4/100;
%X_SiO2 = 49.1/100;
%X_TiO2 = 0.74/100;
%%X_K2O = 1 - X_Na2O - X_CaO - X_MgO - X_FeO - X_Al2O3 - X_SiO2 - X_TiO2;

% sum of all should be 1:
%X_K2O + X_Na2O + X_CaO + X_MgO + X_FeO + X_Al2O3 + X_SiO2 + X_TiO2

%%get molar fractions
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
%Sum_mol = X_K2O/(2*M_K+M_O) + X_Na2O/(2*M_Na+M_O) + X_CaO/(M_Ca+M_O) + X_MgO/(M_Mg+M_O) ...
%   + X_FeO/(M_Fe+M_O) + X_Al2O3/(2*M_Al+3*M_O) + X_SiO2/(M_Si+2*M_O) + X_TiO2/(M_Ti+2*M_O) ...
%   + X_H2O_wt/(2*M_H+M_O);

%X_H2O = X_H2O_wt/Sum_mol/(2*M_H+M_O);
%X_K2O = X_K2O/Sum_mol/(2*M_K+M_O);
%X_Na2O = X_Na2O/Sum_mol/(2*M_Na+M_O);
%X_CaO = X_CaO/Sum_mol/(M_Ca+M_O);
%X_MgO = X_MgO/Sum_mol/(M_Mg+M_O);
%X_FeO = X_FeO/Sum_mol/(M_Fe+M_O);
%X_Al2O3 = X_Al2O3/Sum_mol/(2*M_Al+3*M_O);
%X_SiO2 = X_SiO2/Sum_mol/(M_Si+2*M_O);
%X_TiO2 = X_TiO2/Sum_mol/(M_Ti+2*M_O);

% if already in molar fractions:
Sum_mol = X_K2O + X_Na2O + X_CaO + X_MgO + X_FeO + X_Al2O3 + X_SiO2 + X_TiO2 + X_H2O_wt/(2*M_H+M_O);
89
90

X_H2O = X_H2O_wt/Sum_mol/(2*M_H+M_O);
91
92
93
94
95
96
97
98
X_K2O = X_K2O/Sum_mol;
X_Na2O = X_Na2O/Sum_mol;
X_CaO = X_CaO/Sum_mol;
X_MgO = X_MgO/Sum_mol;
X_FeO = X_FeO/Sum_mol;
X_Al2O3 = X_Al2O3/Sum_mol;
X_SiO2 = X_SiO2/Sum_mol;
X_TiO2 = X_TiO2/Sum_mol;
99
100

% sum of all should be 1:
101
X_K2O + X_Na2O + X_CaO + X_MgO + X_FeO + X_Al2O3 + X_SiO2 + X_TiO2;
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164

%% M_melt normalized to one oxygen (???)
%M_melt = (X_K2O*(2*M_K+M_O) + X_Na2O*(2*M_Na+M_O) + X_CaO*(M_Ca+M_O) ...
%       + X_MgO*(M_Mg+M_O) + X_FeO*(M_Fe+M_O) ...
%       + X_Al2O3*(2*M_Al+3*M_O) + X_SiO2*(M_Si+2*M_O) + X_TiO2*(M_Ti+2*M_O))
%%       + X_Al2O3*(2*M_Al+3*M_O)/3 + X_SiO2*(M_Si+2*M_O)/2 + X_TiO2*(M_Ti+2*M_O)/2)
%%M_melt = M_melt / M_O

X_AI = X_Al2O3 / (X_CaO+X_K2O+X_Na2O);


if hydrous==0
  % non-bridging oxygen divided by ixygen
  NBO = 2 * (X_K2O + X_Na2O + X_CaO + X_MgO + X_FeO - X_Al2O3);
  NBO_O = NBO / ( 2 * X_SiO2 + 2 * X_TiO2 + 3 * X_Al2O3 + X_MgO + X_FeO + X_CaO + X_Na2O + X_K2O );
else
  NBO = 2 * (X_H2O + X_K2O + X_Na2O + X_CaO + X_MgO + X_FeO - X_Al2O3);
  NBO_O = NBO / ( 2 * X_SiO2 + 2 * X_TiO2 + 3 * X_Al2O3 + X_MgO + X_FeO + X_CaO + X_Na2O + X_K2O + X_H2O );
end

if hydrous==1
  % fit parameters for solubility
  d_H2O = -16.4;
  d_AI = 4.4;
  d_FeMg = -17.1;
  d_NaK = 22.8;
  a_CO2 = 1.0;
  b_CO2 = 17.3;
  C_CO2 = 0.12;
  B_CO2 = -6.0;

  a_H2O = 0.53;
  b_H2O = 2.35;
  B_H2O = -3.37;
  C_H2O = -0.02;
else
  d_H2O = 2.3;
  d_AI = 3.8;
  d_FeMg = -16.3;
  d_NaK = 20.1;
  a_CO2 = 1.0;
  b_CO2 = 15.8;
  C_CO2 = 0.14;
  B_CO2 = -5.3;

  a_H2O = 0.54;
  b_H2O = 1.24;
  B_H2O = -2.95;
  C_H2O = 0.02;
end

% partial pressures from H2O_fl_mole (water mole fraction); P_i = mole_i*m_av*g/A
P_CO2 = (1.0-H2O_fl_mole)*P_bar;
P_H2O = H2O_fl_mole*P_bar;

ln_H2O = a_H2O*log(P_H2O) + b_H2O*NBO_O + B_H2O + C_H2O*P_bar/T_K;     % in wt-%

if hydrous==1
ln_CO_3_2min = X_H2O*d_H2O + X_AI*d_AI + (X_FeO+X_MgO)*d_FeMg + (X_Na2O+X_K2O)*d_NaK ... % in ppm
             + a_CO2*log(P_CO2) + b_CO2*NBO_O + B_CO2 + C_CO2*P_bar/T_K;
else
ln_CO_3_2min = X_AI*d_AI + (X_FeO+X_MgO)*d_FeMg + (X_Na2O+X_K2O)*d_NaK ... % in ppm
             + a_CO2*log(P_CO2) + b_CO2*NBO_O + B_CO2 + C_CO2*P_bar/T_K;
165
end           
166
167
             
H2O_dissolved = exp(ln_H2O); % in wt-%
168
if (H2O_fl_mole==0), H2O_dissolved = 0;, end
169
170
CO32m_dissolved = exp(ln_CO_3_2min); % in wt-ppm
CO2_dissolved = CO32m_dissolved/M_carb*M_CO2; % in wt-ppm
171
if (H2O_fl_mole==1), CO2_dissolved = 0;, end
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188

%fwm=36.594; % KILAUEA THOLEIITIC BASALT, HOLLOWAY 1992
%%CO2_dissolved = CO32m_dissolved*1.202656 / (1+CO32m_dissolved*0.202656)
%Mfac = M_CO2/fwm;
%CO32m_dissolved = CO32m_dissolved*10^-6; % weigth fraction
%CO2_dissolved = CO32m_dissolved/fwm / ( (1.0-CO32m_dissolved)/M_CO2 + CO32m_dissolved/fwm ) * 10^6 % ppm
%CO2_dissolved = CO32m_dissolved * Mfac / (1 + CO32m_dissolved*(Mfac-1)) * 10^6 % ppm
%fwm = M_melt%/M_O
%CO2_dissolved = CO32m_dissolved/fwm / ( (1.0-CO32m_dissolved)/M_CO2 + CO32m_dissolved/fwm ) * 10^6 % ppm
%CO2_dissolved = CO32m_dissolved * Mfac / (1 + CO32m_dissolved*(Mfac-1)) * 10^6 % ppm
%CO2_dissolved = CO32m_dissolved * M_CO2/M_carb * 10^6 % ppm
%
%CO2_dissolved = 10^4 * 4400*CO32m_dissolved/(36.6-44*CO32m_dissolved)
%
%%X_CO23m_melt = X_CO23m_melt - CO32m_dissolved

end %function