Viscosity profile of dry and wet matle(based on strain rate invariant)(compose dis and dif to eff).

% Figuring and comparing viscosity map for dry and wet mantle
% Based on No.2 stress invariant.
clear all;clf;
%set parameters of dislocation
A_dis = [3.5e22,2e18]; %unit :( 1/(Pa.^n.*S) )
n_dis = [3.5,3];
m_dis = [0,0];
Ea_dis = [540000,430000];%unit:(J/mol)

%set parameters of diffusion
A_dif = [8.7e15,5.3e15]; %unit :( 1/(Pa.^n.*S) )
n_dif = [1,1];
m_dif = [-2.5,-2.5];
Ea_dif = [300000,240000];%unit:(J/mol)

h = 0.001;%unit :(m)
R = 8.314;%unit :( J/(K*mol) ) 
miu = 8*10.^10; %unit:(Pa)
b = 5.*10.^(-10);%unit:(m)
top_tempeture = 400;%unit:(℃)
bottom_tempeture = 1400;%unit:(℃)
top_stress = 3;
bottom_stress = 9;

% set stress condition unit:(Pa)
stop = top_stress;
sbottom = bottom_stress;
snumber = 51;
S = stop:(sbottom-stop)/(snumber-1):sbottom;

% set T condition unit(K)
ttop = top_tempeture;
tbottom = bottom_tempeture;
tnumber = 51;
T = ttop:(tbottom-ttop)/(tnumber-1):tbottom;

% Setting V_dry and V_wet
V_dry = zeros(snumber,tnumber);
V_wet = zeros(snumber,tnumber);
diff = zeros(snumber,tnumber);
% Computing viscosity by Eq.6.4 and Eq.6.16
for j = 1:tnumber%j = 1:Nx
    for i = 1:snumber
        % Dry
        % Why??
        eterm = Ea_dis(1)/R/(T(j)+273);
%         if (eterm>100) 
%             eterm=100; 
%         end
        eterm=exp(-eterm);
        V_dis = (10.^S(i))/2/( A_dis(1)*((h/b).^m_dis(1))*((10^S(i)/miu).^n_dis(1))*eterm ) ;
        eterm = Ea_dif(1)/R/(T(j)+273);
%         if (eterm>100) 
%             eterm=100; 
%         end
        eterm=exp(-eterm);
        V_dif = (10.^S(i))/2/( A_dif(1)*((h/b).^m_dif(1))*((10^S(i)/miu).^n_dif(1))*eterm ) ;
        V_dry(i,j) = 1/(1/V_dis+1/V_dif);
        % Wet
        % Why??
        eterm = Ea_dis(2)/R/(T(j)+273);
%         if (eterm>100) 
%             eterm=100; 
%         end
        eterm=exp(-eterm);
        V_dis = (10.^S(i))/2/( A_dis(2)*((h/b).^m_dis(2))*((10^S(i)/miu).^n_dis(2))*eterm ) ;
        eterm = Ea_dif(2)/R/(T(j)+273);
%         if (eterm>100) 
%             eterm=100; 
%         end
        eterm=exp(-eterm);
        V_dif = (10.^S(i))/2/( A_dif(2)*((h/b).^m_dif(2))*((10^S(i)/miu).^n_dif(2))*eterm ) ;
        V_wet(i,j) = 1/(1/V_dis+1/V_dif);
        % Compare log(V) map for dry and wet condition.
        diff(i,j) =  log10(V_dry(i,j))-log10(V_wet(i,j));
    end
end

% Log10(viscosity in dry)
figure(1)
subplot(1,3,1);
pcolor(S,T,log10(V_dry));
xlabel('S');
ylabel('T');
title('Log10(viscosity in dry)');
set(gca,'xaxislocation','top');
set (gca,'YDir','reverse')
shading flat
colorbar

% Log10(viscosity in wet)
subplot(1,3,2);
pcolor(S,T,log10(V_wet));
xlabel('S');
ylabel('T');
title('Log10(viscosity in wet)');
set(gca,'xaxislocation','top');
set (gca,'YDir','reverse')
shading flat
colorbar

% Subtract Log10(wet) from Log10(dry)
subplot(1,3,3);
pcolor(S,T,diff);
xlabel('S');
ylabel('T');
title('Different between dry and wet.');
set(gca,'xaxislocation','top');
set (gca,'YDir','reverse')
shading flat
colorbar

The third picture​ look like the subduction.​​​​​​