%use sigmaY which is C*QVol/omega/N, this also effects tc
R=8.3 ;%gas constant in J/mol/K
T_room=273;
Tmelt=1270; %melting of quartz in kelvin, or use Tw=1000C as in Rice 2006?
N=6*10^23 ;% avogadro # in 1/mol %quartz
G=31*10^9; %shear modulus in Pa
Cs=3750; %shear wavespeed in m/s
a_therm0=1*10^(-7);%1.2*10^(-6) m2/s is thermal diffusivity for quartz in m2/s from V&S
%from hanley fig 4. from 28*10^-7 to 6*10^-7
%Gibert and Mainprice show that if
%you assume porosity, their fig 8, this reduces diffusivity by 1 OM.
%from vosteen &Schellschmidt fig 5, T=0 alpha=1.2* 10^-6, and this is
%starting point from which we reduce by porosity alpha by 50%
Q_vol=240000; % Rice2001 says 180KJ, Nakatanis says 90, 200, and 500 at different points in his paper. 200-250 is trypical for quartz
Q_GB=270000; %weakening occurs if QGB>QVol*B
%B=0.82*(4+273/T0)/5; % to comply with about 58kJ asemictivation energy as measured in Nakatani and scholz 2004.<1
B0=0.89; % in general should decay with inreasing N, since then first contact is larger.
rho=2650;%
c0=730; %specific heat capacity in J/kg/K
omega_vol=5.0*10^(-29); % activation m^3. for quartzite from Rice 2001. ??
omega_GB=7.8*10^(-29); % activation m^3. for quartzite from Rice 2001. ??
tauY=Q_GB/omega_GB/N; %2011 Putlet
%
% for quartzite molar mass =0.0064 kg/mol. divide by denisty 2650kg/m3,
% 2.415*10^-6 m3/mol. so omega*N=2.415*10^-6 m3/mol; divide by N=6*10^23--> omega 0.4*10^-29 ???
sigma0=10^6; %normal stress nominal
sigma00=2300*sigma0;
kappa=10*sigma0;
Ar0=0.95; %stress on contcats in saturation: Ar*sigmaYduct=A*400MPa. we measure Ar =0.2A
%in our simulation at slowest strain rate in T=300C
% ratio=1.6*10^(-5);
% dhdtmin=10^(-16);
tc0=2;
imax=12;
Vt=0.5*Cs; % this must be large enouhg, larger thna most V slip, or else the physics breaks.
%Vn0=0.001*Vt;% this gives for Vt=625m/s Vn0=0.6m/s, and for T=273K tc0=1.
%tc=b*(d/Vn)*exp[(1-B)Qvol/RT]. w/o the exponentail term according to Putlet 2011.
% with according to nakatani 2004. so for T=400K,
% it should be 50,000 according to Nakatani and scholz 2004
%fig 7. or about 1 according to dietrich healing. for C=1 if t=0.27s ;
%Vn=5*10^(-11)*Cs, it gives weakning w.no peak. if Vn=5*10^(-12)*Cs, tc=2.7s,
%gives peak at around V=10^-4
d0=10*10^(-6);% this is contcat size. to be used in falsh heating theory. proctor uses d=5micron
u=-1;
q=5;
V0=1*10^(-12);
b_room=R*T_room/Q_vol/B0;
Vn0=b_room*d0*exp((1-B0)*Q_vol/R/T_room)/tc0; %use measured tc0,e.g.=1s
for i=1:imax
V(i)=V0*10^(i);
t(i)=d0/V(i);
dc(i)=d0;
end
% alpha=d0^2*(B0*Q_vol/omega_vol/N)/10^7;
for j= 1:60
depth(j)=j/2; % in kms.
T0(j)=285+j/2*25; %going down in depth , every km, till 20kms.
%sigmaN(j)=sigma0*3/2*(27-9)*j/2; %remove hydrostatic pressure, but assume thrust regime
sigmaN(j)=sigma0*29*j/2; %remove hydrostatic pressure, but assume thrust regime , follow zoback and townend eqn 7c
% this is assuming sigma1=4sigma3, as in T&S eqn 8.35, taking sigma3 as
% lithostatic a finally plugging sigma1 and sigma3 into T&S eqn 8,27 to
% get sigmaN. using theta =60.
D_therm(j)=q*(sigmaN(j)/10^6)^u; % Ditoro nature 2010.
B(j)=B0*exp(-0.0006*(T0(j)-T_room)); %2011 Putlet*B* fit from Evans JGR 84, indentation hardness of quartz
% to allow sigmaY to vary with T.
Ar(j)=Ar0*(0.55+0.45*tanh(4*(sigmaN(j)-sigma00)/sigma00));
if B(j)>1
B(j)=1;
end
a0(j)=R*T0(j)/Q_GB;
b0(j)=R*T0(j)/Q_vol/B(j);
tau_c=tauY*(1+a0(j)*log(V/Vt));
sigmaY(j)=B(j)*Q_vol/omega_vol/N;%this is actually sigma0, the stress at
%time of onset .
% sigma_cmin(j)=B(j)*sigma_cmin0;
Tsurf=T0(j)+sigmaN(j)*tau_c.*sqrt(V.*D_therm(j)/pi./a_therm0)/rho/c0/sigmaY(j); %surface T needs tau not tauc
T=Tsurf+tau_c.*sqrt(V*d0/pi./a_therm0)/rho/c0; % contact T , Proctor JGR 2014, eq 2.
for i=1:imax
if T(i)> 0.99*Tmelt
T(i)=Tmelt; %saturate T for T>Tmelt or T_decomposition.
end
end
for k=2:111 %iterate to converege T and tau for obtaining st-st
for i=1:imax
a(i)=R*T(i)/Q_GB;
b(i)=R*T(i)/Q_vol/B(j);
%tc(i)=b(i)*(d0/Vn0)*exp((1-B(j))*Q_vol/R/T(i)); % the critical time at T
tc(i)=b(i)*(dc(i)/Vn0)*exp((1-B0)*Q_vol/R/T(i)); % the critical time at T
%I add that at ductility sigma is saturated A which is sigmaY at this level
% i know this from runs that show me that this is
% the value of sigmac for v low strain rate (10^-12)
sigma_c(i)=max(sigmaY(j)*(1-b(i)*log(1+dc(i)/V(i)/tc(i))),sigmaN(j)/Ar(j));
%sigma_c(i)=max(sigmaY(j)*(1-b(i)*log(1+dc(i)/V(i)/tc(i))),sigma_cmin0);
dc(i)=d0*sqrt(sigmaY(1)/sigma_c(i));
% dc(i)=sqrt(alpha*sigmaN(j)/sigma_c(i));
tau_c(i)=max(0.1,tauY*(1+a(i)*log(V(i)/Vt)));
%
%recalculate T now with corrected thermal diffusivity
a_therm(i)=1*10^(-4)./T(i)-0.5*10^(-7);%fig 4 from hanley
%a_therm(i)=a_therm0;
c(i)=c0*(1.7-200/T(i)); % from fitting fig 4, Vosteen & Schellschmidt
ff=sqrt(V(i)/pi./a_therm(i))/rho/c(i);
Tsurf(i)=T0(j)+sigmaN(j)*tau_c(i)*ff*sqrt(D_therm(j))/sigma_c(i); %surface T needs tau not tauc
T(i)=Tsurf(i)+tau_c(i)*sqrt(dc(i))*ff;
if T(i)> 0.95*Tmelt
T(i)=Tmelt;
tau_c(i)=0.95*tau_c(i);
a_therm(i)=a_therm(i-1);
c(i)=c(i-1);
elseif T(i)<270
T(i)=300;
tau_c(i)=tauY*(1+a(i)*log(V(i)/Vt));
end
end
end
for i=1:imax
if T(i)> 0.95*Tmelt
T(i)=Tmelt;
% a_therm(i)=(350./(100+T(i)))*a_therm0;%fig 4 from hanley
%a_therm(i)=10^(-4)./T(i-1)-0.5*10^(-7);%fig 4 from hanley
a_therm(i)=a_therm(i)*0.7;%fig 4 from hanley
c(i)=c(i-1); % from Vosteen & Schellschmidt
%c(i)=c0;
ff=sqrt(V(i)/pi./a_therm(i))/rho/c(i);
%sigma_c(i)=sigma/Ar;
%sigma_c(i)=sigmaY(j);
sigma_c(i)=sigma_c(i-1);
%sigma_c(i)=min(sigma_c);
%dc(i)=d0*sqrt(sigmaY(j)/sigma_c(i));
%tau_c(i)=(Tmelt-T0(j))/ff/(sqrt(dc(i)));
dc(i)=d0*sqrt(sigmaY(1)/sigma_c(i));
tau_c(i)=(Tmelt-T0(j))/ff/(sqrt(D_therm(j))*sigmaN(j)/sigma_c(i)+sqrt(dc(i)));
end
end
%end
TCon(:,j)=T;
tauCon(:,j)=tau_c;
sigmaCon(:,j)=sigma_c;
friction(:,j)=tau_c./sigma_c;
tau(:,j)=sigmaN(j)*friction(:,j);
power(:,j)=tau(:,j).*V(:);
tcCon(:,j)=tc;
aCon(:,j)=a;
dcCon(:,j)=dc;
bCon(:,j)=b;
a_thermCon(:,j)=a_therm;
end
figure(77) %plots area as function of depth
plot(sigmaN./sigmaCon(1,:),depth,'r',Ar,depth,'b')
ylabel('depth')
xlabel('Ar')
hold on
figure(87) %plots area as function of depth
semilogx(V,sigmaN(12)./sigmaCon(:,12),'k+-',V,sigmaN(18)./sigmaCon(:,18),'g+-',V,sigmaN(22)./sigmaCon(:,22),'c+-',V,sigmaN(26)./sigmaCon(:,26),'b+-',V,sigmaN(30)./sigmaCon(:,30),'k+-',V,sigmaN(32)./sigmaCon(:,32),'r+-')
%semilogx(V,sigmaCon(:,6),'k+-',V,sigmaCon(:,12),'c+-',V,sigmaCon(:,14),'b+-',V,sigmaCon(:,16),'+-y',V,sigmaCon(:,18),'+-g',V,sigmaCon(:,20),'^-k')
ylabel('Ar')
xlabel('V')
legend(int2str(depth(12)),int2str(depth(18)),int2str(depth(22)),int2str(depth(26)),int2str(depth(30)),int2str(depth(32)))
hold on
figure(2) %plots shear strenth as function of depth
semilogx(V,friction(:,12),'*-r',V,friction(:,18),'*-g',V,friction(:,22),'*-c',V,friction(:,26),'*-b',V,friction(:,30),'o-y',V,friction(:,32),'o-k')
ylabel('friction')
xlabel('V')
legend(int2str(depth(12)),int2str(depth(18)),int2str(depth(22)),int2str(depth(26)),int2str(depth(30)),int2str(depth(32)))
hold on
figure(3) %plots shear strenth as function of depth
plot(tau(1,1:32)/10^6,depth(1:32),'r-',tau(2,1:34)/10^6,depth(1:34),'g-',tau(3,1:36)/10^6,depth(1:36),'k-')
xlabel('shear strength (MPa)')
ylabel('depth (Km)')
hold on
% figure(1) %plots shear strenth as function of depth
% plot(sigmaCon(1,1:31), depth(1:31),'r-',sigmaCon(2,1:33),depth(1:33),'g-',sigmaCon(3,1:35),depth(1:35),'k-')
% xlabel('contact stress (MPa)')
% ylabel('depth (Km)')
% hold on
%
%from joual et al
% A=12.4; % to add powerlaw creep also from 'as is ' in table 1 of Jaoul et al 1984.
% n=2.3; % those numbers are converted for correct units.
% E=171544;
% for i=1:3
% eps_dot=10^(-12+i);
% stress_power(i,:)=(eps_dot/A)^(1/n)*exp(E/R/n./T0);
% end
%from Hirth et al 2001
figure(3)
A=10^(-11.2); %
n=4; % those numbers are converted for correct units.
E=135000;
fH20=37*10^6;
for i=1:3
width=100*3.1^(i-1)
slip_rate=V(i)
eps_dot=slip_rate/width %this gives a 100m wide for 0.3, 310 meter for 3, 960 m for 30.
stress_power(i,:)=(eps_dot/A/fH20)^(1/n)*exp(E/R/n./T0);
end
plot(stress_power(1,32:60)*10^2,depth(32:60),'r-',stress_power(2,34:60)*10^2,depth(34:60),'g-',stress_power(3,36:60)*10^2,depth(36:60),'k-')
R=8.3 ;%gas constant in J/mol/K
T_room=273;
Tmelt=1270; %melting of quartz in kelvin, or use Tw=1000C as in Rice 2006?
N=6*10^23 ;% avogadro # in 1/mol %quartz
G=31*10^9; %shear modulus in Pa
Cs=3750; %shear wavespeed in m/s
a_therm0=1*10^(-7);%1.2*10^(-6) m2/s is thermal diffusivity for quartz in m2/s from V&S
%from hanley fig 4. from 28*10^-7 to 6*10^-7
%Gibert and Mainprice show that if
%you assume porosity, their fig 8, this reduces diffusivity by 1 OM.
%from vosteen &Schellschmidt fig 5, T=0 alpha=1.2* 10^-6, and this is
%starting point from which we reduce by porosity alpha by 50%
Q_vol=240000; % Rice2001 says 180KJ, Nakatanis says 90, 200, and 500 at different points in his paper. 200-250 is trypical for quartz
Q_GB=270000; %weakening occurs if QGB>QVol*B
%B=0.82*(4+273/T0)/5; % to comply with about 58kJ asemictivation energy as measured in Nakatani and scholz 2004.<1
B0=0.89; % in general should decay with inreasing N, since then first contact is larger.
rho=2650;%
c0=730; %specific heat capacity in J/kg/K
omega_vol=5.0*10^(-29); % activation m^3. for quartzite from Rice 2001. ??
omega_GB=7.8*10^(-29); % activation m^3. for quartzite from Rice 2001. ??
tauY=Q_GB/omega_GB/N; %2011 Putlet
%
% for quartzite molar mass =0.0064 kg/mol. divide by denisty 2650kg/m3,
% 2.415*10^-6 m3/mol. so omega*N=2.415*10^-6 m3/mol; divide by N=6*10^23--> omega 0.4*10^-29 ???
sigma0=10^6; %normal stress nominal
sigma00=2300*sigma0;
kappa=10*sigma0;
Ar0=0.95; %stress on contcats in saturation: Ar*sigmaYduct=A*400MPa. we measure Ar =0.2A
%in our simulation at slowest strain rate in T=300C
% ratio=1.6*10^(-5);
% dhdtmin=10^(-16);
tc0=2;
imax=12;
Vt=0.5*Cs; % this must be large enouhg, larger thna most V slip, or else the physics breaks.
%Vn0=0.001*Vt;% this gives for Vt=625m/s Vn0=0.6m/s, and for T=273K tc0=1.
%tc=b*(d/Vn)*exp[(1-B)Qvol/RT]. w/o the exponentail term according to Putlet 2011.
% with according to nakatani 2004. so for T=400K,
% it should be 50,000 according to Nakatani and scholz 2004
%fig 7. or about 1 according to dietrich healing. for C=1 if t=0.27s ;
%Vn=5*10^(-11)*Cs, it gives weakning w.no peak. if Vn=5*10^(-12)*Cs, tc=2.7s,
%gives peak at around V=10^-4
d0=10*10^(-6);% this is contcat size. to be used in falsh heating theory. proctor uses d=5micron
u=-1;
q=5;
V0=1*10^(-12);
b_room=R*T_room/Q_vol/B0;
Vn0=b_room*d0*exp((1-B0)*Q_vol/R/T_room)/tc0; %use measured tc0,e.g.=1s
for i=1:imax
V(i)=V0*10^(i);
t(i)=d0/V(i);
dc(i)=d0;
end
% alpha=d0^2*(B0*Q_vol/omega_vol/N)/10^7;
for j= 1:60
depth(j)=j/2; % in kms.
T0(j)=285+j/2*25; %going down in depth , every km, till 20kms.
%sigmaN(j)=sigma0*3/2*(27-9)*j/2; %remove hydrostatic pressure, but assume thrust regime
sigmaN(j)=sigma0*29*j/2; %remove hydrostatic pressure, but assume thrust regime , follow zoback and townend eqn 7c
% this is assuming sigma1=4sigma3, as in T&S eqn 8.35, taking sigma3 as
% lithostatic a finally plugging sigma1 and sigma3 into T&S eqn 8,27 to
% get sigmaN. using theta =60.
D_therm(j)=q*(sigmaN(j)/10^6)^u; % Ditoro nature 2010.
B(j)=B0*exp(-0.0006*(T0(j)-T_room)); %2011 Putlet*B* fit from Evans JGR 84, indentation hardness of quartz
% to allow sigmaY to vary with T.
Ar(j)=Ar0*(0.55+0.45*tanh(4*(sigmaN(j)-sigma00)/sigma00));
if B(j)>1
B(j)=1;
end
a0(j)=R*T0(j)/Q_GB;
b0(j)=R*T0(j)/Q_vol/B(j);
tau_c=tauY*(1+a0(j)*log(V/Vt));
sigmaY(j)=B(j)*Q_vol/omega_vol/N;%this is actually sigma0, the stress at
%time of onset .
% sigma_cmin(j)=B(j)*sigma_cmin0;
Tsurf=T0(j)+sigmaN(j)*tau_c.*sqrt(V.*D_therm(j)/pi./a_therm0)/rho/c0/sigmaY(j); %surface T needs tau not tauc
T=Tsurf+tau_c.*sqrt(V*d0/pi./a_therm0)/rho/c0; % contact T , Proctor JGR 2014, eq 2.
for i=1:imax
if T(i)> 0.99*Tmelt
T(i)=Tmelt; %saturate T for T>Tmelt or T_decomposition.
end
end
for k=2:111 %iterate to converege T and tau for obtaining st-st
for i=1:imax
a(i)=R*T(i)/Q_GB;
b(i)=R*T(i)/Q_vol/B(j);
%tc(i)=b(i)*(d0/Vn0)*exp((1-B(j))*Q_vol/R/T(i)); % the critical time at T
tc(i)=b(i)*(dc(i)/Vn0)*exp((1-B0)*Q_vol/R/T(i)); % the critical time at T
%I add that at ductility sigma is saturated A which is sigmaY at this level
% i know this from runs that show me that this is
% the value of sigmac for v low strain rate (10^-12)
sigma_c(i)=max(sigmaY(j)*(1-b(i)*log(1+dc(i)/V(i)/tc(i))),sigmaN(j)/Ar(j));
%sigma_c(i)=max(sigmaY(j)*(1-b(i)*log(1+dc(i)/V(i)/tc(i))),sigma_cmin0);
dc(i)=d0*sqrt(sigmaY(1)/sigma_c(i));
% dc(i)=sqrt(alpha*sigmaN(j)/sigma_c(i));
tau_c(i)=max(0.1,tauY*(1+a(i)*log(V(i)/Vt)));
%
%recalculate T now with corrected thermal diffusivity
a_therm(i)=1*10^(-4)./T(i)-0.5*10^(-7);%fig 4 from hanley
%a_therm(i)=a_therm0;
c(i)=c0*(1.7-200/T(i)); % from fitting fig 4, Vosteen & Schellschmidt
ff=sqrt(V(i)/pi./a_therm(i))/rho/c(i);
Tsurf(i)=T0(j)+sigmaN(j)*tau_c(i)*ff*sqrt(D_therm(j))/sigma_c(i); %surface T needs tau not tauc
T(i)=Tsurf(i)+tau_c(i)*sqrt(dc(i))*ff;
if T(i)> 0.95*Tmelt
T(i)=Tmelt;
tau_c(i)=0.95*tau_c(i);
a_therm(i)=a_therm(i-1);
c(i)=c(i-1);
elseif T(i)<270
T(i)=300;
tau_c(i)=tauY*(1+a(i)*log(V(i)/Vt));
end
end
end
for i=1:imax
if T(i)> 0.95*Tmelt
T(i)=Tmelt;
% a_therm(i)=(350./(100+T(i)))*a_therm0;%fig 4 from hanley
%a_therm(i)=10^(-4)./T(i-1)-0.5*10^(-7);%fig 4 from hanley
a_therm(i)=a_therm(i)*0.7;%fig 4 from hanley
c(i)=c(i-1); % from Vosteen & Schellschmidt
%c(i)=c0;
ff=sqrt(V(i)/pi./a_therm(i))/rho/c(i);
%sigma_c(i)=sigma/Ar;
%sigma_c(i)=sigmaY(j);
sigma_c(i)=sigma_c(i-1);
%sigma_c(i)=min(sigma_c);
%dc(i)=d0*sqrt(sigmaY(j)/sigma_c(i));
%tau_c(i)=(Tmelt-T0(j))/ff/(sqrt(dc(i)));
dc(i)=d0*sqrt(sigmaY(1)/sigma_c(i));
tau_c(i)=(Tmelt-T0(j))/ff/(sqrt(D_therm(j))*sigmaN(j)/sigma_c(i)+sqrt(dc(i)));
end
end
%end
TCon(:,j)=T;
tauCon(:,j)=tau_c;
sigmaCon(:,j)=sigma_c;
friction(:,j)=tau_c./sigma_c;
tau(:,j)=sigmaN(j)*friction(:,j);
power(:,j)=tau(:,j).*V(:);
tcCon(:,j)=tc;
aCon(:,j)=a;
dcCon(:,j)=dc;
bCon(:,j)=b;
a_thermCon(:,j)=a_therm;
end
figure(77) %plots area as function of depth
plot(sigmaN./sigmaCon(1,:),depth,'r',Ar,depth,'b')
ylabel('depth')
xlabel('Ar')
hold on
figure(87) %plots area as function of depth
semilogx(V,sigmaN(12)./sigmaCon(:,12),'k+-',V,sigmaN(18)./sigmaCon(:,18),'g+-',V,sigmaN(22)./sigmaCon(:,22),'c+-',V,sigmaN(26)./sigmaCon(:,26),'b+-',V,sigmaN(30)./sigmaCon(:,30),'k+-',V,sigmaN(32)./sigmaCon(:,32),'r+-')
%semilogx(V,sigmaCon(:,6),'k+-',V,sigmaCon(:,12),'c+-',V,sigmaCon(:,14),'b+-',V,sigmaCon(:,16),'+-y',V,sigmaCon(:,18),'+-g',V,sigmaCon(:,20),'^-k')
ylabel('Ar')
xlabel('V')
legend(int2str(depth(12)),int2str(depth(18)),int2str(depth(22)),int2str(depth(26)),int2str(depth(30)),int2str(depth(32)))
hold on
figure(2) %plots shear strenth as function of depth
semilogx(V,friction(:,12),'*-r',V,friction(:,18),'*-g',V,friction(:,22),'*-c',V,friction(:,26),'*-b',V,friction(:,30),'o-y',V,friction(:,32),'o-k')
ylabel('friction')
xlabel('V')
legend(int2str(depth(12)),int2str(depth(18)),int2str(depth(22)),int2str(depth(26)),int2str(depth(30)),int2str(depth(32)))
hold on
figure(3) %plots shear strenth as function of depth
plot(tau(1,1:32)/10^6,depth(1:32),'r-',tau(2,1:34)/10^6,depth(1:34),'g-',tau(3,1:36)/10^6,depth(1:36),'k-')
xlabel('shear strength (MPa)')
ylabel('depth (Km)')
hold on
% figure(1) %plots shear strenth as function of depth
% plot(sigmaCon(1,1:31), depth(1:31),'r-',sigmaCon(2,1:33),depth(1:33),'g-',sigmaCon(3,1:35),depth(1:35),'k-')
% xlabel('contact stress (MPa)')
% ylabel('depth (Km)')
% hold on
%
%from joual et al
% A=12.4; % to add powerlaw creep also from 'as is ' in table 1 of Jaoul et al 1984.
% n=2.3; % those numbers are converted for correct units.
% E=171544;
% for i=1:3
% eps_dot=10^(-12+i);
% stress_power(i,:)=(eps_dot/A)^(1/n)*exp(E/R/n./T0);
% end
%from Hirth et al 2001
figure(3)
A=10^(-11.2); %
n=4; % those numbers are converted for correct units.
E=135000;
fH20=37*10^6;
for i=1:3
width=100*3.1^(i-1)
slip_rate=V(i)
eps_dot=slip_rate/width %this gives a 100m wide for 0.3, 310 meter for 3, 960 m for 30.
stress_power(i,:)=(eps_dot/A/fH20)^(1/n)*exp(E/R/n./T0);
end
plot(stress_power(1,32:60)*10^2,depth(32:60),'r-',stress_power(2,34:60)*10^2,depth(34:60),'g-',stress_power(3,36:60)*10^2,depth(36:60),'k-')