diff --git a/calc_rheology.m b/calc_rheology.m
new file mode 100644
index 0000000..0237d7a
--- /dev/null
+++ b/calc_rheology.m
@@ -0,0 +1,60 @@
+function [differential_stress] = calc_rheology(mu,rho_uc,rho_lc,pore_fluid_factor,z,uc_thickness,crustal_thickness,crustal_strain_rate,mantle_strain_rate,n_uc,Q_uc,log_A_uc,n_lc,Q_lc,log_A_lc,Q_ml,n_ml,m_ml,d_ml,log_A_ml,T)
+% mu = coefficient of friction
+% rho_uc = density of material in upper crust
+% rho_lc = density of material in lower crust
+% pore fluid factor = 
+% z = depth in m
+% uc_thickness = upper crustal thickness (in m)
+% crustal_thickness = crustal thickness in m
+% crustal_strain_rate = in s^-1
+
+%% Constants
+g = 9.8; %accelleration due to gravity
+R = 8.3145; %universal gas constant
+
+%% SETUP PARAMETERS
+frictional_strength = NaN(size(z));
+viscous_strength = NaN(size(z));
+
+%% Upper Crust
+alpha = 2 * mu/((1 + mu^2)^0.5 + mu); %From Brace and Kohlstead, 1980 (JGR) and replicated elsewhere
+
+frictional_strength(z<=uc_thickness) = alpha * rho_uc * g .* z(z<=uc_thickness) * (1 - pore_fluid_factor);
+
+%% Lower Crust
+fs_update = max(frictional_strength);
+frictional_strength(z>uc_thickness) = fs_update + (alpha * rho_lc * g .* (z(z>uc_thickness) - uc_thickness) * (1 - pore_fluid_factor));
+
+%% Viscous Regime (Tempterature dependent)
+
+%crust - dislocation creep
+A_uc = (10^log_A_uc) * 1e6^-n_uc; %convert to Pa (*1e6) - original unit MPa^-n µm^m s^-1 but m=0
+A_lc = (10^log_A_lc) * 1e6^-n_lc;
+
+%mantle lithosphere - diffusion creep
+A_ml = (10^log_A_ml * 1e6^-n_ml) / 1e6^m_ml;
+
+%upper crust - dislocation creep
+viscous_strength(z<=uc_thickness) = (crustal_strain_rate./(A_uc * exp(-Q_uc./(R.*T(z<=uc_thickness))))).^(1/n_uc);
+
+%lower crust - dislocation creep
+viscous_strength(z>uc_thickness & z<=crustal_thickness) = (crustal_strain_rate./(A_lc * exp(-Q_lc./(R.*T(z>uc_thickness & z<=crustal_thickness))))).^(1/n_lc);
+
+%mantle lithosphere - diffusion creep
+viscous_strength(z>crustal_thickness) = (mantle_strain_rate./(A_ml*(d_ml^-m_ml) * exp(-Q_ml./(R.*T(z>crustal_thickness))))).^(1/n_ml);
+
+%% Differential Stress
+%calculate minimum by combining frictional and viscous curves and convert
+%to MPa
+differential_stress=min(frictional_strength,viscous_strength)./1e6;
+
+% plot(frictional_strength./1e6,z./1e3)
+% hold on
+% plot(viscous_strength./1e6,z./1e3)
+% set(gca,'YDir','reverse')
+% xlabel('\sigma_{1} - \sigma_{3} (MPa)')
+% ylabel('Depth (km)')
+% xlim([0 550])
+% ylim([0 50])
+
+
diff --git a/calc_temperature.m b/calc_temperature.m
new file mode 100644
index 0000000..e54a4f7
--- /dev/null
+++ b/calc_temperature.m
@@ -0,0 +1,25 @@
+function [T] = calc_temperature(z,Q0,H_uc,K_uc,T0,H_lc,K_lc,uc_thickness,lc_thickness)
+% z = depth in m
+% Q0 = surface heat flow value (W m^-2)
+% H_uc = surface heat production in the upper crust (Wm^-3)
+% K_uc = thermal conductivity in the upper crust (Wm^-1 K^-1)
+% T0 = initial temperature
+% H_lc = surface heat production in the lower crust (Wm^-3)
+% K_lc = thermal conductivity in the lower crust (Wm^-1 K^-1)
+% uc_thickness = thickness of the upper crust (in m)
+% lc_thickness = thickness of the lower crust (in m)
+
+T = zeros(size(z));
+
+%% Calculate Upper crust temperature
+
+T(z<=uc_thickness) = T0 + (Q0.*z(z<=uc_thickness)/K_uc) - (H_uc.*z(z<=uc_thickness).^2)/(2*K_uc);
+
+%update values for the lower crust
+z_up = z(z>uc_thickness) - max(z(z<=uc_thickness));
+T0 = max(T);
+Q0 = Q0 - H_uc*uc_thickness;
+
+T(z>uc_thickness) = T0 + (Q0.*z_up/K_lc) - (H_lc.*z_up.^2)/(2*K_lc);
+
+% 
\ No newline at end of file
diff --git a/post_seismic_strain_rates.m b/post_seismic_strain_rates.m
new file mode 100644
index 0000000..cb0916f
--- /dev/null
+++ b/post_seismic_strain_rates.m
@@ -0,0 +1,284 @@
+% Postseismic deformation in a power law shear zone from Montesi (2004)
+clear all
+close all
+
+% Vs(t) = V0[1+(1-1/n)t/tau]^(-1/(1-1/n))
+%let g=(1-1/n)
+% Vs(t) = V0[1+g*t/tau]^(-1/g)
+% differentiate to fine the strain rate
+% epsilon_dot=V0[-1/g - t/tau]^(-1/g)-1
+
+%scenario - initial velocity of 1 mm/yr with a maxwell decay time of 50
+%years 
+
+V0=0.0005; %1 mm/yr
+V0=0.5; 
+n=3;
+tau=0.09; %in years
+t=0:0.1:10000;
+
+vs = V0*(1+(1-1/n).*t./tau).^(-1/(1-1/n));
+n=5;
+vsn5=V0*(1+(1-1/n).*t./tau).^(-1/(1-1/n));
+n=1.000000001;
+vsn1=V0*(1+(1-1/n).*t./tau).^(-1/(1-1/n));
+
+%strain rate - divide by a distance varying from 100 m to 1m
+epsilon_dot_100=vs/100;
+epsilon_dot_10 = vs/10;
+epsilon_dot_1 = vs/1;
+
+epsilon_dot_n5_100=vsn5/100;
+epsilon_dot_n5_10 = vsn5/10;
+epsilon_dot_n5_1 = vsn5/1;
+
+epsilon_dot_n1_100=vsn1/100;
+epsilon_dot_n1_10 = vsn1/10;
+epsilon_dot_n1_1 = vsn1/1;
+
+%%
+z = [0:100:100000];
+%Calculate temperature
+%Q0 = 63e-3; %mW m^-2 from Nyblade 1990
+Q0 = 68e-3; %the mean value from the southwestern and western rift branches
+H=2.1e-6;
+K=2.8;
+T0=273.15; %zero degrees celsius in kelvin
+
+T_rift=zeros(1,401);
+T_rift(z<=9000) = T0 + ((Q0.*z(z<=9000))/K) - ((H.*z(z<=9000).^2)/(2*K));
+
+%update values for lower crust
+z_up = z(z>9000)-max(z(z<=9000));
+T0=max(T_rift);
+Q0=Q0-H*9000;
+H=0.4e-6;
+K=2.0;
+
+T_rift(z>9000) = T0 + ((Q0.*z_up)./K) - ((H.*(z_up.^2))/(2*K));
+
+%%Calculate differential stress
+
+%Frictional regime
+mu = 0.60;
+%alpha = (sqrt(1+mu^2))^-2; %original equation given to me by Ake who
+%claimed it's adapted from Sibson 1974 (Nature) - but not seen replicated
+%elsewhere and weirdly results in an increase in fault strength with a
+%reduction in friction
+alpha = 2*mu/((1+mu^2)^0.5 + mu); %From Brace and Kohlstead, 1980 (JGR) and replicated elsewhere
+rho_uc = 2650;
+rho_lc = 2800;
+g = 9.8;
+lambda_v = 0.4;
+
+frictional_strength = nan(size(z));
+
+%upper crust
+frictional_strength(z<=16000) = alpha * rho_uc * g .* z(z<=16000) * (1-lambda_v);
+
+%lower crust
+fs_update = max(frictional_strength);
+frictional_strength(z>16000) = fs_update + (alpha * rho_lc * g .* (z(z>16000)-16000) * (1-lambda_v));
+
+%viscous regime
+
+%--------------------------------
+%Dislocation Creep
+%--------------------------------
+strain_rate_background=1e-16;
+strain_rate_1_week = epsilon_dot_10(t==0.1)/31557600;
+strain_rate_1_yr = epsilon_dot_10(t==1)/31557600;
+strain_rate_10_yrs = epsilon_dot_10(t==10)/31557600;
+strain_rate_100_yrs = epsilon_dot_10(t==100)/31557600;
+strain_rate_1000_yrs = epsilon_dot_10(t==1000)/31557600;
+
+epsilon_dot_mantle_lithosphere = 1e-15;
+
+%universal gas constant
+R =  8.3145;
+
+crustal_thickness=41000; %median crustal thickness from receiver function measurements of 3 transects
+
+
+%-------------------
+%Quartz upper crust
+%from Rutter and Brodie 2004
+n_uc = 3;
+Q_uc = 242e3;
+A_uc = 10^-4.9 * 1e6^-n_uc;
+     
+%-------------------
+%FELSIC LOWER CRUST
+%Values for synthetic anorthite (feldspar) from Rybacki and Dresen, 2000;
+%JGR, Table 2.
+%dry
+n_flc = 3.0; %±0.4
+Q_flc = 648e3; %± 20; kJ mol^-1; convert to J.
+A_flc = 10^12.7 * 1e6^-n_flc; %LogA = 12.7 ± 0.8 MPa^-n µm^m s^-1; convert to Pa s-1 (div by 1e6^-n, m=0).
+
+%wet
+n_flc_w = 3.0; %±0.2
+Q_flc_w = 356e3; %± 9; kJ mol^-1; convert to J.
+A_flc_w = 10^2.6 * 1e6^-n_flc_w; %LogA = 2.6 ± 0.3 MPa^-n µm^m s^-1; convert to m (div by 1e6).
+%-------------------
+%DRY OLIVINE MANTLE LITHOSPHERE
+%from Hirth and Kohlstedt, 2003 - taken from Burgmann and Dresen 2009
+n_doml_disc = 3.5;
+Q_doml_disc = 530e3; %±4; kJ mol^-1; convert to J.
+A_doml_disc = 10^5 * 1e6^-n_doml_disc; %LogA in MPa^-n µm^m s^-1; convert to m (div by 1e6).
+
+%--------------------------------
+%Diffusion Creep
+%--------------------------------
+%for use in the mantle lithosphere below the rifts.
+
+%-------------------
+%DRY OLIVINE MANTLE LITHOSPHERE
+%from Faul and Jackson, 2007, JGR. see paragraph 32.
+Q_doml_difc = 484e3; %±30
+n_doml_difc = 1.37; %±0.06
+m_doml_difc = 3;
+d_doml_difc = 50/1e6; %grain size range given as 2.7-5.6µm
+A_doml_difc = (10^10.3 * 1e6^-n_doml_difc) / 1e6^m_doml_difc; % ± 4 MPa^-n µm^m s^-1; convert to m (div by 1e6).
+
+%--------------------------------
+%calculate viscous strength curves
+%--------------------------------
+
+viscous_stress_background = NaN(size(z));
+viscous_stress_1_week = NaN(size(z));
+viscous_stress_1_yr = NaN(size(z));
+viscous_stress_10_yrs = NaN(size(z));
+viscous_stress_100_yrs = NaN(size(z));
+viscous_stress_1000_yrs = NaN(size(z));
+
+%upper crust
+%dislocation creep of Quartz
+viscous_stress_background(z<=9000) = (strain_rate_background./(A_uc*exp(-Q_uc./(R.*T_rift(z<=9000))))).^(1/n_uc);
+viscous_stress_1_week(z<=9000) = (strain_rate_1_week./(A_uc*exp(-Q_uc./(R.*T_rift(z<=9000))))).^(1/n_uc);
+viscous_stress_1_yr(z<=9000) = (strain_rate_1_yr./(A_uc*exp(-Q_uc./(R.*T_rift(z<=9000))))).^(1/n_uc);
+viscous_stress_10_yrs(z<=9000) = (strain_rate_10_yrs./(A_uc*exp(-Q_uc./(R.*T_rift(z<=9000))))).^(1/n_uc);
+viscous_stress_100_yrs(z<=9000) = (strain_rate_100_yrs./(A_uc*exp(-Q_uc./(R.*T_rift(z<=9000))))).^(1/n_uc);
+viscous_stress_1000_yrs(z<=9000) = (strain_rate_1000_yrs./(A_uc*exp(-Q_uc./(R.*T_rift(z<=9000))))).^(1/n_uc);
+
+%mid- to lower-curst
+%dislocation creep of feldspar (anorthite)
+viscous_stress_background(z>9000 & z<=crustal_thickness) = (strain_rate_background./(A_flc*exp(-Q_flc./(R.*T_rift(z>9000 & z<=crustal_thickness))))).^(1/n_flc);
+viscous_stress_1_week(z>9000 & z<=crustal_thickness) = (strain_rate_1_week./(A_flc*exp(-Q_flc./(R.*T_rift(z>9000 & z<=crustal_thickness))))).^(1/n_flc);
+viscous_stress_1_yr(z>9000 & z<=crustal_thickness) = (strain_rate_1_yr./(A_flc*exp(-Q_flc./(R.*T_rift(z>9000 & z<=crustal_thickness))))).^(1/n_flc);
+viscous_stress_10_yrs(z>9000 & z<=crustal_thickness) = (strain_rate_10_yrs./(A_flc*exp(-Q_flc./(R.*T_rift(z>9000 & z<=crustal_thickness))))).^(1/n_flc);
+viscous_stress_100_yrs(z>9000 & z<=crustal_thickness) = (strain_rate_100_yrs./(A_flc*exp(-Q_flc./(R.*T_rift(z>9000 & z<=crustal_thickness))))).^(1/n_flc);
+viscous_stress_1000_yrs(z>9000 & z<=crustal_thickness) = (strain_rate_1000_yrs./(A_flc*exp(-Q_flc./(R.*T_rift(z>9000 & z<=crustal_thickness))))).^(1/n_flc);
+
+%mantle lithosphere
+%diffusion creep of dry olivine (rift)
+viscous_stress_background(z>crustal_thickness) = (epsilon_dot_mantle_lithosphere./(A_doml_difc*(d_doml_difc^-m_doml_difc)*exp(-Q_doml_difc./(R.*T_rift(z>crustal_thickness))))).^(1/n_doml_difc);
+viscous_stress_1_week(z>crustal_thickness) = (epsilon_dot_mantle_lithosphere./(A_doml_difc*(d_doml_difc^-m_doml_difc)*exp(-Q_doml_difc./(R.*T_rift(z>crustal_thickness))))).^(1/n_doml_difc);
+viscous_stress_1_yr(z>crustal_thickness) = (epsilon_dot_mantle_lithosphere./(A_doml_difc*(d_doml_difc^-m_doml_difc)*exp(-Q_doml_difc./(R.*T_rift(z>crustal_thickness))))).^(1/n_doml_difc);
+viscous_stress_10_yrs(z>crustal_thickness) = (epsilon_dot_mantle_lithosphere./(A_doml_difc*(d_doml_difc^-m_doml_difc)*exp(-Q_doml_difc./(R.*T_rift(z>crustal_thickness))))).^(1/n_doml_difc);
+viscous_stress_100_yrs(z>crustal_thickness) = (epsilon_dot_mantle_lithosphere./(A_doml_difc*(d_doml_difc^-m_doml_difc)*exp(-Q_doml_difc./(R.*T_rift(z>crustal_thickness))))).^(1/n_doml_difc);
+viscous_stress_1000_yrs(z>crustal_thickness) = (epsilon_dot_mantle_lithosphere./(A_doml_difc*(d_doml_difc^-m_doml_difc)*exp(-Q_doml_difc./(R.*T_rift(z>crustal_thickness))))).^(1/n_doml_difc);
+
+%combine different curves and convert to MPa
+stress_background = min(viscous_stress_background./1e6,frictional_strength./1e6);
+stress_1_week = min(viscous_stress_1_week./1e6,frictional_strength./1e6);
+stress_1_yr = min(viscous_stress_1_yr./1e6,frictional_strength./1e6);
+stress_10_yrs = min(viscous_stress_10_yrs./1e6,frictional_strength./1e6);
+stress_100_yrs = min(viscous_stress_100_yrs./1e6,frictional_strength./1e6);
+stress_1000_yrs = min(viscous_stress_1000_yrs./1e6,frictional_strength./1e6);
+
+%--------------------------------------------------------------------------
+%%Plot
+figure('Position',[200,200,700,500])
+subplot(2,2,1)
+loglog(t,vsn1.*1e3,'black')
+hold on
+loglog(t,vs.*1e3,'blue')
+loglog(t,vsn5.*1e3,'red')
+
+% plot data from Mozambique earthquake (data in Ingleby and Wright, 2017).
+time = [0.25,1.85,3.74];
+velocity = [104.29,17.36,2.83]; %in mm/yr
+loglog(time,velocity,'+')
+%loglog(t,vt100.*1e3)
+
+txt = '(a)';
+text(0,0,txt,'Position',[0.05 0.95 0],'Units','normalized','FontSize',12,'FontName','helvetica')
+
+
+xlabel('time (years)')
+ylabel('velocity (mm/yr)')
+title('Postseismic deformation in shear zones')
+txt = {'V_0 = 500 mm/yr','\tau = 0.09 yrs^{-1}'};
+text(300,0.08,txt)
+txt = {'$$ V_s(t) = V_{0} \left[ 1 + \left( 1-\frac{1}{n} \right)\frac{1}{\tau} \right]^{ \left( \frac{-1}{1-1/n} \right)} $$'};
+text(2,100,txt,'Interpreter','latex')
+
+legend('n=1','n=3','n=5','Mozambique','Location','east')
+ylim([0.01 600])
+xlim([0.1 10000])
+
+subplot(2,2,3)
+loglog(t,epsilon_dot_100./31557600,'-.blue'); %divide to get strain rate per second
+hold on
+loglog(t,epsilon_dot_10./31557600,'-blue');
+loglog(t,epsilon_dot_1./31557600,'--blue');
+
+loglog(t,epsilon_dot_n5_100./31557600,'-.red');
+loglog(t,epsilon_dot_n5_10./31557600,'-red');
+loglog(t,epsilon_dot_n5_1./31557600,'--red');
+
+loglog(t,epsilon_dot_n1_100./31557600,'-.black');
+loglog(t,epsilon_dot_n1_10./31557600,'-black');
+loglog(t,epsilon_dot_n1_1./31557600,'--black');
+
+xlabel('time (years)')
+ylabel('strain rate (s^{-1})')
+legend('100 m wide shear zone','10 m wide shear zone','1 m wide shear zone','Location','northeast')
+title('Postseismic strain rate with time')
+%text(0.1,3e-11,'1 m wide shear zone')
+%text(0.1,3e-12,'10 m wide shear zone')
+%text(0.1,3e-13,'100 m wide shear zone')
+txt = '(b)';
+text(0,0,txt,'Position',[0.05 0.95 0],'Units','normalized','FontSize',12,'FontName','helvetica')
+
+
+xlim([0.1 10000])
+ylim([1e-16 1e-7])
+
+subplot(2,2,[2,4])
+plot(stress_background,z./1e3,'k','LineWidth',2)
+hold on
+plot(stress_1_yr,z./1e3,'LineWidth',1)
+plot(stress_10_yrs,z./1e3,'LineWidth',1)
+plot(stress_100_yrs,z./1e3,'LineWidth',1)
+plot(stress_1000_yrs,z./1e3,'LineWidth',1)
+set(gca,'YDir','reverse')
+xlabel('\sigma_{1} - \sigma_{3} (MPa)')
+ylabel('Depth (km)')
+xlim([0 500])
+ylim([0 50])
+title('Crustal strength changes during postseismic deformation')
+legend('background strain rate','1 yr','10 yrs','100 yrs','1000 yrs')
+
+txt = '(c)';
+text(0,0,txt,'Position',[0.05 0.98 0],'Units','normalized','FontSize',12,'FontName','helvetica')
+
+time_ps=[1,10,100,1000];
+
+%calculate BDT
+[~,y1yr] = max(stress_1_yr(z<=crustal_thickness));
+[~,y10yr] = max(stress_10_yrs(z<=crustal_thickness));
+[~,y100yr] = max(stress_100_yrs(z<=crustal_thickness));
+[~,y1000yr] = max(stress_1000_yrs(z<=crustal_thickness));
+bdt(1) = z(y1yr);
+bdt(2) = z(y10yr);
+bdt(3) = z(y100yr);
+bdt(4) = z(y1000yr);
+
+% figure
+% semilogx(time_ps,bdt./1e3,'+-')
+% xlabel('time (years)')
+% ylabel('brittle ductile transition (km)')
+% 
+% ylim([30 40])
diff --git a/run_rheology_calculator.m b/run_rheology_calculator.m
new file mode 100644
index 0000000..d3c4584
--- /dev/null
+++ b/run_rheology_calculator.m
@@ -0,0 +1,200 @@
+% run scripts
+close all
+clear all
+
+z = [0:100:100000];
+
+%temperature parameters
+Q0=68e-3;
+T0 = 273.15;
+
+H_uc = 2.1e-6;
+K_uc = 2.8;
+
+H_lc = 0.4e-6;
+K_lc = 2.0;
+
+uc_thickness = 9000;
+crustal_thickness = 41000;
+
+%rheological parameters
+mu = 0.4;
+rho_uc = 2650;
+rho_lc = 2800;
+pore_fluid_factor = 0.4;
+
+%strain rate parameters
+crustal_strain_rate = 1e-15;
+mantle_strain_rate = 1e-15;
+
+%dislocation creep parameters
+%for Quartz Upper crust from Rutter and Brodio 2004
+n_uc = 3;
+Q_uc = 242e3;
+log_A_uc = -4.9;
+
+%for dry Felsic lower crust
+%Values for synthetic anorthite (feldspar) from Rybacki and Dresen, 2000;
+%JGR, Table 2.
+n_lc = 3.0;
+Q_lc = 648e3;
+log_A_lc = 12.7;
+
+%Diffusion Creep Parameters
+%DRY OLIVINE MANTLE LITHOSPHERE
+%from Faul and Jackson, 2007, JGR. see paragraph 32.
+Q_ml = 484e3;
+n_ml = 1.37;
+m_ml = 3;
+d_ml = 50/1e6; %grain size range given as 2.7-5.6µm
+log_A_ml = 10.3;
+
+
+
+
+T = calc_temperature(z,Q0,H_uc,K_uc,T0,H_lc,K_lc,uc_thickness,crustal_thickness);
+diff_stress = calc_rheology(mu,rho_uc,rho_lc,pore_fluid_factor,z,uc_thickness,crustal_thickness,crustal_strain_rate,mantle_strain_rate,n_uc,Q_uc,log_A_uc,n_lc,Q_lc,log_A_lc,Q_ml,n_ml,m_ml,d_ml,log_A_ml,T);
+
+figure
+subplot(1,2,1)
+plot(T-273.15,z./1e3)
+hold on
+set(gca,'YDir','reverse')
+ylim([0 50])
+
+subplot(1,2,2)
+plot(diff_stress,z./1e3)
+hold on
+%plot(fric_strength./1e6,z./1e3)
+set(gca,'YDir','reverse')
+ylim([0 50])
+
+%vary strain rate and pore fliud factor and analyse depth of Brittle Ductil
+%Transition
+
+figure('Position',[200,200,600,1200])
+subplot(3,2,1)
+plot(T-273.15,z./1e3,'LineWidth',2) 
+hold on
+plot([0 1200],[9 9],'k--')
+plot([0 1200],[41 41],'k--')
+hold on
+set(gca,'YDir','reverse')
+ylabel('Depth (km)')
+xlabel('Temperature ({\circ}C)')
+ylim([0 50])
+title('Temperature profile')
+
+txt='upper crust: 0-9 km';
+text(225,7,txt)
+txt='mid-lower crust: 9-41 km';
+text(100,39,txt)
+txt='lithospheric mantle: 41+ km';
+text(100,48,txt)
+
+text(0,0,'(a)','Position',[0.93 0.96 0],'Units','normalized','FontSize',12,'FontName','helvetica')
+
+
+subplot(3,2,2)
+clear T
+T_2 = calc_temperature(z,Q0,H_uc,K_uc,T0,H_lc,K_lc,uc_thickness,25000);
+diff_stress_25 = calc_rheology(mu,rho_uc,rho_lc,pore_fluid_factor,z,uc_thickness,25000,crustal_strain_rate,mantle_strain_rate,n_uc,Q_uc,log_A_uc,n_lc,Q_lc,log_A_lc,Q_ml,n_ml,m_ml,d_ml,log_A_ml,T_2);
+plot(diff_stress_25,z./1e3,'--','LineWidth',2)
+hold on
+T = calc_temperature(z,Q0,H_uc,K_uc,T0,H_lc,K_lc,uc_thickness,crustal_thickness);
+diff_stress_41 = calc_rheology(mu,rho_uc,rho_lc,pore_fluid_factor,z,uc_thickness,crustal_thickness,crustal_strain_rate,mantle_strain_rate,n_uc,Q_uc,log_A_uc,n_lc,Q_lc,log_A_lc,Q_ml,n_ml,m_ml,d_ml,log_A_ml,T);
+plot(diff_stress_41,z./1e3,'r','LineWidth',2)
+set(gca,'YDir','reverse')
+xlabel('\sigma_{1} - \sigma_{3} (MPa)')
+ylabel('Depth (km)')
+xlim([0 550])
+ylim([0 50])
+title('Crustal thickness and lithospheric strength')
+l1 = sprintf('25 km: strength = %3.1f GPa', (trapz(diff_stress_25)/1e3));
+l2 = sprintf('41 km: strength = %3.1f GPa', (trapz(diff_stress_41)/1e3));
+
+legend(l1,l2);
+legend('boxoff')
+text(0,0,'(c)','Position',[0.93 0.96 0],'Units','normalized','FontSize',12,'FontName','helvetica')
+
+
+
+subplot(3,2,3)
+T = calc_temperature(z,Q0,H_uc,K_uc,T0,H_lc,K_lc,uc_thickness,crustal_thickness);
+hold on
+set(gca,'YDir','reverse')
+xlabel('\sigma_{1} - \sigma_{3} (MPa)')
+ylabel('Depth (km)')
+xlim([0 550])
+ylim([0 50])
+title('Change in crustal strain rate')
+text(0,0,'(c)','Position',[0.93 0.96 0],'Units','normalized','FontSize',12,'FontName','helvetica')
+
+
+%vary strain rate
+count = 1;
+for i=[logspace(-15,-10,6)]
+    crustal_strain_rate = i;
+    diff_stress = calc_rheology(mu,rho_uc,rho_lc,pore_fluid_factor,z,uc_thickness,crustal_thickness,crustal_strain_rate,mantle_strain_rate,n_uc,Q_uc,log_A_uc,n_lc,Q_lc,log_A_lc,Q_ml,n_ml,m_ml,d_ml,log_A_ml,T);
+    plot(diff_stress,z./1e3);
+    
+    %calculate the BDT
+    [~,y] = max(diff_stress(z<=crustal_thickness));
+    bdt_strain_rate(count) = z(y);
+    strain_rate(count) = i;
+    
+    count = count+1;
+end
+
+lgd = legend('1e-15 s^{-1}','1e-14 s^{-1}','1e-13 s^{-1}','1e-12 s^{-1}','1e-11 s^{-1}','1e-10 s^{-1}','Location','northeast');
+title(lgd,'strain rate')
+legend('boxoff')
+
+subplot(3,2,4)
+hold on
+set(gca,'YDir','reverse')
+xlabel('\sigma_{1} - \sigma_{3} (MPa)')
+ylabel('Depth (km)')
+xlim([0 550])
+ylim([0 50])
+title('Change in pore fluild factor')
+text(0,0,'(d)','Position',[0.93 0.96 0],'Units','normalized','FontSize',12,'FontName','helvetica')
+
+
+crustal_strain_rate = 1e-15;
+count = 1;
+for i = [0.4:0.1:0.9]
+    pore_fluid_factor = i;
+    diff_stress = calc_rheology(mu,rho_uc,rho_lc,pore_fluid_factor,z,uc_thickness,crustal_thickness,crustal_strain_rate,mantle_strain_rate,n_uc,Q_uc,log_A_uc,n_lc,Q_lc,log_A_lc,Q_ml,n_ml,m_ml,d_ml,log_A_ml,T);
+    plot(diff_stress,z./1e3);
+    
+    %calculate the BDT
+    [~,y] = max(diff_stress(z<=crustal_thickness));
+    bdt_pore_fluid(count) = z(y);
+    pore_fluid(count) = i;
+    
+    count = count+1;    
+end
+
+lgd=legend('0.4','0.5','0.6','0.7','0.8','0.9','Location','northeast');
+title(lgd,'pore fluid factor')
+legend('boxoff')
+
+subplot(3,2,5)
+semilogx(strain_rate,bdt_strain_rate./1e3,'-+')
+ylim([30 40])
+set(gca,'YDir','reverse')
+xlabel('crustal strain rate (s^{-1})')
+ylabel('brittle ductile transition (km)')
+text(0,0,'(e)','Position',[0.93 0.96 0],'Units','normalized','FontSize',12,'FontName','helvetica')
+
+
+subplot(3,2,6)
+plot(pore_fluid,bdt_pore_fluid./1e3,'-+')
+ylim([30 40])
+set(gca,'YDir','reverse')
+xlabel('pore fluild factor')
+ylabel('brittle ductile transition (km)')
+text(0,0,'(f)','Position',[0.93 0.96 0],'Units','normalized','FontSize',12,'FontName','helvetica')
+
+exportgraphics(gcf,'Figure_4.png','Resolution',300)