solubility_call.m 3.42 KB
Newer Older
Lena Noack's avatar
Lena Noack committed
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
37
38
39
40
41
42
43
44
45
46
47
48
49
M_H = 1.0079;
M_C = 12.0107;
M_O = 15.9994;
M_CO2 = M_C + 2*M_O;
M_H2O = 2*M_H + M_O;

%X_H2O_mantle = 500; % ppm -> wt or mole?
%X_CO2_mantle = 1000; % ppm

p_lid = 100; % in bar
T_lid = 1300 + 273.15; % in K

% solubility calculated from Iacono-Marzian et al., 2012; H2O in wt-% and CO2 in wt-ppm
X_H2O = 0.01;%0.1;%0.4;%0.1; %weight percent of water in the melt
X_CO2 = 10000;%600;%4000; %600; %weight ppm of CO2 in the melt

if (X_H2O<4) % personal communication with Fabrice Gaillard
  hydrous=0;
else
  hydrous=1;
end

H2O_dissolved = 0;
CO2_dissolved = 0;

tol = 10e-10;
i_max = 200;

%fprintf('Start with CO2_melt = %f wt-ppm and H2O_melt = %f wt-%%; p=%f, T=%f\n',...
%        X_CO2,X_H2O,p_lid,T_lid)

ln_p = [-3:0.01:3];
CO2_diss = zeros(size(ln_p));
H2O_diss = zeros(size(ln_p));

for p_i=1:length(ln_p)

p_lid = 10.^(ln_p(p_i));

H2O_fl_mole = 0.5; % = P_H2O/P_lid
no_water = 0; % if water pressure goes to zero, set no_water to one, otherwise wrong oscillating solution
i=1;
fac=1;
do
  H2O_dissolved_old = H2O_dissolved;
  CO2_dissolved_old = CO2_dissolved;
  [H2O_dissolved,CO2_dissolved] = solubility(p_lid,T_lid,H2O_fl_mole,X_H2O,hydrous);
  dX_H2O = max(0,X_H2O-H2O_dissolved); % what would go into gas bubble in wt-%
  dX_CO2 = max(0,X_CO2-CO2_dissolved); % what would go into gas bubble in wt-ppm
50
51
  if (H2O_fl_mole==1), dX_CO2=0;, CO2_dissolved=CO2_dissolved_old;, end
  if (H2O_fl_mole==0), dX_H2O=0;, H2O_dissolved=H2O_dissolved_old;, end
Lena Noack's avatar
Lena Noack committed
52
53
54
55
56
57
58
59
  if (dX_H2O+dX_CO2>0)&&(no_water==0)
    dX_H2O_mole = dX_H2O/M_H2O/10^2 / (dX_H2O/M_H2O/10^2 + dX_CO2/M_CO2/10^6);
    dX_CO2_mole = dX_CO2/M_CO2/10^6 / (dX_H2O/M_H2O/10^2 + dX_CO2/M_CO2/10^6);
    fac=fac*0.9;
    H2O_fl_mole = (1-fac)*H2O_fl_mole+fac*(dX_H2O_mole / (dX_H2O_mole+dX_CO2_mole)); % update only by 10% to avoid oscillations between 0 and 100%
  else
    no_water=1;
    H2O_fl_mole = 0;%0.5*H2O_fl_mole; % neither CO2 nor H2O go out, so set to smaller value to see if at least some CO2 goes out?
60
  end
Lena Noack's avatar
Lena Noack committed
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
%  if (i>90)
%  fprintf('%d H-C-mole = %f CO2_diss = %f H2O_diss = %f CO2_deg = %f%% H2O_deg = %f%%\n',i,...
%    H2O_fl_mole,CO2_dissolved,H2O_dissolved,dX_CO2/X_CO2*100,dX_H2O/X_H2O*100)
%  end
  dX_H2O = max(0,X_H2O-H2O_dissolved);
  dX_CO2 = max(0,X_CO2-CO2_dissolved);
  i=i+1;
until ((abs(H2O_dissolved-H2O_dissolved_old)<tol && abs(CO2_dissolved-CO2_dissolved_old)<tol) || (i>i_max))

%if (i>i_max) 
%  fprintf('Warning, iteration stopped after %d iterations at p=%f bar!\n',i,p_lid)
%endif

CO2_diss(p_i) = 1-dX_CO2/X_CO2;
H2O_diss(p_i) = 1-dX_H2O/X_H2O;

end

graphics_toolkit gnuplot

81
%f = figure('visible','off');
Lena Noack's avatar
Lena Noack committed
82
83
84
85
86
87
88
89
90

semilogx(10.^(ln_p),CO2_diss,'r-',10.^(ln_p),H2O_diss,'b-')
set (findall (gcf (), '-property', 'fontsize'), 'fontsize', 20);
set (findall (gcf (), '-property', 'linewidth'), 'linewidth', 5);
title(sprintf(strcat('Melt (ppm): X_{H2O}=',num2str(X_H2O*10000,"%4.0f"),'; X_{CO2}=',num2str(X_CO2,"%4.0f"))),'FontSize',30)
%title(sprintf(strcat('Melt: X_{H2O}=',num2str(X_H2O,"%4.3f"),'; X_{CO2}=',num2str(X_CO2/10000.0,"%4.3f"),' wt-%%')),'FontSize',24)
xlabel('Total pressure in bar','FontSize',28)
ylabel('Abundance of volatiles dissolved in melt','FontSize',28)

91
92
93
94
95
%filename=sprintf(strcat('C:/OwnCloud_FUB/Literature_Projects/Eigene_Projekte/Project_Ortenzi_Dorn/Solubility/',...
%                 'Solubility_H2O_',num2str(X_H2O*10000,'%4.0f'),'_CO2_',num2str(X_CO2,'%4.0f'),'.png'),f);                                                               
%print(filename,'-dpng','-S800,800');
%
%clear f;