%-------------------------------------------------------
% Matlab code to do Backus averaging as in Kumar (2013),
% "Applying Backus averaging for deriving seismic anisotropy
% of long-wavelength equivalent medium from well-log data"
% Backus averaging of isotropic/VTI mixed layers
% This program is used to compute data and plot Figure 3
% Author: Dhananjay Kumar, Email: dhananjaykumar@gmail.com
% Input: layer properties in terms of Thomsen parameters (VTI)
% Output: Thomsen's parameters (VTI) of the Effective Medium
% How to run: save this file and type file.m on the matlab prompt
% Workflow:
% 1. Input parameters (need user input)
% 2. Computation as described in paper
% 3. Matlab plotting
%-------------------------------------------------------
clear all;
%-------------------------------------------------------
%1. Input parameters (Table 1)
%-------------------------------------------------------
no_rock = 2; %no of rock types present
% rock type 1 - shale (VTI)
vp1 = 3060; %m/s
vs1 = 1490; %m/s
rho1 = 2.42; %g/cc
eps1 = 0.256; %Thomsen's parameters
gam1 = 0.481;
del1 = -0.051;
% rock type 2 - sand (Isotropic)
vp2 = 2950; %m/s
vs2 = 1480; %m/s
rho2 = 2.0; %g/cc
eps2 = 0.0; %Thomsen's parameters
gam2 = 0.0;
del2 = 0.0;
%-------------------------------------------------------
% 2. Computation
%-------------------------------------------------------
%Layer's properties (Cij for VTI model)
c33_1 = vp1*vp1*rho1; % rock type 1
c44_1 = vs1*vs1*rho1;
c11_1 = (1+2*eps1)*c33_1;
c66_1 = (1+2*gam1)*c44_1;
c13_1 = sqrt(2*del1*c33_1*(c33_1-c44_1)+(c33_1-c44_1)^2)-c44_1;
c33_2 = vp2*vp2*rho2; % rock type 2
c44_2 = vs2*vs2*rho2;
c11_2 = (1+2*eps2)*c33_2;
c66_2 = (1+2*gam2)*c44_2;
c13_2 = sqrt(2*del2*c33_2*(c33_2-c44_2)+(c33_2-c44_2)^2)-c44_2;
% Build one set of Cij matrix for all rock types
for i = 1:no_rock
if rem(i,2) == 1 %layer1 (shale)
c11(i) = c11_1;
c33(i) = c33_1;
c44(i) = c44_1;
c13(i) = c13_1;
c66(i) = c66_1;
elseif rem(i,2) == 0 %layer2 (sand)
c11(i) = c11_2;
c33(i) = c33_2;
c44(i) = c44_2;
c13(i) = c13_2;
c66(i) = c66_2;
end
end
j = 0;
for rock1_fraction = 0:0.01:1
j = j+1;
rock1_vol_fraction(j) = rock1_fraction;
v(1) = rock1_fraction; %volume fraction: rock1
v(2) = 1-v(1); %volume fraction: rock2
c33_e = 1/sum((v./c33)); %Cij of effective model (eqns 8a-8e)
c44_e = 1/sum((v./c44));
c13_e = sum(v.*(c13./c33))*c33_e;
c66_e = sum(v.*c66);
c11_e = sum(v.*c11)+c13_e*sum(v.*(c13./c33))-sum(v.*(c13.*c13./c33));
epsilon(j) = 0.5*(c11_e-c33_e)/c33_e;
delta(j) = 0.5/(c33_e*(c33_e-c44_e))*((c13_e+c44_e)^2-(c33_e-c44_e)^2);
gamma(j) = 0.5*(c66_e - c44_e)/c44_e;
rho_e = v(1)*rho1 + v(2)*rho2;
% Note that vertical Vp and Vs are independent of intrinsic VTI
vp_e(j) = sqrt(c33_e/rho_e); %m/s
vs_e(j) = sqrt(c44_e/rho_e); %m/s
end % end of for loop
%-------------------------------------------------------
% 3. Plot output: Thomsen parameters as a function of rock1 volume
%-------------------------------------------------------
plot(rock1_vol_fraction,epsilon,'k.','linewidth',2);
hold on;
plot(rock1_vol_fraction,delta,'k-','linewidth',2);
hold on;
plot(rock1_vol_fraction,gamma,'k--','linewidth',2);
hold on;
ylim ([-0.1 0.5]);
xlabel ('Rock type 1 (volume fraction)');
ylabel ('VTI Anisotropy');
legend ('epsilon','delta','gamma');
%-------------------------END---------------------------
