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

Add Earth profile

parent 51b0d1b6
% solubility calculated from Iacono-Marzian et al., 2012
% valied only for T in 900-1500 °C and p in 0 - 10 000 bar (1 GPa)
function [H2O_dissolved,CO2_dissolved] = solubility(P_bar,T_K,H2O_fl_mole,X_H2O_wt,hydrous) % X_H2O here in wt-fraction, not molar!!!
% molar masses, g/mol
M_H = 1.0079;
M_C = 12.0107;
M_O = 15.9994;
M_Na = 22.9897;
M_Mg = 24.305;
M_Al = 26.9815;
M_Si = 28.0855;
M_K = 39.0983;
M_Ca = 40.078;
M_Ti = 47.867;
M_Fe = 55.845;
M_CO2 = M_C + 2*M_O;
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 = 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;
%else
% hydrous = 1;
%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;
% 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;
% basaltic example composition, all in weight fractions
%X_K2O = 0.0772;
%X_Na2O = 0.0195;
%X_CaO = 0.114;
%X_MgO = 0.0575;
%X_FeO = 0.0782;
%X_Al2O3 = 0.1557;
%X_SiO2 = 0.4989;
%X_TiO2 = 0.0089;
%X_Fe2O3 = 0;
%% Dixon 1997, MORB Thol. Basalt, wt-fractions:
%X_K2O = 0.07/100;
%X_Na2O = 2.13/100;
%X_CaO = 11.7/100;
%X_MgO = 10.2/100;
%X_FeO = 8.86/100;
%X_Al2O3 = 16.4/100;
%X_SiO2 = 49.1/100;
%X_TiO2 = 0.74/100;
%X_Fe2O3 = 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_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);
%
%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;
%% 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);
X_AI = X_Al2O3 / (X_CaO+X_K2O+X_Na2O);
if hydrous==0
% non-bridging oxygen divided by ixygen
NBO = 2 * (X_K2O + X_Na2O + X_CaO + X_MgO + X_FeO - X_Al2O3);
NBO_O = NBO / ( 2 * X_SiO2 + 2 * X_TiO2 + 3 * X_Al2O3 + X_MgO + X_FeO + 3*X_Fe2O3 + X_CaO + X_Na2O + X_K2O );
else
NBO = 2 * (X_H2O + X_K2O + X_Na2O + X_CaO + X_MgO + X_FeO - X_Al2O3);
NBO_O = NBO / ( 2 * X_SiO2 + 2 * X_TiO2 + 3 * X_Al2O3 + X_MgO + X_FeO + 3*X_Fe2O3 + X_CaO + X_Na2O + X_K2O + X_H2O );
end
if hydrous==1
% fit parameters for solubility
d_H2O = -16.4;
d_AI = 4.4;
d_FeMg = -17.1;
d_NaK = 22.8;
a_CO2 = 1.0;
b_CO2 = 17.3;
C_CO2 = 0.12;
B_CO2 = -6.0;
a_H2O = 0.53;
b_H2O = 2.35;
B_H2O = -3.37;
C_H2O = -0.02;
else
d_H2O = 2.3;
d_AI = 3.8;
d_FeMg = -16.3;
d_NaK = 20.1;
a_CO2 = 1.0;
b_CO2 = 15.8;
C_CO2 = 0.14;
B_CO2 = -5.3;
a_H2O = 0.54;
b_H2O = 1.24;
B_H2O = -2.95;
C_H2O = 0.02;
end
% partial pressures from H2O_fl_mole (water mole fraction); P_i = mole_i*m_av*g/A
P_CO2 = (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-%
% 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
+ 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
CO2_dissolved = CO32m_dissolved/M_carb*M_CO2; % in wt-ppm
if (H2O_fl_mole==1), 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
\ No newline at end of file
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
%Legend=cell(6,1); -> call in main matlab window
% solubility calculated from Iacono-Marzian et al., 2012; H2O in wt-% and CO2 in wt-ppm
X_H2O = 10000;% in wt-ppm
X_CO2 = 10000;% in wt-ppm
X_H2O = X_H2O / 10000; % from wt-ppm to wt-%
%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;
H2O_dissolved_old = 1;
CO2_dissolved_old = 1;
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 = [-5:0.01:5];
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;
H2O_dissolved = 0;
CO2_dissolved = 0;
H2O_dissolved_old = 1;
CO2_dissolved_old = 1;
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);
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;, end
if (H2O_fl_mole==0), dX_H2O=0;, H2O_dissolved=H2O_dissolved_old;, end
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?
end
% 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;
end
if (i==i_max)
fprintf('Warning, iteration stopped after %d iterations at p=%f bar!\n',i,p_lid)
end
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');
h = semilogx(10.^(ln_p),CO2_diss,'r-.', 10.^(ln_p),H2O_diss,'b-.');
hold on
set(h(1),'linewidth',3);
set(h(2),'linewidth',3);
set (findall (gcf (), '-property', 'fontsize'), 'fontsize', 14);
%set (findall (gcf (), '-property', 'linewidth'), 'linewidth', 4);
%title(sprintf(strcat('Melt (ppm): X_{H2O}=',num2str(X_H2O*10000,"%4.0f"),'; X_{CO2}=',num2str(X_CO2,"%4.0f"))),'FontSize',20)
%%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',20)
ylabel('Abundance of volatiles dissolved in melt','FontSize',20)
%legend(sprintf(strcat('CO_2: X_{H2O}=',num2str(X_H2O*10000,"%4.0f"),'ppm; X_{CO2}=',num2str(X_CO2,"%4.0f"),' ppm')),...
% sprintf(strcat('H_2O: X_{H2O}=',num2str(X_H2O*10000,"%4.0f"),'ppm; X_{CO2}=',num2str(X_CO2,"%4.0f"),' ppm')),...
% 'location','northwest','fontsize',10)
Legend{1}=sprintf(strcat('CO_2: X_{H2O}=',num2str(X_H2O*10000,"%4.0f"),'ppm; X_{CO2}=',num2str(X_CO2,"%4.0f"),' ppm'));
Legend{2}=sprintf(strcat('H_2O: X_{H2O}=',num2str(X_H2O*10000,"%4.0f"),'ppm; X_{CO2}=',num2str(X_CO2,"%4.0f"),' ppm'));
%legend(Legend,'location','northwest','fontsize',10) -> call in main Matlab window
%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;
! Geometry 0
geometry = 1 ! 0-box, 1-sph. annulus, 2-cylinder
dim = 2 ! 2-2D, 3-3D
resolution = 20000
rmin = 0
rmax = 0
phi_factor = 0.25 ! only used when geometry>0, 1.0: full annulus, 0.5: half annulus, ...
CorrRot = 0
! Outgassing 0
use_atmos = 0 !1 ! 0 -> don't use atmos_outgassing, 2 -> include escape (implemented for Mars)
gas_spec = 0 !1
CO2_ini = 500 ! initial amount of CO2 in the mantle
dH2O = 0.01 ! partition coefficient of H20
dCO2 = 1.0 ! if 1 then no specific partitioning of C into melt considered, CO2 depends only on fugacity
fO2 = 0 !-2 ! correspond to the shift of the value IW of fO2. For 1 we calculate for IW+1
mantle_w_H2O = 0 !100 ! initial amount of H2O in the mantle in ppm
mantle_w_r = 0 !1 ! water exponent, if 0 then no effect of water on viscosity considered
sol_hyd = 0.043 ! hydrous solidus change following Katz et al., 2003:
sol_hyd_exp = 0.75 ! Delta_T_sold(w [wt-ppm]) = sol_hyd * (w[wt-ppm])^sol_hyd_exp = 43 C * (w[wt%])^0.75
dehydr_melt = 1.0 !0.05 ! if depletion reaches this value, then material is completely dehydrated (not used here, instead dH2O)
Chi_extr = 0.1 !0.4 ! 40% extrusive volatile loss (following Grott et al., 2011)
alpha_melt = 3.0e-5
Cp_melt = 1793
density_melt = 3000
! Continents 0
add_crust = 1
crust_ini = 50000 !100000 ! initial basaltic crust in m
add_crust_depth = 200000
add_crust_incr = 1
Crust_Lambda = 1 !2
dU238 = 0.002535
dU235 = 0.002535
dTh232 = 0.002042
dK40 = 0.003658
cont_OC_rho = 2900 ! kg/m^3
cont_OC_rhoRef = 2900
cont_OC_TRef = 200 ! C
cont_OC_eta = 1.0e+25
cont_OC_cp = 790
cont_OC_k = 2.3 ! (k/cp/rho = 1e-6)
cont_OC_h = -1 ! is set via Crust_Lambda and mantle heat sources
cont_OC_w = 0
cont_OC_wr = 1
cont_OC_fr_angl = 0
cont_OC_fr_coh = 2e+7 ! Pa
! Bnd_cond 0
Ttop = 288.0 ! temperature at surface, until now: has to be 0 or dimensional value (e.g. 288 for Earth)
Tbottom = 2500.0 ! temperature at CMB, until now: has to be 1 or dimensional value (e.g. 3900 for Earth)
DeltaTc = 1500.0
noslip_t = 0 ! 0: free-slip, 1: no-slip at surface boundary
noslip_b = 0 ! 0: free-slip, 1: no-slip at CMB boundary
inner_bound = 0 ! 1: isolated
outer_bound = 0 ! 1: isolated
open_b_bound = 0 ! 1: open boundary at bottom
periodic = 0 !1 ! 0: reflective boundary, 1: periodic boundary, is currently tested, normally: 0 for box and 1 for sphere
! Restart_TS 0
use_snap = 0 ! insert number of TS from the snapshot, from which the simulation shall re-start
last_snap = 0 ! set to 1, then the last snapshot is automatically used (use_snap is ignored)
Cdt = 1 ! Courant factor when positive value (e.g. 0.1), Delta criterium factor when negative (e.g. -10.0), dt_ini if 0
dt_ini = 1.0e+3 ! Initial time step or fixed time step if CourantCoeff=0
dt_max = 1.0e+6 ! Maximal time step
dt_min = 1.0e+0 ! Minimal time step
ReduceMaxTS = 0
nbiter = 0 ! max number of time steps, not used if tmax>0
iter_out = 0 ! output every ... time steps
tmax = 4.5e+9 ! maximal time, stop the simulation when reached; if tmax=0 then either nbiter or conv is used to stop the simulation (ini 4.5e+9)
t_output = 1e+8 ! output time
binary = 1
! Initial 0
ampl = 0.01 ! amplitude of initially applied sphercial harmonics
sph_l = 0 ! mode of spherical harmonics, chose +/-1.0 times the aspect ratio for box geometry and number of preferred plumes for polar geometry
sph_y = 0
linear = 7 ! initial temperature profile: 0: T=Tini everywhere; 1: linear, 2: conductive profile, 3: convective profile with TBLs
Tini = 1800 ! initial mantle temperature for profile 3
s_bot = 100000 ! lower thermal boundary layer (TBL) for profile 3
s_top = 100000 ! upper thermal boundary layer (TBL) for profile 3
isotherm = 1800
CutSolidus = 1
! Viscosity 0
n_fac = -1 !-3.5 ! Arrhenius AND FK: non-Newtonian factor (1.0 - Newtonian, >1.0 non-Newtonian); negative values: use prefactors A instead or ref viscosity
p_fac = 2.5 ! grain size exponent
iniStrRate = 0.0 ! Initial Strain rate
StartNN = 0 ! Number of timesteps after which non-Newtonian law is applied
Tref = 0 ! reference temperature for Rayleigh number
pref = 0 ! reference pressure in GPa for Rayleigh number
e_gamT = 0 ! FK Viscosity: exp(gamma_T)
e_gamP = 0 ! FK Viscosity: exp(gamma_p)
rheol_vec = -3 ! take pv/ppv diffusion rheology from Tackley et al., 2013
E = 240000.0 ! Arrhenius: actvation energy
V = 5.0 ! Arrhenius: activation volume
E_dis = 0 ! Arrhenius, non-Newtonian: activation energy
V_dis = 0 ! Arrhenius, non-Newtonian: activation volume
A_dif = 3.7e-19 !1.17e-11 ! Arrhenius, mixed Newtonian: prefactor Newtonian
A_dis = 0 ! Arrhenius, mixed Newtonian: prefactor non-Newtonian: A' = A * eta_ref^n * (kappa/D^2)^(n-1)
T0 = 288.0 ! Arrhenius: surface temperature (Tsurf/DeltaT)
AdVisc = 0 ! adapt viscosity at boundaries, such that viscosity at ri/ro is correct
ViscMin = 1.0e+21
ViscMax = 1.0e+30
DeltaVisc = 1
DeltaVisc_LM = 0.01 ! re-scale Tackley's lower mantle viscosity (in 2013 paper an artifical factor of 100 was multiplied to viscosity law)
! Compress 0
Di = 1
DiT = 1
TinclTref = 1
Compressible = 1 ! 0=BA/EBA, 1=TALA, 2=ALA
Compr_Gr = 1.0
output_compr = 0
Psurf = 1.0e-4 ! in GPa
read_profs = 5
ReadNondimValues = 1
! Melt 0
use_melt = 1
LatHeat = 600000.0 ! use either LatH or Entropy; LatH = Entropy(=300)*DeltaT
use_lat_heat = 5
update_solidus = 1
Crust_Lambda = 1
melt1s_0 = 1409.15
melt1s_1 = 134.2
melt1s_2 = -6.581
melt1s_3 = 0.1054
melt1l_0 = 2035.15
melt1l_1 = 57.46
melt1l_2 = -3.487
melt1l_3 = 0.0769
melt_depth = 7.0 ! in GPa
melt_z_change = 10.0 ! in GPa
melt2s_0 = 1835
melt2s_1 = 36.918
melt2s_2 = -0.065444
melt2s_3 = 0.000076686
melt2s_4 = -3.09272E-08
melt2l_0 = 1980
melt2l_1 = 36.918
melt2l_2 = -0.065444
melt2l_3 = 0.000076686
melt2l_4 = -3.09272E-08
max_depl = 0.3
! MantleRefVal 0 ! Here only used for non-dimsensionalization, are re-set via profiles
rho = 4500 ! Reference density mantle
rho_c = 11000 ! Reference density core
Tm_ref = 0 ! Reference mantle temperature: Density = rho*(1-alpha(T-Tm_ref))
k = 4 ! Thermal conductivity mantle
Cp = 1250 ! Specific heat capacity mantle
Cp_c = 700 ! Specific heat capacity core
g = 9.81 ! Surface gravitational acceleration
alpha = 2e-5 ! Heat expansion coefficient
kappa = 1.0e-6
grain_size = 0.001
! Ra_H0 0
Ra = 0.0 ! Rayleigh number at Tref/zref
eta_ref = 1.0e+21 ! reference viscosity, used for example Rayleigh number
H0 = 1.0 ! 1.0 times C_U/C_Th/C_K heating effect in mantle
lambda = 1.0 ! 1.0 times standard heat source decrease over time
H_from_surf = 0
C_U = 20.3 ! present-day values after Wänke and Dreibus 1994 (is re-calculated to initial values in the code)
C_Th = 79.5
C_K = 240
CoreCooling = 1 ! heat flux at core-mantle boundary
! Chemical 0
B = 1.0 ! Buoyancy ratio for chemical diffusion, needed in momentum equation
drho = 0 !733.33 ! density diff pyrolite to harzburgite; is divided by average mantle density and multiplied to B/(alpha*T)
Le = 0 ! Lewis number for chemical diffusion, if 0, then no double-diffusive solver is used
LeV = 0
comp_init = 0 ! +-1 = linear, +-2 = rho value; negative: add cylinder, positive: add random perturbation
comp_ampl = 0 ! Amplitude of composition perturbation
C_ref = 0 ! reference composition (should correspond to reference density)
! Particles 0
use_part = 1 ! use particles instead of compositional field approach
part_n = 0 ! total number of particles, ignored if part_n_cell>0
part_n_cell = 5 ! number of particles per cell
min_part = 3
max_part = 8
part_depth = 1
restart_p = 1 !0
part_rk = 4
bnd_l = 1
bnd_r = 1
fieldsAll = 1
! Numerics 0
debug = 0 !2 !0 ! set between 1 and 5 to obtain detailed information
max_outer = 50 ! number of max outer iterations (coupled energy-momentum solver), used if e_solver>0
max_inner = 3 ! number of max inner iterations (coupled pressure-velocity solver), used if m_solver>0
rout_solver = 0 ! 0: Pardiso solver, 1: linBCG solver, -1: linBCG solver in thermal and Pardiso in stokes
itmax_bcg = 10000
tol_bcg = 1.0e-6
itol_bcg = 1
relax = 0.95 ! relaxation factor for pressure correction (if m_solver not 1)
relax_v = 0.05
noBouss = 0 ! if set to one then no Boussinesq is used in momentum solver
e_solver = 3 ! -1: no e-solver, 0: direct first-order explicit, 1: second-order iterative explicit, 2: first-oder iterative implicit, 3: second-order iterative implicit
CN_theta = 1 ! if set inbetween [0,1] then e_solver is ignored and generalized Euler (including Crank-Nicolson) is used
UseUpwImpl = 1 !-1
UseUpwExpl = 1
m_solver = 1 ! 1: coupled, 2: p and u-v-w, 3: p, u, v, w, -3: explicit, 0: no m-solver
c_solver = 0
t_method = 1 ! 0=FDM, 1=FVM1 (+FDM for 2nd deriv.), 2=FVM2
s_method = 1 ! 0=FDM, 2=FVM2
no_m_solver = 0
convT = 1e-4 ! criterion for convergence of outer iteration, used if e_solver>0 and number of outer iterations < max_outer
convV = 1e-4 ! criterion for convergence of inner iteration, used if m_solver>0 and number of inner iterations < max_inner
conv = 0 ! stops the simulation if convergence of time steps is reached, not used if tmax>0 or nbiter is reached
test = 0 ! output of pressure, when used IDL keyword /test has to be used, as well
beta = 1.0 ! for flux limiter; 1: minmod, 2: superbee
average = 2 !3 ! averaging scheme to be used for viscosity: 0=old, 1=harmonic, 2=arithmetic, 3=geometric
p_average = 2 ! averaging scheme to be used for particles: 1=harmonic, 2=arithmetic, 3=geometric
p_average2 = 2 ! for density averaging; if 0 then set to p_average
p_penalty = 0 !1.0e-7 ! this value shoud be 1.0e-7 and is divided by the maximal viscosity in the stokes equation
nmesures = 0 ! currently not used... should be included again for the cases when quasi-steady-state is reached to gain average values for the last X timesteps
EOF = 1
This source diff could not be displayed because it is too large. You can view the blob instead.
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