Add new files

Octave/Plot_mole_frac_depth.m

+
+T = zeros(4,1); % four depths
+p = zeros(4,1); % four depths
+
+T(1) = 1400; %K
+p(1) = 1; %bar
+T(2) = 1500; %K
+p(2) = 4500; %bar -> 15 km depth
+T(3) = 1600; %K
+p(3) = 9000; %bar -> 30 km depth (end of crust)
+T(4) = 1800; %K
+p(4) = 30000; %bar -> 100 km depth (below lithosphere)
+z(1) = 0; %km
+z(2) = 15; %km
+z(3) = 30; %km
+z(4) = 100; %km
+
+H2 = zeros(4,4); % four depths, four redox cases
+H2O = zeros(4,4); % four depths, four redox cases
+CO = zeros(4,4); % four depths, four redox cases
+CO2 = zeros(4,4); % four depths, four redox cases
+
+for i=1:4  
+  for j=1:4 
+    [H2(i,j), H2O(i,j), CO(i,j), CO2(i,j)] = ortenzi_simultaneous_eq_C_O_H(T(i), p(i), 0.5, 0.5, j, 0);
+  end
+end
+
+subplot(1,4,1)
+plot(H2(:,1),z,'b--','LineWidth',2)
+hold on
+plot(H2O(:,1),z,'b','LineWidth',2)
+plot(CO(:,1),z,'r--','LineWidth',2)
+plot(CO2(:,1),z,'r','LineWidth',2)
+plot([0,0.5],[15,15],'k:','LineWidth',1)
+plot([0,0.5],[30,30],'k:','LineWidth',1)
+hold off
+%legend('H2','H2O','CO','CO2','location','northeast')
+xlabel('Mole fraction','FontSize',16)
+ylabel('Depth in km','FontSize',16)
+title('QIF','FontSize',18)
+xlim([0,0.5])
+h = gca;  % Handle to currently active axes
+set(h, 'YDir', 'reverse');
+
+subplot(1,4,2)
+plot(H2(:,2),z,'b--','LineWidth',2)
+hold on
+plot(H2O(:,2),z,'b','LineWidth',2)
+plot(CO(:,2),z,'r--','LineWidth',2)
+plot(CO2(:,2),z,'r','LineWidth',2)
+plot([0,0.5],[15,15],'k:','LineWidth',1)
+plot([0,0.5],[30,30],'k:','LineWidth',1)
+hold off
+%legend('H2','H2O','CO','CO2','location','northeast')
+xlabel('Mole fraction','FontSize',16)
+ylabel('Depth in km','FontSize',16)
+title('IW','FontSize',18)
+xlim([0,0.5])
+h = gca;  % Handle to currently active axes
+set(h, 'YDir', 'reverse');
+
+subplot(1,4,3)
+plot(H2(:,3),z,'b--','LineWidth',2)
+hold on
+plot(H2O(:,3),z,'b','LineWidth',2)
+plot(CO(:,3),z,'r--','LineWidth',2)
+plot(CO2(:,3),z,'r','LineWidth',2)
+plot([0,0.5],[15,15],'k:','LineWidth',1)
+plot([0,0.5],[30,30],'k:','LineWidth',1)
+text(0.15,2,'1400K, 1bar')
+text(0.14,16.5,'1500K, 45kbar')
+text(0.14,31.5,'1600K, 90kbar')
+text(0.12,97.5,'1800K, 300kbar')
+hold off
+%legend('H2','H2O','CO','CO2','location','northeast')
+xlabel('Mole fraction','FontSize',16)
+ylabel('Depth in km','FontSize',16)
+title('NiNiO','FontSize',18)
+xlim([0,0.5])
+h = gca;  % Handle to currently active axes
+set(h, 'YDir', 'reverse');
+
+subplot(1,4,4)
+plot(H2(:,4),z,'b--','LineWidth',2)
+hold on
+plot(H2O(:,4),z,'b','LineWidth',2)
+plot(CO(:,4),z,'r--','LineWidth',2)
+plot(CO2(:,4),z,'r','LineWidth',2)
+plot([0,0.5],[15,15],'k:','LineWidth',1)
+plot([0,0.5],[30,30],'k:','LineWidth',1)
+hold off
+legend('H2','H2O','CO','CO2','location','north')
+xlabel('Mole fraction','FontSize',16)
+ylabel('Depth in km','FontSize',16)
+title('QFM','FontSize',18)
+xlim([0,0.5])
+h = gca;  % Handle to currently active axes
+set(h, 'YDir', 'reverse');
+
+% subplot(2,2,1)
+% plot(p,H2(:,1),'b--','LineWidth',2)
+% hold on
+% plot(p,H2O(:,1),'b','LineWidth',2)
+% plot(p,CO(:,1),'r--','LineWidth',2)
+% plot(p,CO2(:,1),'r','LineWidth',2)
+% hold off
+% %legend('H2','H2O','CO','CO2')
+% xlabel('Pressure in bar','FontSize',16)
+% ylabel('Mole fraction','FontSize',16)
+% title('QIF','FontSize',18)
+% 
+% subplot(2,2,2)
+% plot(p,H2(:,2),'b--','LineWidth',2)
+% hold on
+% plot(p,H2O(:,2),'b','LineWidth',2)
+% plot(p,CO(:,2),'r--','LineWidth',2)
+% plot(p,CO2(:,2),'r','LineWidth',2)
+% hold off
+% %legend('H2','H2O','CO','CO2')
+% xlabel('Pressure in bar','FontSize',16)
+% ylabel('Mole fraction','FontSize',16)
+% title('IW','FontSize',18)
+% 
+% subplot(2,2,3)
+% semilogy(p,H2(:,3),'b--','LineWidth',2)
+% hold on
+% semilogy(p,H2O(:,3),'b','LineWidth',2)
+% semilogy(p,CO(:,3),'r--','LineWidth',2)
+% semilogy(p,CO2(:,3),'r','LineWidth',2)
+% hold off
+% legend('H_2','H_2O','CO','CO_2','location','east','FontSize',12)
+% xlabel('Pressure in bar','FontSize',16)
+% ylabel('Mole fraction','FontSize',16)
+% title('NiNiO','FontSize',18)
+% 
+% subplot(2,2,4)
+% semilogy(p,H2(:,4),'b--','LineWidth',2)
+% hold on
+% semilogy(p,H2O(:,4),'b','LineWidth',2)
+% semilogy(p,CO(:,4),'r--','LineWidth',2)
+% semilogy(p,CO2(:,4),'r','LineWidth',2)
+% hold off
+% %legend('H2','H2O','CO','CO2')
+% xlabel('Pressure in bar','FontSize',16)
+% ylabel('Mole fraction','FontSize',16)
+% title('QFM','FontSize',18)
\ No newline at end of file

Octave/heat_sources_time.m

+t = 8.17e+6; % in yr
+
+% specific heat rates in W/kg
+H0U235 = 5.69e-4;
+H0U238 = 9.46e-5; 
+H0Th232= 2.64e-5;
+H0K40  = 2.92e-5;
+
+% initial concentrations (calc. backward 4.5 Gyr)
+X0U235 = 1.2104634497665052E-008; 
+X0U238 = 4.0495624371087973E-008;
+X0Th232= 9.9340593961753404E-008;
+X0K40  = 3.7250244956595340E-007;
+%X0U238 = 0.9928_dp * pI%C_U * 1.0e-9 ! C in ppb
+%X0U235 = 0.0071_dp * pI%C_U * 1.0e-9 ! C in ppb
+%X0Th232 =            pI%C_Th * 1.0e-9!Th in ppb
+%X0K40 = 0.000128_dp* pI%C_K * 1.0e-6 ! K in ppm
+
+% half-life 
+t12U238 = 4.47e+9; % in yr
+t12U235 = 7.04e+8;
+t12Th232 = 1.40e+10;
+t12K40 = 1.25e+9;
+
+
+HtU238 = H0U238*exp(-log(2.0)*t/t12U238);
+HtU235 = H0U235*exp(-log(2.0)*t/t12U235);
+HtTh232 = H0Th232*exp(-log(2.0)*t/t12Th232);
+HtK40 = H0K40*exp(-log(2.0)*t/t12K40);
+
+radiogenic_nuclei = (X0U238*HtU238 + X0U235*HtU235 + X0Th232*HtTh232 + X0K40*HtK40)*10^12 % pW/kg
+

Octave/read_therm_evol_dim0_new.m

+%
+% 
+%
+function read_therm_evol_dim0_new()
+  
+  [fname, fpath, fltidx] = uigetfile("data_char_val_ts.res"); 
+
+  program_dir = cd();  % store octave path to be able to return
+  
+  cd(fpath);
+
+  %%%%%%%%%%%
+  % profile %
+  %%%%%%%%%%%
+
+       
+  %t [yr]	dt [yt]	Tr [K]	Tcr [K]	Tl [K]	Tm [K]	Tb [K]	Tc [K]	Dreg [m]	Dcr [m]	Dl [m]	...
+  % 1      2        3       4      5       6       7       8       9         10      11
+  %Delta_u [m]	Delta_c [m]	Qm [W/kg]	Qcr [W/kg]	qs [W/m^2]	ql [W/m^2]	qc [W/m^2]	mz_min [m]	...
+  %   12            13      14        15          16          17          18            19
+  %mz_max [m]	XmH2O [ppm]	cr_prod [m^3/s]	ma	dr_th [m]	dr_cf [m]	dr_md [m]	u [m/s]	Ri [m]	XiS
+  %  20         21          22            23    24       25         26      27      28      29
+
+  data = dlmread(fname, "", 0, 0);
+
+  
+  %%%%%%%%%%%%%%%%
+  % heat sources %
+  %%%%%%%%%%%%%%%%
+  
+  data_HS = dlmread("data_heat_sources.res", "", 0, 0);
+  t_Gyr_HS = data_HS(:,1)/1e+9;
+  
+  
+  Rp = 3390; % Mars' radius
+  t_Gyr = data(:,1)/1e+9;
+  Dl_km = data(:,11)/1000;
+  Dcr_km = data(:,10)/1000;
+  v_cm_yr = data(:,27)*100*3600*24*365.25;
+  q_s = data(:,16)*1000;
+  q_l = data(:,17)*1000;
+  q_c = data(:,18)*1000;
+  Qm = data(:,14)*1e8;
+  cr_pr = data(:,22)*10^-9*3600*24*365.25;
+
+  nr = size(q_c,1);
+  
+  f = figure();%'visible','off');  
+  
+  subplot(3,2,1)
+  plot(t_Gyr,Dl_km);
+  hold on
+  plot(t_Gyr,Dcr_km);
+  hold off
+  xlabel("Time (Gyr)");
+  ylabel("Thickness (km)");
+  legend('Dl','Dcr')
+
+  subplot(3,2,2)
+  plot(t_Gyr,data(:,8)); % ,"linewidth",3);
+  hold on
+  plot(t_Gyr,data(:,7));
+  plot(t_Gyr,data(:,6));
+  plot(t_Gyr,data(:,5));
+  plot(t_Gyr,data(:,4));
+  hold off
+  xlabel("Time (Gyr)");
+  ylabel("Temperature (K)");
+  legend('Tc','Tb','Tm','Tl','Tcr')
+
+  subplot(3,2,3)
+  plot(t_Gyr,q_s);
+  hold on
+  plot(t_Gyr,q_l);
+  plot(t_Gyr,q_c);
+  hold off
+  axis([0,5,0,140])
+  xlabel("Time (Gyr)");
+  ylabel("Heat flux (mW/m^2)");
+  legend('qs','ql','qc','Location','East')
+
+  subplot(3,2,4)
+  plot(t_Gyr_HS,data_HS(:,11))
+  hold on
+  plot(t_Gyr_HS,data_HS(:,12))
+  hold off
+%  plot(t_Gyr,data(:,21));
+  xlabel("Time (Gyr)");
+%  ylabel("Residual H_2O concentration mantle");
+  ylabel("Heat production pW/kg");
+  legend('Qm','Qcr')
+  
+  subplot(3,2,5)
+  plot(t_Gyr,cr_pr);
+  axis([0,5,0,60])
+  xlabel("Time (Gyr)");
+  ylabel("Crust producation rate (km^3/yr)");
+  
+  subplot(3,2,6)
+  plot(t_Gyr,Dl_km);
+  hold on
+  plot(t_Gyr,Dcr_km);
+  index = data(:,19)>0;
+  plot(t_Gyr(index),Rp-data(index,19)/1000,'k');
+  plot(t_Gyr(index),Rp-data(index,20)/1000,'k');
+  hold off
+  set (gca(), "ydir", "reverse")
+  xlabel("Time (Gyr)");
+  ylabel("Depth (km)");
+  
+  
+  cd(program_dir)
\ No newline at end of file

Octave/solubility_call.m

+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
+  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
+
+graphics_toolkit gnuplot
+
+f = figure('visible','off');
+
+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)
+
+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;