Commit 91161f01 authored by Lena Privat Laptop's avatar Lena Privat Laptop
Browse files

Update speciation post-processing routines Octave

parent d563d5c2
......@@ -50,7 +50,8 @@ function run_speciation(fpath=' ',redox="IW")
program_dir = cd(); % store octave path to be able to return
cd('/home/noacklen/Arbeit/SOURCE/CHIC/Octave');
%cd('/home/noacklen/Arbeit/SOURCE/CHIC/Octave');
cd('C:/Users/medhe/GIT/CHIC/Octave');
load('gianluigi_H20_500ppm_CO2_1000ppm.mat')
X_H2O_mantle = 500; %ppm % Needs to fit to speciation file grom Gianluigi or wrong pressures will come out
X_CO2_mantle = 1000; %ppm
......@@ -187,16 +188,16 @@ M_H2O = 2*M_H + M_O;
%X_CO2_mantle = 1000; % ppm
% solubility calculated from Iacono-Marzian et al., 2012; H2O in wt-% and CO2 in wt-ppm
X_H2O = 10 %weight percent of water in the melt
X_CO2 = 10000 %weight ppm of CO2 in the melt
hydrous = 0;
H2O_fl_mole = 0.5; % = P_H2O/P_lid
[H2O_dissolved,CO2_dissolved] = solubility(p_lid*10^-5,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
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)
H2O_fl_mole = dX_H2O_mole / (dX_H2O_mole+dX_CO2_mole)
%X_H2O = 10 %weight percent of water in the melt
%X_CO2 = 10000 %weight ppm of CO2 in the melt
%hydrous = 0;
%H2O_fl_mole = 0.5; % = P_H2O/P_lid
%[H2O_dissolved,CO2_dissolved] = solubility(p_lid*10^-5,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
%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)
%H2O_fl_mole = dX_H2O_mole / (dX_H2O_mole+dX_CO2_mole)
%%%%%%%%%%%%%%%%%%%%
......@@ -349,6 +350,11 @@ for i=2:length(P_H2)
% end
end
P_tot_H2(nt)
P_tot_H2O(nt)
P_tot_CO(nt)
P_tot_CO2(nt)
%f = figure('visible','off');
%plot(time_plot*1e-9,P_tot_H2)
%hold on
......@@ -364,17 +370,20 @@ end
% write all output files in one main folder
path_write = 'C:\OwnCloud_FUB\Literature_Projects\Eigene_Projekte\Project_Ortenzi_Dorn\OctaveResultsNew2\';
%path_write = sprintf('/home/noacklen/ownCloud/Literature_Projects/Eigene_Projekte/Project_Ortenzi_Dorn/OctaveResultsNew/%s/',redox);
path_write = sprintf('/home/noacklen/Arbeit/Ortenzi_Dorn_Paper/%s/',redox); % to avoid update conflicts with NextCloud
name_file = substr(fpath,64); % /scratch-2/PublishedProjects/Project_CaroDorn/Project_CaroDorn_
[name1,name2] = strtok(name_file,'/'); % Review /New_Case10_Visc10/New_Run_1.5M_LowAcc/RC0809_FS0.5_MS1.0
[name2,name3] = strtok(name2,'/'); % New_Case10_Visc10 /New_Run_1.5M_LowAcc/RC0809_FS0.5_MS1.0
[name3,name4] = strtok(name3,'/'); % New_Run_1.5M_LowAcc /RC0809_FS0.5_MS1.0
[name4,residual] = strtok(name4,'/'); % RC0809_FS0.5_MS1.0
%path_write = sprintf('/home/noacklen/Arbeit/Ortenzi_Dorn_Paper/%s/',redox); % to avoid update conflicts with NextCloud
%name_file = substr(fpath,64); % /scratch-2/PublishedProjects/Project_CaroDorn/Project_CaroDorn_
%[name1,name2] = strtok(name_file,'/'); % Review /New_Case10_Visc10/New_Run_1.5M_LowAcc/RC0809_FS0.5_MS1.0
%[name2,name3] = strtok(name2,'/'); % New_Case10_Visc10 /New_Run_1.5M_LowAcc/RC0809_FS0.5_MS1.0
%[name3,name4] = strtok(name3,'/'); % New_Run_1.5M_LowAcc /RC0809_FS0.5_MS1.0
%[name4,residual] = strtok(name4,'/'); % RC0809_FS0.5_MS1.0
name2=name2(9:end); % delete "New_Case"
name3=name3(9:end); % delete "New_Run_"
new_name = sprintf('Case%s--%s--%s.txt',name2,name3,name4);
%name2=name2(9:end); % delete "New_Case"
%name3=name3(9:end); % delete "New_Run_"
%new_name = sprintf('Case%s--%s--%s.txt',name2,name3,name4);
name1=' ';,name2=' ';,name3=' ';,name4=' ';,new_name = 'Test_old';
cd(path_write)
......
This diff is collapsed.
......@@ -176,9 +176,11 @@ switch redox_buffer
H2O_wt_percent_QIF=H2O_mole_fraction_QIF.*100;
QIF_ratio_H2_H2O=log10(H2_mole_fraction_QIF./H2O_mole_fraction_QIF);
if (H2O_mole_fraction_QIF>0)
QIF_ratio_H2_H2O=log10(H2_mole_fraction_QIF./H2O_mole_fraction_QIF);
else
QIF_ratio_H2_H2O=0;
endif
% oxygen fugacity at QIF buffer
......@@ -249,8 +251,11 @@ switch redox_buffer
H2O_wt_percent_IW=H2O_mole_fraction_IW.*100;
IW_ratio_H2_H2O=log10(H2_mole_fraction_IW./H2O_mole_fraction_IW);
if (H2O_mole_fraction_IW>0)
IW_ratio_H2_H2O=log10(H2_mole_fraction_IW./H2O_mole_fraction_IW);
else
IW_ratio_H2_H2O=0;
endif
% oxygen fugacity at IW buffer
......@@ -317,8 +322,11 @@ switch redox_buffer
H2O_wt_percent_QFM=H2O_mole_fraction_QFM.*100;
QFM_ratio_H2_H2O=log10(H2_mole_fraction_QFM./H2O_mole_fraction_QFM);
if (H2O_mole_fraction_QFM>0)
QFM_ratio_H2_H2O=log10(H2_mole_fraction_QFM./H2O_mole_fraction_QFM);
else
QFM_ratio_H2_H2O=0;
endif
% oxygen fugacity at QFM buffer
......@@ -385,9 +393,11 @@ switch redox_buffer
H2O_wt_percent_NiNiO=H2O_mole_fraction_NiNiO.*100;
NiNiO_ratio_H2_H2O=log10(H2_mole_fraction_NiNiO./H2O_mole_fraction_NiNiO);
if (H2O_mole_fraction_NiNiO>0)
NiNiO_ratio_H2_H2O=log10(H2_mole_fraction_NiNiO./H2O_mole_fraction_NiNiO);
else
NiNiO_ratio_H2_H2O=0;
endif
% oxygen fugacity at NiNiO buffer
......
......@@ -18,11 +18,11 @@ 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 = 2000; %p_lid;
P_bar = 4000; %p_lid;
T_K = 1200+273.15; %T_lid;
H2O_fl_mole = 0.01; % pre-defined mole fraction of H2O (over H2O+CO2) in gas bubble, calc partial pressures P_CO2/P_H2O from that
hydrous = 0; % unhydrous okay for water contents below ~4 wt-%; in online calculator only anhydrous version used
hydrous = 1; % 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;
......
......@@ -7,67 +7,71 @@ M_H2O = 2*M_H + M_O;
%X_H2O_mantle = 500; % ppm -> wt or mole?
%X_CO2_mantle = 1000; % ppm
p_lid = 10; % in bar
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
hydrous = 0; % determine later depending on X_H2O
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
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 = 25;
i_max = 100;
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)
%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));
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 (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);
H2O_fl_mole = dX_H2O_mole / (dX_H2O_mole+dX_CO2_mole);
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?
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
%% fprintf('%d H-C-mole = %f CO2_diss = %f H2O_diss = %f P_CO2 = %f P_H2O = %f\n',i,...
%% H2O_fl_mole,CO2_dissolved,H2O_dissolved,(1-H2O_fl_mole)*p_lid*10^-5,H2O_fl_mole*p_lid*10^-5)
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)
% H2O_dissolved
% H2O_dissolved_old
% CO2_dissolved
% CO2_dissolved_old
% 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)
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))
%CO2_diss(p_i) = 1-dX_CO2/X_CO2;
%H2O_diss(p_i) = 1-dX_H2O/X_H2O;
%
%end
%
%semilogx(10.^(ln_p),CO2_diss,'r-',10.^(ln_p),H2O_diss,'b-')
\ No newline at end of file
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
semilogx(10.^(ln_p),CO2_diss,'r-',10.^(ln_p),H2O_diss,'b-')
\ No newline at end of file
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