Commit 45fd4449 authored by Lena Noack's avatar Lena Noack
Browse files

Update solubility etc, add new example 1D input.txt

parent fbf6e86b
......@@ -35,26 +35,26 @@ M_carb = M_C + 3*M_O; %molar mass of carbonate CO_2^3-
%endif
% 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;
%X_Fe2O3 = 0;
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;
X_Fe2O3 = 0;
% 1921 Kilauea tholeiitic basalt, Holloway 1992, in mass fraction:
X_K2O = 0.0051;
X_Na2O = 0.0215;
X_CaO = 0.1093;
X_MgO = 0.1041;
X_FeO = 0.0968;
X_Fe2O3 = 0.0172;
X_Al2O3 = 0.0172;
X_SiO2 = 0.4916;
X_TiO2 = 0.1333;
%X_K2O = 0.0051;
%X_Na2O = 0.0215;
%X_CaO = 0.1093;
%X_MgO = 0.1041;
%X_FeO = 0.0968;
%X_Fe2O3 = 0.0172;
%X_Al2O3 = 0.0172;
%X_SiO2 = 0.4916;
%X_TiO2 = 0.1333;
% basaltic example composition, all in weight fractions
%X_K2O = 0.0772;
......@@ -83,34 +83,34 @@ X_TiO2 = 0.1333;
%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_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_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);
% 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_H2O_wt/(2*M_H+M_O);
%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_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_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);
% 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_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;
% sum of all should be 1:
%X_K2O + X_Na2O + X_CaO + X_MgO + X_FeO + X_Fe2O3 + X_Al2O3 + X_SiO2 + X_TiO2;
......
function make_movie(fpath)
[fname, fpath, fltidx] = uigetfile("data_char_val_ts.res");
function make_movie(start=0,fpath=' ')
if fpath==' '
[fname, fpath, fltidx] = uigetfile("data_char_val_ts.res");
endif
filenames = sprintf('%s%s', fpath, '/data_snap*')
files = glob(filenames);
for i=1:numel(files)
[~, name] = fileparts(files{i});
eval(sprintf('plot_snap(30,3000,-1,1,0,"%s", "-ascii");', files{i}));
% eval(sprintf('%s = load("%s", "-ascii");', name, files{i}));
nr=str2num(name(14:end));
if nr>=start
eval(sprintf('plot_snap(30,4500,-1,0,3,"%s", "-ascii");', files{i}));
endif
%% eval(sprintf('%s = load("%s", "-ascii");', name, files{i}));
endfor
%plot_snap(30,4000,-1,0,1,"/scratch-3/Simulations_New/Project_MantleMixing/New4_Tm1800_r50km-dTr100km_drho0/","data_snapshot 334115.res")
\ No newline at end of file
%
% Main function to produce field plots and profiles for one time step
%
function plot_snap(max_depl=0,max_T=0,fullV=-1,part=0,reduc=0,fin=' ')
function plot_snap(max_depl=0,max_T=0,fullV=-1,part=0,reduc=0,fin=' ')
% if reduc=1, then plot only T2D, V2D, Vi2D, D2D, Mradp2D
......@@ -45,7 +45,7 @@ function plot_snap(max_depl=0,max_T=0,fullV=-1,part=0,reduc=0,fin=' ')
% plot fields
if(max_depl==0),max_depl=max(max(data.depl));, endif
min_T=min(min(data.T));
min_T=0;%min(min(data.T));
if(max_T==0),max_T=max(max(data.T));, endif
if (reduc<=2)
......@@ -53,7 +53,9 @@ function plot_snap(max_depl=0,max_T=0,fullV=-1,part=0,reduc=0,fin=' ')
produce_plot('T2D_',"Temperature in K",data.T,input,data,res,ar,50,'jet',min_z=min_T,max_z=max_T,full=fullV)
produce_plot('Vi2D_',"Viscosity in log_{10} Pa\ s",log10(data.visc),input,data,res,ar,50,'jet',min_z=-1,max_z=-1,full=fullV)
produce_plot('V2D_',"Velocity in cm/yr",data.vel,input,data,res,ar,50,'jet',min_z=-1,max_z=-1,full=fullV)
endif
if (reduc==0)
% data.cr_age only surface values, need to be set everywhere where there is crust
% data.oc_vol_tot is total volume per lateral/perpendicular cell, need to get depth of crust
nl = size(data.T,1);
......@@ -81,9 +83,11 @@ function plot_snap(max_depl=0,max_T=0,fullV=-1,part=0,reduc=0,fin=' ')
end
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
endif
if (reduc==0)
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))
......@@ -91,7 +95,9 @@ function plot_snap(max_depl=0,max_T=0,fullV=-1,part=0,reduc=0,fin=' ')
elseif(size(data.output_tracer,1)>1)
produce_plot('O2D_',"Tracer",data.output_tracer,input,data,res,ar,50,'jet',min_z=-1,max_z=-1,full=fullV)
endif
if (reduc==0)
%min_z=-max(abs(min(min(data.vl))),max(max(data.vl)));
%max_z=max(abs(min(min(data.vl))),max(max(data.vl)));
produce_plot('Vlat2D_',"Lateral velocity in cm/yr",data.vl,input,data,res,ar,10,'jet',min_z=-1,max_z=-1,full=fullV)
......@@ -99,8 +105,10 @@ function plot_snap(max_depl=0,max_T=0,fullV=-1,part=0,reduc=0,fin=' ')
%max_z=max(abs(min(min(data.vr))),max(max(data.vr)));
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)
if (reduc==0)
......@@ -111,13 +119,15 @@ function plot_snap(max_depl=0,max_T=0,fullV=-1,part=0,reduc=0,fin=' ')
endif
endif
% 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
% if(size(data.depl,1)>1), produce_plot('D2D_',sprintf(strcat("Depletion at t=",num2str(data.time/data.t_yr/1000000.0,"%04.0f"),"Myr")),data.depl,input,data,res,ar,20,'bgry',min_z=0,max_z=max_depl,full=fullV), endif
% if(size(data.depl,1)>1), produce_plot('D2Db_',sprintf(strcat("Depletion at t=",num2str(data.time/data.t_yr/1000000.0,"%04.0f"),"Myr")),data.depl,input,data,res,ar,20,'hot',min_z=0,max_z=max_depl,full=fullV), endif
if (reduc<=2)
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
if(size(data.depl,1)>1), produce_plot('D2D_',sprintf(strcat("Depletion at t=",num2str(data.time/data.t_yr/1000000.0,"%04.0f"),"Myr")),data.depl,input,data,res,ar,20,'bgry',min_z=0,max_z=max_depl,full=fullV), endif
if(size(data.depl,1)>1), produce_plot('D2Db_',sprintf(strcat("Depletion at t=",num2str(data.time/data.t_yr/1000000.0,"%04.0f"),"Myr")),data.depl,input,data,res,ar,20,'hot',min_z=0,max_z=max_depl,full=fullV), endif
endif
if (reduc==0)
% 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.w,1)>1), produce_plot('W2D_',"Water content in ppm",-data.w,input,data,res,ar,50,'jet',min_z=-100,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.w,1)>1), produce_plot('W2D_',"Water content in ppm",-data.w,input,data,res,ar,50,'jet',min_z=-100,max_z=0,full=fullV), endif %input.H2O), 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
if(size(data.rec,1)>1), produce_plot('Rec2D_',"Recycling",data.rec,input,data,res,ar,50,'jet',min_z=-1,max_z=-1,full=fullV), endif %input.nr_ph
......@@ -171,11 +181,18 @@ function x_profile = plot_x_profile(field,data,input,name,title_plot,orient,dime
x_profile(k,3) = max(field(1:nl,1,k));
if (melt==1)
p = (max(r)-r(k))*input.rho*input.g/10.0^6; %km -> GPa
if (p<12)
DeltaTs = (0.1-input.Fe_Mantle/100.0)*(102.3 + 64.1*p - 3.62*p^2);
else
DeltaTs = (0.1-input.Fe_Mantle/100.0)*360.0;
endif
if (p<input.mz)
m_profile(k,1) = input.m1s0 + input.m1s1*p + input.m1s2*p^2 + input.m1s3*p^3;
m_profile(k,1) = input.m1s0 + input.m1s1*p + input.m1s2*p^2 + input.m1s3*p^3 + DeltaTs;
m_profile(k,2) = input.m1l0 + input.m1l1*p + input.m1l2*p^2 + input.m1l3*p^3;
else
m_profile(k,1) = input.m2s0 + input.m2s1*p + input.m2s2*p^2 + input.m2s3*p^3 + input.m2s4*p^4;
m_profile(k,1) = input.m2s0 + input.m2s1*p + input.m2s2*p^2 + input.m2s3*p^3 + input.m2s4*p^4 + DeltaTs;
m_profile(k,2) = input.m2l0 + input.m2l1*p + input.m2l2*p^2 + input.m2l3*p^3 + input.m2l4*p^4;
endif
endif
......@@ -319,7 +336,7 @@ function produce_plot(name,title_plot,field,input,data,res,ar,nr_contour,cmap,mi
else
colormap(cmap)
endif
%caxis([min_z max_z]);
caxis([min_z max_z]);
colorbar
FN = findall(f,'-property','FontName');
......
......@@ -51,6 +51,7 @@ function input = read_input(fpath)
input.Dlid = value(strcmp("s_top", data));
input.DeltaVisc = value(strcmp("DeltaVisc", data));
input.use_melt = value(strcmp("use_melt", data));
input.Fe_Mantle = value(strcmp("Fe_Mantle", data));
if (isempty(input.use_melt)), input.use_melt=0;, endif
if (input.use_melt==1)
input.m1s0 = value(strcmp("melt1s_0", data)); % K
......
......@@ -32,7 +32,7 @@ function read_therm_evol_dim0_new()
t_Gyr_HS = data_HS(:,1)/1e+9;
Rp = 3390; % Mars' radius
Rp = 6371;%Earth %3390; % Mars' radius
t_Gyr = data(:,1)/1e+9;
Dl_km = data(:,11)/1000;
Dcr_km = data(:,10)/1000;
......@@ -74,7 +74,7 @@ function read_therm_evol_dim0_new()
plot(t_Gyr,q_l);
plot(t_Gyr,q_c);
hold off
axis([0,5,0,140])
%axis([0,5,0,140])
xlabel("Time (Gyr)");
ylabel("Heat flux (mW/m^2)");
legend('qs','ql','qc','Location','East')
......@@ -92,7 +92,7 @@ function read_therm_evol_dim0_new()
subplot(3,2,5)
plot(t_Gyr,cr_pr);
axis([0,5,0,60])
%axis([0,5,0,60])
xlabel("Time (Gyr)");
ylabel("Crust producation rate (km^3/yr)");
......@@ -101,6 +101,8 @@ function read_therm_evol_dim0_new()
hold on
plot(t_Gyr,Dcr_km);
index = data(:,19)>0;
% shade(t_Gyr(index),Rp-data(index,19)/1000,t_Gyr(index),Rp-data(index,20)/1000,'FillType',[1 2;2 1])
% fill([t_Gyr(index),fliplr(t_Gyr(index))],[Rp-data(index,19)/1000,fliplr(Rp-data(index,20)/1000)],'r');
plot(t_Gyr(index),Rp-data(index,19)/1000,'k');
plot(t_Gyr(index),Rp-data(index,20)/1000,'k');
hold off
......
......@@ -35,26 +35,26 @@ M_carb = M_C + 3*M_O; %molar mass of carbonate CO_2^3-
%endif
% 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;
%X_Fe2O3 = 0;
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;
X_Fe2O3 = 0;
% 1921 Kilauea tholeiitic basalt, Holloway 1992, in mass fraction:
X_K2O = 0.0051;
X_Na2O = 0.0215;
X_CaO = 0.1093;
X_MgO = 0.1041;
X_FeO = 0.0968;
X_Fe2O3 = 0.0172;
X_Al2O3 = 0.0172;
X_SiO2 = 0.4916;
X_TiO2 = 0.1333;
%X_K2O = 0.0051;
%X_Na2O = 0.0215;
%X_CaO = 0.1093;
%X_MgO = 0.1041;
%X_FeO = 0.0968;
%X_Fe2O3 = 0.0172;
%X_Al2O3 = 0.0172;
%X_SiO2 = 0.4916;
%X_TiO2 = 0.1333;
% basaltic example composition, all in weight fractions
%X_K2O = 0.0772;
......@@ -83,34 +83,34 @@ X_TiO2 = 0.1333;
%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_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_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);
% 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_H2O_wt/(2*M_H+M_O);
%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_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_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);
% 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_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;
% sum of all should be 1:
%X_K2O + X_Na2O + X_CaO + X_MgO + X_FeO + X_Fe2O3 + X_Al2O3 + X_SiO2 + X_TiO2;
......@@ -118,7 +118,7 @@ X_TiO2 = X_TiO2/Sum_mol/(M_Ti+2*M_O);
%% 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_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
......@@ -128,7 +128,7 @@ M_melt = (X_K2O*(2*M_K+M_O) + X_Na2O*(2*M_Na+M_O) + X_CaO*(M_Ca+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)
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)
+ 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);
X_AI = X_Al2O3 / (X_CaO+X_K2O+X_Na2O);
......
......@@ -22,7 +22,7 @@ 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 = 1; % unhydrous okay for water contents below ~4 wt-%; in online calculator only anhydrous version used
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;
......@@ -154,10 +154,10 @@ for H2O_fl_mole=0:0.01:1
% endif
H2O_dissolved = exp(ln_H2O); % in wt-%
if (H2O_fl_mole==0), H2O_dissolved = 0;, endif
if (H2O_fl_mole==0), H2O_dissolved = 0;, end
CO32m_dissolved = exp(ln_CO_3_2min); % in wt-ppm
CO2_dissolved = CO32m_dissolved/M_carb*M_CO2; % in wt-ppm
if (H2O_fl_mole==1), CO2_dissolved = 0;, endif
if (H2O_fl_mole==1), CO2_dissolved = 0;, end
% fprintf('%f %f %f %f\n',P_bar,H2O_fl_mole,H2O_dissolved,CO2_dissolved)
CO2(i) = CO2_dissolved;
......@@ -165,4 +165,5 @@ for H2O_fl_mole=0:0.01:1
i=i+1;
end
figure
plot(H2O,CO2,'k-')
\ 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.
This diff is collapsed.
This source diff could not be displayed because it is too large. You can view the blob instead.