Commit a0f5bff2 authored by Lena Noack's avatar Lena Noack
Browse files

Update sollubility; add code Nick Wogan

parent 2cac3544
......@@ -49,14 +49,16 @@ H2O_fl_mole_last = 0;
%fprintf('%d H-C-mole = %f CO2_diss = %f H2O_diss = %f CO2_deg = %f%% H2O_deg = %f%% nw=%d H-C-mole = %f\n',0,...
% H2O_fl_mole,CO2_dissolved,H2O_dissolved,dX_CO2/X_CO2*100,dX_H2O/X_H2O*100,no_water,0)
dX_H2O = max(0,X_H2O-H2O_dissolved); %dX_H2O_lid = max(0,X_H2O-H2O_dissolved);
dX_CO2 = max(0,X_CO2-CO2_dissolved); %dX_CO2_lid = max(0,X_CO2-CO2_dissolved);
% dX_H2O = max(0,X_H2O-H2O_dissolved); %dX_H2O_lid = max(0,X_H2O-H2O_dissolved);
% dX_CO2 = max(0,X_CO2-CO2_dissolved); %dX_CO2_lid = max(0,X_CO2-CO2_dissolved);
%while (redo>0 && j<outer_steps)
while ((abs(H2O_dissolved-H2O_dissolved_old)>tol || abs(CO2_dissolved-CO2_dissolved_old)>tol) && (i<i_max))
while ((abs(H2O_dissolved-H2O_dissolved_old)>tol && abs(CO2_dissolved-CO2_dissolved_old)>tol) && (i<i_max))
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);
[H2O_dissolved,CO2_dissolved] = solubility(p_lid,T_lid,H2O_fl_mole,min(H2O_dissolved,X_H2O),hydrous);
% [H2O_dissolved,CO2_dissolved] = solubility(p_lid,T_lid,H2O_fl_mole,H2O_dissolved,hydrous);
% [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
if (H2O_fl_mole==1), dX_CO2=0;, CO2_dissolved=X_CO2;, end %CO2_dissolved_old;, end
......@@ -65,11 +67,11 @@ H2O_fl_mole_last = 0;
dX_H2O_mole = dX_H2O/mass_H2O/10^2 / (dX_H2O/mass_H2O/10^2 + dX_CO2/mass_CO2/10^6);
dX_CO2_mole = dX_CO2/mass_CO2/10^6 / (dX_H2O/mass_H2O/10^2 + dX_CO2/mass_CO2/10^6);
fac=fac*fac_mult; % update only by 10% to avoid oscillations between 0 and 100%
if (dX_H2O_mole+dX_CO2_mole==0)
H2O_fl_mole_pred = 0;
else
%if (dX_H2O_mole+dX_CO2_mole==0)
% H2O_fl_mole_pred = 0;
%else
H2O_fl_mole_pred = dX_H2O_mole /(dX_H2O_mole+dX_CO2_mole);
end
%end
H2O_fl_mole = (1-fac)*H2O_fl_mole+fac*(H2O_fl_mole_pred);
else
no_water=1;
......@@ -86,13 +88,52 @@ H2O_fl_mole_last = 0;
% end
%fprintf('%d H-C-mole = %f CO2_diss = %f H2O_diss = %f CO2_deg = %f%% H2O_deg = %f%% nw=%d H-C-mole = %f\n',i,...
%H2O_fl_mole,CO2_dissolved,H2O_dissolved,dX_CO2/X_CO2*100,dX_H2O/X_H2O*100,no_water,dX_H2O_mole /(dX_H2O_mole+dX_CO2_mole))
%H2O_fl_mole,CO2_dissolved,H2O_dissolved,dX_CO2/X_CO2*100,dX_H2O/X_H2O*100,no_water,H2O_fl_mole_pred)
% H2O_fl_mole,CO2_dissolved,H2O_dissolved,dX_CO2/X_CO2*100,dX_H2O/X_H2O*100,no_water,dX_H2O_mole /(dX_H2O_mole+dX_CO2_mole))
dX_H2O = max(0,X_H2O-H2O_dissolved); %dX_H2O_lid = max(0,X_H2O-H2O_dissolved);
dX_CO2 = max(0,X_CO2-CO2_dissolved); %dX_CO2_lid = 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))
end
%H2O_fl_vec = 0:0.01:1;
%for i=1:length(H2O_fl_vec)
% [H2O_dissolved,CO2_dissolved] = solubility(p_lid,T_lid,H2O_fl_vec(i),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
% fprintf('H2O_fl=%5.2f: deg_CO2=%5.1f ppm, rem_CO2=%5.1f ppm, deg_H2O=%5.1f wt-%%, rem_H2O=%5.1f wt-%%\n',
% H2O_fl_vec(i),dX_CO2,CO2_dissolved,dX_H2O,H2O_dissolved)
%end
% bisection method to find right H2O_fl_mole:
a = 0.0;
c = 1.0;
b = (a+c)/2.0; # b is H2O_fl_mole
i=0;
while (c-a>tol && i<i_max)
i=i+1;
% [H2O_dissolved,CO2_dissolved] = solubility(p_lid,T_lid,H2O_fl_mole,H2O_dissolved,hydrous);
[H2O_dissolved,CO2_dissolved] = solubility(p_lid,T_lid,b,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
if (dX_H2O+dX_CO2 == 0)
c=a;
else
dX_H2O_mole = dX_H2O/mass_H2O/10^2 / (dX_H2O/mass_H2O/10^2 + dX_CO2/mass_CO2/10^6);
dX_CO2_mole = dX_CO2/mass_CO2/10^6 / (dX_H2O/mass_H2O/10^2 + dX_CO2/mass_CO2/10^6);
if (b<dX_H2O_mole /(dX_H2O_mole+dX_CO2_mole)) %
a=b;
else
c=b;
end
end
b = (a+c)/2.0; # b is H2O_fl_mole
fprintf("a=%5.5f, b=%5.5f, c=%5.5f, H_deg=%8.6f, C_deg=%6.4f, H_diss=%6.4f (%4.2f%%), C_diss=%5.2f (%4.2f%%)\n",
a,b,c,dX_H2O,dX_CO2,H2O_dissolved,H2O_dissolved/X_H2O*100,...
CO2_dissolved,CO2_dissolved/X_CO2*100)
end
% fprintf('%d H-C-mole = %f CO2_diss = %f H2O_diss = %f CO2_deg = %f%% H2O_deg = %f%% nw=%d H-C-mole = %f\n',i,...
% H2O_fl_mole,CO2_dissolved,H2O_dissolved,dX_CO2/X_CO2*100,dX_H2O/X_H2O*100,no_water,H2O_fl_mole_pred)
%
......
......@@ -4,8 +4,8 @@
function plot_snap(max_depl=30,max_T=0,fullV=-1,part=0,reduc=1,fin=' ',print_plot=' ')
% if reduc=1, then plot only T2D, V2D, Vi2D, D2D, Mradp2D
%addpath("/home/noacklen/Arbeit/SOURCE/CHIC/Octave")
addpath("C:/Arbeit/NextCloud_FUB/CHIC/Code_Backup_March2020/Octave")
addpath("/home/noacklen/Arbeit/SOURCE/CHIC/Octave")
%addpath("C:/Arbeit/NextCloud_FUB/CHIC/Code_Backup_March2020/Octave")
[data,input,program_dir,fname,fpath] = read_snap(fin); % read data from snapshot
% cd(program_dir) % return to octave folder
......@@ -25,8 +25,8 @@ function plot_snap(max_depl=30,max_T=0,fullV=-1,part=0,reduc=1,fin=' ',print_plo
data = get_data(data,input); % get fields for plotting and dimensionalize if needed
cd(fpath);
mkdir "plots";
graphics_toolkit gnuplot
%graphics_toolkit gnuplot
graphics_toolkit("fltk")
%size(data.depl)
......@@ -73,9 +73,9 @@ if print_plot~=' '
else
t_profile = plot_x_profile(data.T,data,input,"PrT","Temperature","southwest"," [K]");
if (reduc==0)
% plot all relevant radial profiles (min/mean/max profiles)
t_profile = plot_x_profile(data.T,data,input,"PrT","Temperature","southwest"," [K]");
if (input.use_melt==1), t_profile = plot_x_profile(data.T,data,input,"PrTM","Temperature and Melting temp.","southwest"," [K]",log=0,melt=1);, endif
v_profile = plot_x_profile(data.vel,data,input,"PrVel","Velocity","northeast"," [cm/yr]");
vi_profile = plot_x_profile(data.visc,data,input,"PrVisc","Viscosity","northeast"," [log_{10} Pa\ s]",log=1);
......@@ -126,7 +126,7 @@ else
produce_plot('Cr2D_',sprintf(strcat("Crust age at t=",num2str(data.time/data.t_yr/1000000.0,"%04.0f"),"Myr")),crust,input,data,res,ar,50,'bgry',min_z=-1,max_z=-1,full=fullV)
endif
if (reduc==0)
if (reduc<=2)
if(part==1)
produce_plot('OP2D_',"Tracer",data.output_tracer,input,data,res,ar,50,'jet',min_z=-1,max_z=-1,full=fullV,partX=A(:,1),partY=A(:,2),partV=A(:,21))
%% produce_plot('PD2D_',"Tracer",data.output_tracer,input,data,res,ar,50,'jet',min_z=-1,max_z=-1,full=fullV,partX=A(:,1),partY=A(:,2),partV=A(:,17))
......@@ -145,8 +145,8 @@ else
produce_plot('Vrad2D_',"Radial velocity in cm/yr",data.vr,input,data,res,ar,10,'jet',min_z=-1,max_z=-1,full=fullV)
if(min(min(data.strR))>0), produce_plot('StR2D_',"Strain rate in log_{10} 1/s",log10(data.strR),input,data,res,ar,50,'autumn',min_z=-1,max_z=-1,full=fullV), endif
if(size(data.c,1)>1), produce_plot('C2D_',"Composition",data.c,input,data,res,ar,50,'jet',min_z=-1,max_z=-1,full=fullV), endif
endif
if(size(data.c,1)>1), produce_plot('C2D_',"Composition",data.c,input,data,res,ar,50,'jet',min_z=-1,max_z=-1,full=fullV), endif
if(size(data.ref_r,1)>1)
......@@ -166,8 +166,8 @@ else
if(size(data.depl,1)>1), produce_plot('D2Dc_',sprintf(strcat("Depletion at t=",num2str(data.time/data.t_yr/1000000.0,"%04.0f"),"Myr")),-data.depl,input,data,res,ar,50,'jet',min_z=-max_depl,max_z=0,full=fullV), endif
endif
if(size(data.w,1)>1), produce_plot('W2D_',"Water content in ppm",-data.w,input,data,res,ar,50,'jet',min_z=-150,max_z=0,full=fullV), endif %input.H2O), endif
if (reduc==0)
if(size(data.w,1)>1), produce_plot('W2D_',"Water content in ppm",-data.w,input,data,res,ar,50,'jet',min_z=-150,max_z=0,full=fullV), endif %input.H2O), endif
if(size(data.DeltaF,1)>1), produce_plot('DF2D_',"Melt fraction",data.DeltaF,input,data,res,ar,20,'jet',min_z=-1,max_z=-1,full=fullV), endif
if(size(data.ref_a,1)>1), produce_plot('Al2D_',"Heat expansion",data.ref_a,input,data,res,ar,50,'jet',min_z=-1,max_z=-1,full=fullV), endif
......@@ -180,7 +180,7 @@ else
if(size(data.entropy,1)>1), produce_plot('En2D_',"Entropy",data.entropy,input,data,res,ar,50,'jet',min_z=-1,max_z=-1,full=fullV), endif
endif
if (reduc<=2)
if (reduc==0)
% lateral temperature variations
% data.T_var = data.T; %initialization size
% %[tprof,tprof] = meshgrid(t_profile(:,2)); % ToDo: speed-up if t_profile matrix set up, and then just subract matrix from matrix (no loop needed)? 2D/3D - how to do that?
......@@ -203,7 +203,7 @@ else
endif
if (reduc<=2)
if (reduc==0)
if(size(data.heat,1)>1), produce_plot('H2D_',"Heat souce redistribution",data.heat,input,data,res,ar,50,'jet',min_z=-1,max_z=-1,full=fullV), endif
endif
......@@ -313,6 +313,7 @@ function x_profile = plot_x_profile(field,data,input,name,title_plot,orient,dime
filename=sprintf(strcat('plots/',name,num2str(data.iter,"%07.0f"),'.png'),f);
% graph=strcat("-S", num2str((ar+0.4)*res),",", num2str(res));
% print(filename,'-deps','-S800,720');%graph);
print(filename,'-dpng','-S800,720');%graph);
clear f;
......
......@@ -24,14 +24,6 @@ M_H2O = 2*M_H + M_O;
M_P = 30.973762;
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;
......@@ -50,7 +42,7 @@ M_carb = M_C + 3*M_O; %molar mass of carbonate CO_2^3-
%X_TiO2 = 0.0167;
%X_P2O5 = 0.0051;
%X_Fe2O3 = 0;
%X_H2O = X_H2O_wt; %molar fraction of water in the melt
%in_mole = 0;
% basaltic example composition, aparently in molar fractions, Saal et al 2012, D-6-1 sample
X_K2O = 0.00023;
......@@ -63,6 +55,7 @@ X_SiO2 = 0.4957;
X_TiO2 = 0.0084;
X_Fe2O3 = 0;
X_P2O5 = 0;
in_mole = 1;
% 1921 Kilauea tholeiitic basalt, Holloway 1992, in mass fraction:
%X_K2O = 0.0051;
......@@ -75,6 +68,7 @@ X_P2O5 = 0;
%X_SiO2 = 0.4916;
%X_TiO2 = 0.1333;
%X_P2O5 = 0;
%in_mole = 0;
% basaltic example composition, all in weight fractions
%X_K2O = 0.0772;
......@@ -87,6 +81,7 @@ X_P2O5 = 0;
%X_TiO2 = 0.0089;
%X_Fe2O3 = 0;
%X_P2O5 = 0;
%in_mole = 0;
%% Dixon 1997, MORB Thol. Basalt, wt-fractions:
%X_K2O = 0.07/100;
......@@ -99,60 +94,49 @@ X_P2O5 = 0;
%X_TiO2 = 0.74/100;
%X_Fe2O3 = 0;
%X_P2O5 = 0;
%%X_K2O = 1 - X_Na2O - X_CaO - X_MgO - X_FeO - X_Fe2O3 - X_Al2O3 - X_SiO2 - X_TiO2;
% sum of all should be 1:
%X_K2O + X_Na2O + X_CaO + X_MgO + X_FeO + X_Fe2O3 + X_Al2O3 + X_SiO2 + X_TiO2
%get molar fractions
%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_Fe2O3/(2*M_Fe+3*M_O) + X_P2O5/(2*M_P+5*M_O) + X_H2O/(2*M_H+M_O);
%
%X_H2O = X_H2O/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_Fe2O3 = X_Fe2O3/Sum_mol/(2*M_Fe+3*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);
%X_P2O5 = X_P2O5/Sum_mol/(2*M_P+5*M_O);
% if already in molar fractions:
Sum_mol = X_K2O + X_Na2O + X_CaO + X_MgO + X_FeO + X_Fe2O3 + X_Al2O3 + X_SiO2 + X_TiO2 + X_P2O5 + 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;
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_Fe2O3 = X_Fe2O3/Sum_mol;
X_Al2O3 = X_Al2O3/Sum_mol;
X_SiO2 = X_SiO2/Sum_mol;
X_TiO2 = X_TiO2/Sum_mol;
X_P2O5 = X_P2O5/Sum_mol;
% sum of all should be 1:
%X_K2O + X_Na2O + X_CaO + X_MgO + X_FeO + X_Fe2O3 + X_Al2O3 + X_SiO2 + X_TiO2;
%% 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_Fe2O3*(2*M_Fe+3*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_MgO*(M_Mg+M_O) + X_FeO*(M_Fe+M_O) + X_Fe2O3*(2*M_Fe+3*M_O)/3 ...
%%% + 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
%
%%fwm = (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_Fe2O3*(2*M_Fe+3*M_O)/3 ...
%% + 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)
%fwm = (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_Fe2O3*(2*M_Fe+3*M_O)/3 ...
% + 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);
%in_mole = 0;
include_H2O = 1;
if include_H2O==1
X_H2O = X_H2O_wt/100; %molar fraction of water in the melt
else
X_H2O = 0.0;
endif
if (in_mole==0) %need to make mole fractions first
%get molar fractions from weight fractions
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_Fe2O3/(2*M_Fe+3*M_O) + X_P2O5/(2*M_P+5*M_O) + X_H2O/(2*M_H+M_O);
X_H2O = X_H2O/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_Fe2O3 = X_Fe2O3/Sum_mol/(2*M_Fe+3*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);
X_P2O5 = X_P2O5/Sum_mol/(2*M_P+5*M_O);
else
% already in molar fractions:
Sum_mol = X_K2O + X_Na2O + X_CaO + X_MgO + X_FeO + X_Fe2O3 + X_Al2O3 + X_SiO2 + X_TiO2 + X_P2O5 + 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;
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_Fe2O3 = X_Fe2O3/Sum_mol;
X_Al2O3 = X_Al2O3/Sum_mol;
X_SiO2 = X_SiO2/Sum_mol;
X_TiO2 = X_TiO2/Sum_mol;
X_P2O5 = X_P2O5/Sum_mol;
end
X_AI = X_Al2O3 / (X_CaO+X_K2O+X_Na2O);
......@@ -202,38 +186,31 @@ P_CO2 = CO2_fl_mole*P_bar; %(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-%
X_H2O_mol = exp(ln_H2O)/100/Sum_mol/(2*M_H+M_O); % in wt-amount
if isinf(ln_H2O) | isnan(ln_H2O)
X_H2O_mol = X_H2O; % everything stays dissolved
H2O_dissolved = X_H2O_wt; % in wt-%
else
X_H2O_mol = exp(ln_H2O)/100/Sum_mol/(2*M_H+M_O); % in mole-amount
H2O_dissolved = max(0,exp(ln_H2O)); % in wt-%
endif
if (iscomplex(H2O_dissolved)), H2O_dissolved=X_H2O_wt;, end
% in formula here Fe2O3 missing since not in original paper
%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
%ln_CO_3_2min = X_H2O_mol*d_H2O + X_AI*d_AI + (X_FeO+X_MgO)*d_FeMg + (X_Na2O+X_K2O)*d_NaK ... % in ppm
ln_CO_3_2min = min(X_H2O,X_H2O_mol)*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;
%end
H2O_dissolved = exp(ln_H2O); % in wt-%
if (H2O_fl_mole==0), H2O_dissolved = 0;, end
CO32m_dissolved = exp(ln_CO_3_2min); % in wt-ppm
%H2O_dissolved = max(0,exp(ln_H2O)) % in wt-%
%if (H2O_fl_mole==0), H2O_dissolved = 0;, end ???
%if (iscomplex(H2O_dissolved)), H2O_dissolved=0, end
CO32m_dissolved = max(0,exp(ln_CO_3_2min)); % in wt-ppm
CO2_dissolved = CO32m_dissolved/M_carb*M_CO2; % in wt-ppm
if (iscomplex(CO2_dissolved)), CO2_dissolved=0, end
if ((H2O_fl_mole==1) || (CO2_fl_mole==0)), CO2_dissolved = 0;, end
%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
......@@ -20,7 +20,7 @@ else
hydrous=1;
end
H2O_dissolved = 0;
H2O_dissolved = 0;%X_H2O;%0;
CO2_dissolved = 0;
tol = 10e-10;
......@@ -46,6 +46,7 @@ CO2_dissolved_old=1.0;
while ((abs(H2O_dissolved-H2O_dissolved_old)>tol && abs(CO2_dissolved-CO2_dissolved_old)>tol) && (i<i_max))
H2O_dissolved_old = H2O_dissolved;
CO2_dissolved_old = CO2_dissolved;
% [H2O_dissolved,CO2_dissolved] = solubility(p_lid,T_lid,H2O_fl_mole,H2O_dissolved,hydrous);
[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
......@@ -76,6 +77,7 @@ end
CO2_diss(p_i) = 1-dX_CO2/X_CO2;
H2O_diss(p_i) = 1-dX_H2O/X_H2O;
%fprintf('I %d: CO2_diss=%8.4f H2O_diss=%8.4f\n',p_i,CO2_diss(p_i),H2O_diss(p_i))
end
%graphics_toolkit gnuplot
......