Seismic Reflection/Transmission Coefficients with Matlab

A common computation in seismology is the calculation of reflection and transmission coefficients that describe the partitioning of energy when a seismic wave strikes a boundary between elastic materials. The coefficients depend on the seismic wave speeds and densities on either side of the boundary and on the incident wave horizontal slowness, or ray parameter. The formulas for the values can be found in Table 3.1 of Lay and Wallace (1995) (note that an error in the sign of the second term in the computation of b) or in equations 5.38 and 5.39 of Aki and Richards (1980). For a welded elastic boundary, an incident P or SV wave results in four waves, two reflected and two transmitted or refracted waves. Thus for each case we need four coefficients – notation for each is illustrated in the figure below.

Notation used to identify reflection and transmission coefficients. We must consider two cases, that for an incident P wave and that for an incident SV wave. R identifies “Reflected” waves, T represents “Transmitted” or “refracted” waves.

The Code

A program to compute the coefficient equations is straight forward so I won’t go into details (but always watch for typos). The implementation in the script listed below includes one twist. The arguments to the script include the velocity and density values for each material and the ray parameter – which can be a vector. You can call the function with a single value of the ray parameter, or a range of ray parameters stored in a vector. If you pass the function a vector of ray parameters, it will return a complex matrix containing the R/T coefficients.

function [RTmatrix] = PSVRTmatrix(p,mi,mt)
%
%    p = ray parameter (a scalar or row vector)
%
%    mi = model of incident    wave [Vp Vs Rho]
%    mt = model of transmitted wave [Vp Vs Rho]
%

% vertical slownesses
etaai = sqrt(1/(mi(1)*mi(1)) - p.*p);
etaat = sqrt(1/(mt(1)*mt(1)) - p.*p);
etabi = sqrt(1/(mi(2)*mi(2)) - p.*p);
etabt = sqrt(1/(mt(2)*mt(2)) - p.*p);
%
a = mt(3)*(1-2*mt(2)*mt(2)*p.*p)-mi(3)*(1-2*mi(2)*mi(2)*p.*p);
b = mt(3)*(1-2*mt(2)*mt(2)*p.*p)+2*mi(3)*mi(2)*mi(2)*p.*p;
c = mi(3)*(1-2*mi(2)*mi(2)*p.*p)+2*mt(3)*mt(2)*mt(2)*p.*p;
d = 2*(mt(3)*mt(2)*mt(2)-mi(3)*mi(2)*mi(2));
%
E = b .* etaai + c .* etaat;
F = b .* etabi + c .* etabt;
G = a - d * etaai .* etabt;
H = a - d * etaat .* etabi;
D = E.*F + G.*H.*p.*p;
%
Rpp = ( (b.*etaai-c.*etaat).*F - (a + d*etaai.*etabt).*H.*p.*p)./D;
Rps = -(2 * etaai .* (a .* b + d * c .* etaat .* etabt) .* p * mi(1)/mi(2) )./D;
Rss = -((b.*etabi-c.*etabt).*E-(a+d.*etaat.*etabi).*G.*p.*p)./D;
Rsp = -(2*etabi.*(a.*b+d*c.*etaat.*etabt).*p*(mi(2)/mi(1)))./D;
Tpp =  (2*mi(3)*etaai.*F*(mi(1)/mt(1)))./D;
Tps =  (2*mi(3)*etaai.*H.*p*(mi(1)/mt(2)))./D;
Tss =   2*mi(3)*etabi.*E*(mi(2)/mt(2))./D;
Tsp =  -2*(mi(3)*etabi.*G.*p*(mi(2)/mt(1)))./D;
%
RTmatrix = [ Rpp' Rps' Rss' Rsp' Tpp' Tps' Tss' Tsp'];
%

An Example

Here is an example call of the script that also generates plots similar to those in Lay and Wallace (1995), Figure 3.28 on page 105.

% set up the velocity model for two half-spaces
%
%     Vp   Vs   Density
mt = [8.00, 4.6, 3.38 ];  % mt is the half space with transmitted waves
mi = [4.98, 2.9, 2.667];  % mi is the half space with incident wave
%
% generate the ray parameter array
%   the value goes from 0 to 1/Vp(incident)
%
n  = 200;
p  = (0:n)*(1/mi(1))/n;
%
R = PSVRTmatrix(p,mi,mt);
%
Rpp = R(:,1); % P-to-P reflection
Rps = R(:,2); % P-to-S reflection
Tpp = R(:,5); % P-to-P transmission
Tps = R(:,6); % P-to-S transmission
%
%  The rest is graphics
%
% compute incidence angle in degrees
%
angle = (180/pi) * (asin(p*mi(1)));
%
ymin = -.1; ymax = 2.0;
xmin = 0; xmax = 90;
%
subplot('position',[0.1 0.6, 0.35 0.3])';
plot(angle,abs(Rpp),'k-')
xlabel('Incidence Angle (i)');
ylim([ymin ymax]); xlim([xmin xmax]); grid on;
title('Reflected P');
%
subplot('position',[0.55 0.6, 0.35 0.3])';
plot(angle,abs(Rps),'k-')
xlabel('Incidence Angle (i)');
ylim([ymin ymax]); xlim([xmin xmax]); grid on;
title('Reflected S');
%
subplot('position',[0.1 0.1, 0.35 0.3])';
plot(angle,abs(Tpp),'k-')
xlabel('Incidence Angle (i)');
ylim([ymin ymax]); xlim([xmin xmax]); grid on;
title('Transmitted P');
%
subplot('position',[0.55 0.1, 0.35 0.3])';
plot(angle,abs(Tps),'k-')
xlabel('Incidence Angle (i)');
ylim([ymin ymax]); xlim([xmin xmax]); grid on;
title('Transmitted S');
h = gca;
%

The results of running the script are shown below. Note the large changes in the coefficients when the P-wave becomes critical.

Results of running the PSVRTmatrix script with the parameters used in Lay and Wallace (1995) Figure 3.28. The critical angle of the P-wave is 38.5°. Beyond this angle the transmission coefficient is complex.

Summary

Matlab, like Mathematica and Python Jupyter Notebooks, is an effective tool for exploring relatively simple yet often important scientific computations. Most practicing scientists use a number of these computations on a daily basis to explore data and gain insight into the theoretical basis for much scientific analysis.