clc
clear
MkSize=12;
LdWidth=2;
FtSize=20;
FtName='Arial';
FtSize2=16;
%%%%%% read the laboratory data of Fontainebleau sandstone
v = xlsread('C:\Users\DELL\Desktop\桌面文献及材料\Fontainebleau-Vosges sandstone data\Fontainebleau sandstone\Fontain_data.xlsx');
vpw_linear = xlsread('vp_linear.xls');  %% saturated vp computed by Liu's model (Liu et al., 2021)
vs_linear =  xlsread('vs_linear.xls');  %% saturated vs computed by Liu's model (Liu et al., 2021)
Kw_linear =  xlsread('Kw_linear.xls');  %% saturated bulk modulus computed by Liu's model (Liu et al., 2021)
muw_linear = xlsread('Muw_linear.xls'); %% saturated shear modulus computed by Liu's model (Liu et al., 2021)
Kw_DZ =  xlsread('K_Fon_DZ_model.xls'); %% saturated bulk modulus computed by D-Z model (David and Zimmerman, 2012)
muw_DZ = xlsread('Mu_Fon_DZ_model.xls');%% saturated shear modulus computed by D-Z model (David and Zimmerman, 2012)

p_lab = v(:,1)*1e6;     %% effective pressure
vpdry_lab = v(:,2)*1e3; %% laboratory measured dry vp
vpw_lab = v(:,4)*1e3;   %% laboratory measured saturated vp
vsdry_lab = v(:,5)*1e3; %% laboratory measured dry vs
vsw_lab = v(:,7)*1e3;   %% laboratory measured saturated vs
den_d     = 2544;       %% dry density
rho_total = 2570;       %% total(saturated) density

%% DZ-velocity
vsw_DZ = sqrt(muw_DZ(:,121)./rho_total);                       %% saturated vs from D-Z model
vpw_DZ = sqrt((Kw_DZ(:,121)+(4/3).*muw_DZ(:,121))./rho_total); %% saturated vp from D-Z model
%%%%%%%%% moduli_dry_g_w
Kdry_lab = den_d.*(vpdry_lab.^2-(4/3)*vsdry_lab.^2);%% laboratory measured dry bulk modulus
mudry_lab = den_d.*vsdry_lab.^2;                    %% laboratory measured dry shear modulus
Kw_lab = rho_total.*(vpw_lab.^2-(4/3)*vsw_lab.^2);  %% laboratory measured saturated bulk modulus
muw_lab = rho_total.*vsw_lab.^2;                    %% laboratory measured saturated shear modulus

%% Basic parameters for modelling
Kgr =    37e9;   %% solid bulk modulus
mu_gr =  44e9;   %% solid shear modulus
fai_s0 = 0.04;   %% Initial stiff porosity
etagly = 1e-3;   %% Fluid viscosity 
Kf =   2.24e9;   %% Fluid bulk modulus
rho_f =  1000;   %% Fluid density
k0 =     1*10e-15; %% Permeability
at = 2;            %% Tortuosity
K_drs =  32.8e9;   %% Dry bulk modulus without cracks
mu_drs = 39.2e9;   %% Dry shear modulus without cracks
lamda_drs = K_drs-2*mu_drs/3; %% Dry Lame parameter without cracks
pos_drs = (3*K_drs-2*mu_drs)./(6*K_drs+2*mu_drs);%% Dry Poisson's ratio without cracks
j = sqrt(-1);        %% Imaginary unit
fai_c0 = 1.69e-4;    %% Crack porosity
theta_c = 2755;      %% Pressure sensitivity coefficient
Cdrs = 1/K_drs;      %% Dry bulk compressibility with out cracks
Cgr = 1/Kgr;         %% Solid bulk compressibility
theta_cK = -20100e9; %% Pressure sensitivity coefficient
theta_cmu = -45500e9;%% Pressure sensitivity coefficient  
theta_cpasai1s = 0e9;%% Pressure sensitivity coefficient, no effect    
theta_cpsai2s =  0e9;%% Pressure sensitivity coefficient, no effect
psai_1s = 240e9;     %% Combined third-order elastic moduli, 3oEM
psai_2s = 3000e9;    %% 3oEM     
gamma = -450e9;      %% 3oEM    
gamma1 = 0e9;        %% 3oEM, no effect
gamma2 = -60e9;      %% 3oEM,
gamma3 =  72e9;      %% 3oEM,

f = 1e6;             %% Frequency
w = 2*pi*f;          %% Angular frequency
p_close = 100e6;      %% Close pressure
arf0 = (4-4*pos_drs^2)./(3*pi*K_drs*(1-2*pos_drs))*p_close; %% Initial aspect ratio

i = 1;
for pc = (10:5:95)*1e6  %% Confining pressure 
pf = 5e6;    %% Pore pressure              
p = pc-pf;   %% Effective pressure 
arf(i) = arf0 - (4-4*pos_drs^2)./(3*pi*K_drs*(1-2*pos_drs))*p; %% Pressure-dependent aspect ratio

rho_wave = at.*rho_f./(fai_s0-p.*(Cdrs-Cgr))-j*etagly./(k0.*w);

varib = 1/(K_drs+theta_cK.*fai_c0*exp(-theta_c.*p.*Cdrs))-1/K_drs;

ka = sqrt(3*j*w*etagly/Kf)/arf(i);
Kf_x = (1-2*besseli(1,ka)/(ka*besseli(0,ka)))*Kf;      %% Modified fluid modulus from Equ.(20) in article

Kmf(i) = 1./(1./K_drs+1./(1/varib...
    +1/((1/Kf_x-1/Kgr)*fai_c0*exp(-theta_c*p*Cdrs)))); %% Modified bulk modulus from Equ.(19-1) in article

mu_mf(i) = 1./(1/(mu_drs+theta_cmu.*fai_c0*exp(-theta_c.*p.*Cdrs))...
    -(4/15)*(1/(K_drs+theta_cK.*fai_c0*exp(-theta_c.*p.*Cdrs))-1/Kmf(i)));%% Modified shear modulus from Equ.(19-2) in article

lamda_mf = Kmf(i)-(2/3)*mu_mf(i); %% Modified Lame parameter from Equ.(23)
n = 1-Kmf(i)/Kgr; %% Biot coefficient from Equ.(28-1) in article

M  = 1./((n-fai_s0+p.*(Cdrs-Cgr))./Kmf(i)+(fai_s0-p.*(Cdrs-Cgr))./Kf);%% Equ.(28-2) in article
es = (pc-(fai_s0-(pc-pf).*(Cdrs-Cgr))*pf)./((3-3*fai_s0+3*(pc-pf)*(Cdrs-Cgr)).*Kgr); %% Equ.(28-3) in article
theta = (pc-(1-Kmf(i)./Kgr).*pf)./Kmf(i); %% Equ.(28-5) in article
ef = (theta-3*(1-fai_s0+(pc-pf).*(Cdrs-Cgr)).*es)./(3*fai_s0-3.*(pc-pf).*(Cdrs-Cgr));%% Equ.(28-4) in article

E1313 = mu_mf(i)-(n*M+gamma/2)*ef...
    +(3*lamda_mf+4*mu_mf(i)+psai_1s+theta_cpasai1s*fai_c0...
    *exp(-theta_c*p*Cdrs)+3*n^2*M)*es; %% Equ.(27-1) in article
E3333 = lamda_mf+2*mu_mf(i)+n.^2.*M...
    -(n.*M-2*gamma2).*ef+(5*lamda_mf+6*mu_mf(i)+psai_2s...
    +theta_cpsai2s.*fai_c0.*exp(-theta_c.*p.*Cdrs)+5*n^2.*M).*es; %% Equ.(27-2) in article
R33 = n.*M+(n.*M-2*gamma-6*gamma2).*es-(2*gamma3+n.*M).*ef; %% Equ.(27-3) in article
R = M+(6*gamma3+3*n.*M).*es-(3*M-6*gamma1).*ef; %% Equ.(27-4) in article

vp_fast(i) = sqrt((rho_total.*R+E3333.*rho_wave-2*rho_f.*R33)./(2*(rho_total.*rho_wave-rho_f.^2))...
    .*(1+sqrt(1-(4*(rho_total.*rho_wave-rho_f.^2).*(R.*E3333-R33.^2))./(rho_total.*R...
    +rho_wave.*E3333-2*rho_f.*R33).^2))); %% Fast P-wave velocity from Equ.(5) in article

vs(i) = sqrt(E1313./(rho_total-rho_f.^2/rho_wave)); %% S-wave velocity from Equ.(6) in article
i = i+1;

end

Kw_pre = rho_total.*(vp_fast.^2-(4/3)*vs.^2); %% Computed saturated bulk modulus 
muw_pre = rho_total.*vs.^2;  %% Computed saturated shear modulus 


p1 = 5:5:90;
pp=0:0.1:100;

figure(1)
set(gcf,'Position',[500 150 800 600])
plot(p1,vpw_lab/1e3,'ko','MarkerSize',10,'MarkerFaceColor','k');
hold on
plot(p1,vpdry_lab/1e3,'ko','MarkerSize',10);
hold on
plot(p1,vpw_linear/1e3,'k--','linewidth',LdWidth);
hold on
plot(pp,vpw_DZ/1e3,'m--','linewidth',LdWidth);
hold on
plot(p1,vp_fast/1e3,'k-','linewidth',LdWidth);
xlabel('Effective pressure (MPa)','FontName',FtName,'FontSize',FtSize);
ylabel('P-wave velocity (km/s)','FontName',FtName,'FontSize',FtSize);
axis([0 90 5.1 5.9])
set(gca,'xtick',0:10:90,'ytick',5.1:0.2:5.9,'FontName',FtName,'FontSize',FtSize)
lgd=legend(' Measured Vp, saturated',' Measured Vp, dry',' Predicted  Vp, saturated, Liu model',...
    ' Predicted  Vp, saturated, DZ model',' Predicted  Vp, saturated, new model','Location','Southeast');
legend('boxoff')
title(lgd,'Fontainebleau \phi_{\its} = 4%','FontName',FtName,'FontSize',FtSize,'FontWeight','Normal')

figure(2)
set(gcf,'Position',[500 150 800 600])
plot(p1,vsw_lab/1e3,'ko','MarkerSize',10,'MarkerFaceColor','k');
hold on
plot(p1,vsdry_lab/1e3,'ko','MarkerSize',10);
hold on
plot(p1,vs_linear/1e3,'k--','linewidth',LdWidth);
hold on
plot(pp,vsw_DZ/1e3,'m--','linewidth',LdWidth);
hold on
plot(p1,vs/1e3,'k-','linewidth',LdWidth);
xlabel('Effective pressure (MPa)','FontName',FtName,'FontSize',FtSize);
ylabel('S-wave velocity (km/s)','FontName',FtName,'FontSize',FtSize);
axis([0 90 3.5 4])
set(gca,'xtick',0:10:90,'ytick',3.5:0.1:4,'FontName',FtName,'FontSize',FtSize)
lgd=legend(' Measured Vs, saturated',' Measured Vs, dry',' Predicted  Vs, saturated, Liu model',...
    ' Predicted  Vs, saturated, DZ model',' Predicted  Vs, saturated, new model','Location','Southeast');
legend('boxoff')
title(lgd,'Fontainebleau \phi_{\its} = 4%','FontName',FtName,'FontSize',FtSize,'FontWeight','Normal')


figure(3)
set(gcf,'Position',[500 150 800 600])
plot(p1,Kw_lab/1e9,'ko','MarkerSize',10,'MarkerFaceColor','k');
hold on
plot(p1,Kdry_lab/1e9,'ko','MarkerSize',10);
hold on
plot(p1,Kw_linear/1e9,'k--','linewidth',LdWidth);
hold on
plot(pp,Kw_DZ(:,121)/1e9,'m--','linewidth',LdWidth);
hold on
plot(p1,Kw_pre/1e9,'k-','linewidth',LdWidth);
xlabel('Effective pressure (MPa)','FontName',FtName,'FontSize',FtSize);
ylabel('Bulk modulus (GPa)','FontName',FtName,'FontSize',FtSize);
axis([0 90 22 42])
set(gca,'xtick',0:10:90,'ytick',22:5:42,'FontName',FtName,'FontSize',FtSize)
lgd=legend(' Measured K, saturated',' Measured K, dry',' Predicted  K, saturated, Liu model',...
    ' Predicted  K, saturated, DZ model',' Predicted  K, saturated, new model','Location','Southeast');
legend('boxoff')
title(lgd,'Fontainebleau \phi_{\its} = 4%','FontName',FtName,'FontSize',FtSize,'FontWeight','Normal')

figure(4)
set(gcf,'Position',[500 150 800 600])
plot(p1,muw_lab/1e9,'ko','MarkerSize',10,'MarkerFaceColor','k');
hold on
plot(p1,mudry_lab/1e9,'ko','MarkerSize',10);
hold on
plot(p1,muw_linear/1e9,'k--','linewidth',LdWidth);
hold on
plot(pp,muw_DZ(:,121)/1e9,'m--','linewidth',LdWidth);
hold on
plot(p1,muw_pre/1e9,'k-','linewidth',LdWidth);
xlabel('Effective pressure (MPa)','FontName',FtName,'FontSize',FtSize);
ylabel('Shear modulus (GPa)','FontName',FtName,'FontSize',FtSize);
axis([0 90 22 42])
set(gca,'xtick',0:10:90,'ytick',22:5:42,'FontName',FtName,'FontSize',FtSize)
lgd=legend(' Measured \mu, saturated',' Measured \mu, dry',' Predicted  \mu, saturated, Liu model',...
    ' Predicted  \mu, saturated, DZ model',' Predicted  \mu, saturated, new model','Location','Southeast');
legend('boxoff')
title(lgd,'Fontainebleau \phi_{\its} = 4%','FontName',FtName,'FontSize',FtSize,'FontWeight','Normal')

p2 = 0:5:90;
arf2(1) = arf0;
arf2(2:19) = arf;
figure(5)
set(gcf,'Position',[500 150 800 600])
plot(p2,arf2,'k--','MarkerSize',10,'MarkerFaceColor','k');
xlabel('Effective pressure (MPa)','FontName',FtName,'FontSize',FtSize);
ylabel('Crack aspect ratio (-)','FontName',FtName,'FontSize',FtSize);
axis([0 90 1e-4 15e-4])
set(gca,'xtick',0:10:90,'ytick',1e-4:2e-4:15e-4,'FontName',FtName,'FontSize',FtSize)
lgd=legend(' ','Location','northeast');
legend('boxoff')
title(lgd,'Fontainebleau \phi_{\its} = 4%','FontName',FtName,'FontSize',FtSize,'FontWeight','Normal')
