Commit 50a5804d authored by Lena Noack's avatar Lena Noack
Browse files

Update 1D model and several visualization files

parents 3fbc443d 31360ea3
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.1;%0.4;%0.1; %weight percent of water in the melt
X_CO2 = 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
if (H2O_fl_mole==1), dX_CO2=0;, CO2_dissolved=CO2_dissolved_old;, endif
if (H2O_fl_mole==0), dX_H2O=0;, H2O_dissolved=H2O_dissolved_old;, endif
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?
endif
% 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
figure
semilogx(10.^(ln_p),CO2_diss,'r-',10.^(ln_p),H2O_diss,'b-')
\ No newline at end of file
This diff is collapsed.
This diff is collapsed.
This source diff could not be displayed because it is too large. You can view the blob instead.
......@@ -2,7 +2,7 @@
"cells": [
{
"cell_type": "code",
"execution_count": 3,
"execution_count": 2,
"metadata": {},
"outputs": [
{
......
This diff is collapsed.
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment