Russel D., Acoustics Animations : Longitudinal and Transverse Wave Motion
The site illustrates the 1D propagation of longitudinal and transverse harmonic waves as well as Rayleigh waves guided by a free surface.
Baker S.A. and Wysession M.E., Seismic shear wave propagation : a synthetic animation
Animations of the propagation and reflection of wave fronts in the Earth's mantle,
and examples of ray tracing and recordings.
Tanimoto T., Southern California wavefield reconstruction
Visualisation of the propagation of wave fields recorded by a dense array of seismic stations.
.  elastic modulus  dilatation  shear  volumetric 
strain  .  ε_{xx} = ∂u_{x}/∂x  ε_{xy} = (∂u_{x}/∂y+∂u_{y}/∂x)/2  divu = ε_{xx}+ε_{yy}+ε_{zz} 
uniaxial compression  Young modulus E
Poisson coefficient ν  σ_{xx} = Eε_{xx}
ε_{yy} = νε_{xx}  .  divu = (12ν)ε_{xx} 
stressstrain relation  coef. Lamé λ , μ
bulk modulus K  σ_{xx} = λdivu+2με_{xx}  σ_{xy} = 2με_{xy}  divu = (σ_{xx}+σ_{yy}+σ_{zz})/3K 
relationships among moduli  E = 10^{10}  10^{11} Pa
ν = 0.2  0.4 (rocks) ; 0.5 (water)  λ = Eν/((1+ν)(12ν))  μ = E/(2(1+ν))  K = λ+2μ/3 = E/(3(12ν)) 
Values of elastic moduli et densities for air, water and some minerals are the following :
air  water  ice  clays  calcite CaCO3  quartz SiO2  olivine (Mg,Fe)SiO4  
K (GPa)  0.0001  2.4  8.4  25  70  37  130 
μ (GPa)  0  0  3.6  9  29  44  80 
ρ (kg/m^{3})  1  1000  920  2550  2700  2700  3320 
In cylindrical coordinates (r,θ,z), the strain and stress components are :
strain  ε_{rr} = ∂u_{r}/∂r  ε_{θθ} = (1/r)(u_{r}+∂u_{θ}/∂θ)  ε_{zz} = ∂u_{z}/∂z  ε_{zr} = (∂u_{r}/∂z+∂u_{z}/∂r)/2  ε_{rθ} = ((1/r)(∂u_{r}/∂θu_{θ})+∂u_{θ}/∂r)/2  ε_{θz} = ((1/r)∂u_{z}/∂θ+∂u_{θ}/∂z)/2 
stress  σ_{rr} = λdivu+2με_{rr}  σ_{θθ} = λdivu+2με_{θθ}  σ_{zz} = λdivu+2με_{zz}  σ_{zr} = 2με_{zr}  σ_{rθ} = 2με_{rθ}  σ_{θz} = 2με_{θz} 
Compressional waves  1D longitudinal waves  3D acoustic waves  P waves 
particle motion  u_{x}(x,t)  u = gradφ(x,t)  u = gradφ(x,t) 
strains  ε_{xx} = ∂u_{x}/∂x  divu = Δφ  ε_{xx} = ∂^{2}φ/∂x^{2} , ε_{yy} , ε_{zz}
ε_{xy} = ∂^{2}φ/∂x∂y , ε_{yz} , ε_{zx} 
stresses  σ_{xx} = E∂u_{x}/∂x  p = Kdivu  σ_{xx} = λΔφ+2μ∂^{2}φ/∂x^{2} , σ_{yy} , σ_{zz}
σ_{xy} = 2μ∂^{2}φ/∂x∂y , σ_{yz} , σ_{zx} 
equilibrium equations  ρ∂^{2}u_{x}/∂t^{2} = ∂σ_{xx}/∂x  ρ∂^{2}u/∂t^{2} = gradp  ρ∂^{2}u_{x}/∂t^{2} = ∂σ_{xx}/∂x+∂σ_{xy}/∂y+∂σ_{xz}/∂z
ρ∂/∂x(∂^{2}φ/∂t^{2}) = (λ+2μ)∂/∂x(Δφ) ρgrad(∂^{2}φ/∂t^{2}) = (λ+2μ)grad(Δφ) 
wave equations  ∂^{2}u_{x}/∂t^{2} = V_{L}^{2}∂^{2}u_{x}/∂x^{2}  ∂^{2}p/∂t^{2} = V_{A}^{2}Δp  ∂^{2}φ/∂t^{2} = V_{P}^{2}Δφ 
propagation velocity  V_{L}^{2} = E/ρ  V_{A}^{2} = K/ρ  V_{P}^{2} = (λ+2μ)/ρ 
For S waves propagating in a vertical plane (it is the case when the propagation velocity is constant or when it depends only on depth, as in a horizontally stratified medium), it is convenient to distinguish two types of transverse polarizations. SH waves have an horizontal polarization orthogonal to the propagation plane : the particle motion has a single non zero component. SV waves have a transverse polarization within the propagation plane : the two non zero components of the particle motion are obtained from a vector potential with a single nonzero component, the one that is orthogonal to the propagation plane.
Shear waves (divu = 0)  1D transverse wave  SH waves (polarization ⊥ propagation plane)  SV waves (polarization within the propagation plane) 
particle motion  u_{y}(x,t)  u_{y}(x,z,t)  u(x,z,t) = curlψ_{y}(x,z,t)
u_{x} = ∂ψ_{y}/∂z , u_{z} = ∂ψ_{y}/∂x 
strains  ε_{xy} = (1/2)∂u_{y}/∂x  ε_{xy} = (1/2)∂u_{y}/∂x , ε_{yz}  ε_{xx} = ∂^{2}ψ_{y}/∂z∂x = ε_{zz}
ε_{xz} = (∂^{2}ψ_{y}/∂z^{2}+∂^{2}ψ_{y}/∂x^{2})/2 
stresses  σ_{xy} = μ∂u_{y}/∂x  σ_{xy} = μ∂u_{y}/∂x , σ_{yz}  σ_{xx} = 2μ∂^{2}ψ_{y}/∂z∂x = σ_{zz}
σ_{xz} = μ(∂^{2}ψ_{y}/∂z^{2}+∂^{2}ψ_{y}/∂x^{2}) 
equilibrium equations  ρ∂^{2}u_{y}/∂t^{2} = ∂σ_{xy}/∂x  ρ∂^{2}u_{y}/∂t^{2} = ∂σ_{xy}/∂x+∂σ_{yz}/∂z  ρ∂^{2}u_{x}/∂t^{2} = ∂σ_{xx}/∂x+∂σ_{xz}/∂z
ρ∂/∂z(∂^{2}ψ_{y}/∂t^{2}) = μ∂/∂z(∂^{2}ψ_{y}/∂x^{2}+∂^{2}ψ_{y}/∂z^{2}) 
wave equations  ∂^{2}u_{y}/∂t^{2} = V_{S}^{2}∂^{2}u_{y}/∂x^{2}  ∂^{2}u_{y}/∂t^{2} = V_{S}^{2}Δu_{y}  ∂^{2}ψ_{y}/∂t^{2} = V_{S}^{2}Δψ_{y} 
propagation velocity  V_{S}^{2} = μ/ρ  V_{S}^{2} = μ/ρ  V_{S}^{2} = μ/ρ 
Values of propagation velocities, Poisson coefficients and densities for some sedimentary and crystalline rocks are the following :
unconsolidated sandstone  consolidated sandstone  limestone  granite  basalt  granulite  dunite  
V_{P} (km/s)  2.7  4.1  4.6  6.2  5.9  7.0  8.3 
V_{S} (km/s)  1.4  2.4  2.4  3.7  3.2  3.8  4.8 
ν = (V_{P}^{2}/V_{S}^{2}2)/(2V_{P}^{2}/V_{S}^{2}2)  0.32  0.24  0.31  0.22  0.29  0.29  0.25 
ρ (kg/m^{3})  2100  2400  2400  2650  2880  3000  3300 
1D wave  P waves  S waves 
u→ = u(x/Vt)
u← = u(x/Vt)  φ(e.x/V_{P}t)
u_{longitudinal} = gradφ = (e/V_{P})φ'  ψ(e.x/V_{S}t)
u_{transverse} = curlψ = (e/V_{S})∧ψ' 
harmonic case : u = Ue^{iω(x/Vt)} = Ue^{i2π(x/λt/T)} = Ue^{i(kxωt)}
u→ : k = ω/V u← : k = ω/V  φ = Φe^{i(k.xωt)} , k = (ω/V_{P})e
u = ikΦe^{i(k.xωt)}  ψ = Ψe^{i(k.xωt)} , k = (ω/V_{S})e
u = ik∧Ψe^{i(k.xωt)} 
For harmonic waves ∂φ/∂x = ik_{x}φ , ∂φ/∂t = iωφ .
The wave equation becomes the dispersion relation : ω^{2} = (k_{x}^{2}+k_{y}^{2}+k_{z}^{2})V^{2} .
It is a sphere of radius 1/V in the coordinates (k_{x}/ω , k_{y}/ω , k_{z}/ω).
This is because in a homogeneous isotropic medium, the propagation velocity is the same in all propagation directions.
An harmonic wave fills the space with periodic oscillations.
If one looks at it or records it along its propagation direction e, one sees its true wavelength λ = 2π/k.
Along a different direction e', one sees an apparent wavelength larger than the true wavelength :
λ_{ae'} = λ/cos(angle(e,e')).
Along the x,y,z axes, one sees the apparent wavelengths λ_{x} = 2π/k_{x} , λ_{y} = 2π/k_{y} , λ_{z} = 2π/k_{z}.



ondePS.m
∂^{2}(rφ)/∂t^{2} = V^{2}∂^{2}(rφ)/∂r^{2}  particle motion  near/far field 
φ = (A/r)f(±r/Vt)  u = (A/r)e_{r}(f/r±f'/V)  . 
φ = (A/r)e^{i(krωt)}  u = (A/r)e_{r}(ik1/r)e^{i(krωt)}  kr>>1 u = (ikA/r)e_{r}e^{i(krωt)}
kr<<1 u = (A/r^{2})e_{r}e^{i(krωt)} A = U_{0}r_{0}^{2} for a compressional source with radius r_{0}<<λ 
Δφ(r) = (1/r)∂^{2}(rφ)/∂r^{2} , divr = ∂x/∂x+∂y/∂y+∂z/∂z = 3 , dive_{r} = div(r/r) = grad(1/r).r+divr/r = 2/r
.  wavefront equation  inverse of apparent velocities along axes  eikonal equation : Pythagore for waves 
plane  T(x) = e.x/V  gradT = e/V  e.e=1
gradT.gradT = (∂T/∂x)^{2}+(∂T/∂y)^{2}+(∂T/∂z)^{2} = 1/V^{2} (1/V_{ax})^{2}+(1/V_{ay})^{2}+(1/V_{az})^{2} = 1/V^{2} 
spherical  T(x) = r/V  gradT = e_{r}/V  
any surface  T(x)  gradT = e_{M}/V 
The propagation of wavefronts obeys Huygens principle : the wavefront at time t+dt is the envelop of spherical wavelets of radii Vdt, centered on the wavefront at time t.
The geometry of rays obeys Fermat principle : the traveltime along the ray is stationary relative to any perturbation of the ray geometry.
These two principles illustrate how a wave finds its direction propagation : it tries them all and chooses those where constructive interferences occur between neighbouring paths.
They can also be numerically implemented to determine the rays geometry of the rays in nonhomogeneous media.
Traveltime curves T(X) of direct, reflected and head waves for sourcereceiver offset SR=X and homogeneous media of velocity V_{1} and V_{2}  
dip α  direct wave  reflected wave
total reflection for incidence θ_{1}>θ_{c}  head wave (only if V_{2}>V_{1})
critical incidence θ_{c} : sinθ_{c} = V_{1}/V_{2} 
direco.m  T = X/V_{1}  T = (4d^{2}+X^{2})^{½}/V_{1}  critical distance : X_{c} = 2dtgθ_{c}
apparent horizontal et vertical velocities : V_{ax} = V_{1}/sinθ_{c} = V_{2} ; V_{az} = V_{1}/cosθ_{c} traveltime : T = X/V_{ax}+2d/V_{az} = X/V_{2}+2dcosθ_{c}/V_{1} 
T = X/V_{1}  virtual source symmetrical of the source with respect to interface
x_{v} = 2dsinα ; z_{v} = 2dcosα T = (4(dcosα)^{2}+(X±2dsinα)^{2})^{½}/V_{1}
 critical distance (triangle virtual source, vertical projection of VS on surface, receiver) :
X_{c} = 2d(cosαtg(θ_{c}±α)±sinα) = 2dsinθ/cos(θ±α) apparent velocities along and orthogonally to interface :
distances along and orthogonally to interface :
traveltime : T = x_{α}/V_{aα}+d_{⊥α}/V_{a⊥α}

If the sourcereceiver direction is not parallel to the dip azimuth, the propagation occurs in the plane that is orthogonal to the interface and that includes the SR line. α is the apparent dip in the propagation plane. It is related to the true dip α_{true} by the equation : sinα = sin(α_{true})cos(azimuth(RS)azimuth(dip)).
In the elastic case, converted waves have to be taken into account. Here are the figures drawn by Cagniard in 1939 in the case where V_{P2}>V_{S2}>V_{P1}>V_{S1} :
For a point source located on the free surface of a layer having a thickness d above a halfspace, the multiple reflections between the free surface and the interface generate reverberations (see also this quiz illustrating how these reverberations give rise to guided waves).
Ocean Acoustics Library shows an animation of the propagation of direct, reflected, transmitted and head waves that illustrates the physics of the interaction of the spherical wave with a plane interface. It includes the evanescent wave in the rapid medium created by the incident wave at incidences corresponding to total reflection.
Heelan P.A., 1953, On the theory of head waves
Mitchell J.F. and R.J. Bolander, 1986, Structural interpretation using refraction velocities from marine seismic surveys
show a nice recording at sea that can be used to determine the velocity V_{1} of water, the velocity V_{2} and tickness d of a sedimentary layer, and the velocity V_{3} beneath
this layer from direct and head waves. The reflected wavefield is unclear because of reverberations in the water layer. A clear reflection on a dipping interface is apparent on the lower part of the record.
Allen, R., Fratta D., Introduction to applied geophysics, offer detailed lecture notes on the use of head waves in the seismic refraction method.
If the propagation velocity depends only on depth (stratified medium with horizontal interfaces or vertical velocity gradient), the geometry of the rays is determined by Snell's law :
the ray parameter p, inverse of the horizontal apparent velocity (also called horizontal slowness), is constant.
The incidence angle θ(z) of the ray relative to the vertical (which is the common normal to all interfaces or isovelocity curves) increases or decreases with increases or decreases in velocity.
The horizontal distance X(p) travelled along a ray and the traveltime T(p) are expressed as a function of the parameter p by summing over depth intervals corresponding to layer thicknesses
in a stratified medium.
It is also useful to express τ(p) = T(p)pX(p), which represents the vertical traveltime in the medium at the apparent vertical velocities corresponding to the parameter p.
In a medium with a velocity gradient, the discrete summation is replaced by an integration over infinitesimal depth intervals dz.
If the velocity gradient a is a constant (V(z) = V_{0}+az), the rays are circles centered at the depth z = V_{0}/a where V(z) = 0.
In this case, the wavefronts are also circles. They are the circles orthogonal to the family of ray circles, all centered on the Oz axis (for a source placed at the origin), at depths increasing with time.
In the limit of small incidence angles, T(X) can be approximated with an hyperbolic expression analogous to that in a medium of constant velocity.
The velocity of the equivalent medium is V_{RMS}. This approximation is commonly used in the processing of seismic reflection data.
Snell's law : p = sinθ/V = (1/V_{ax}) = cte for a ray  ray parameter  horizontal distance X(p) between two points of a same ray of parameter p  traveltime T(p) between two points of a same ray of parameter p
vertical traveltime τ(p) at vertical apparent velocity 
stratified medium with plane horizontal interfaces (layers with thicknesses d_{i} and constant velocities V_{i})
 hodostratif.m  hodostrat.m  p = sinθ_{i}/V_{i}  X(p) = ∑d_{i}tanθ_{i} = p∑d_{i}V_{i}/(1p^{2}V_{i}^{2})^{½}  T(p) = ∑d_{i}/cosθ_{i}V_{i} = ∑d_{i}/(1p^{2}V_{i}^{2})^{½}V_{i}
τ(p) = T(p)pX(p) = ∑d_{i}cosθ_{i}/V_{i} = ∑d_{i}(1p^{2}V_{i}^{2})^{½}/V_{i} 
velocity gradient V(z)
for a positive gradient, the ray reaches a maximum depth Z_{max}(p)  p = sinθ(z)/V(z) = 1/V(Z_{max})  X(p) = p∫dzV(z)/(1p^{2}V(z)^{2})^{½}  T(p) = ∫dz/(1p^{2}V(z)^{2})^{½}V(z)
τ(p) = ∫dz(1p^{2}V(z)^{2})^{½}/V(z) 
case V(z) = V_{0}+az
the rays are circles of radius R(p)=1/(ap) centered at z = V_{0}/a where V(z) = 0  sinθ(z) = (V_{0}/a+z)/R = pV(z)  X(p) = R[cosθ] = R[(1p^{2}V(z)^{2})^{½}]  T(p) = ∫Rdθ/V(z(θ)) = (1/a)∫dθ/sinθ = (1/a)[Log(tan(θ/2))] = (1/a)[Log(pV(z)/(1+(1p^{2}V(z)^{2})^{½}))] 
case V(z) = V_{0}+az
the wavefronts at times T = (1/a)[Log(tan(θ_{0}/2))] are the circles orthogonal to the rays the equation of the wavefronts is x^{2}+(zZ_{max}(T))^{2} = X(T)^{2} hodovz.m  tan(θ_{0}/2) = e^{aT}  X(T) = V_{0}/atanθ_{0} = (V_{0}/2a)(e^{aT}e^{aT})
X(T) = (V_{0}/a)sinh(aT)  Z_{max}(T) = (V_{0}/a)(1/sinθ_{0}1) = (V_{0}/2a)(e^{aT}+e^{aT}2)
Z_{max}(T) = (V_{0}/a)(cosh(aT)1) 
V_{RMS} approximation for small incidence angles
vrms.m  θ ≈ 0 p ≈ θ_{i}/V_{i}  X(p) ≈ p∑d_{i}V_{i}  T(p) ≈ ∑d_{i}(1+p^{2}V_{i}^{2}/2)/V_{i} = ∑d_{i}/V_{i}+(p^{2}/2)∑d_{i}V_{i} = T(X=0)+X^{2}/2∑d_{i}V_{i} 
T(X)^{2} ≈ T(0)^{2}+X^{2}/V_{RMS}^{2}
V_{RMS}^{2} = (∑d_{i}V_{i})/(∑d_{i}/V_{i}) = (∫V(z)dz)/(∫dz/V(z)) V_{RMS}(V_{0}+az) = (a(V_{0}z+az^{2}/2)/Log(1+az/V_{0}))½ 
V(z) generally increases with depth. The rays originating from a point source are curved toward the surface and reach it at increasing horizontal distances
for incidence angles θ_{0} decreasing from 90° (ray departing horizontally from the source) to 0° (ray departing vertically from the source).
The wavefront reaches the surface with an increasing apparent horizontal velocity.
If the velocity gradient has a large increase in a depth interval, the traveltime curve shows a triplication.
If the velocity gradient has a sign change in a depth interval, there is a shadow zone where no geometrical ray penetrates. In the corresponding X interval, no ray reaches the surface.
If the source is situated in a low velocity zone, the rays that bottom and top in the zone never leave it. The low velocity zone acts as a waveguide.
A plane wavefront incident on a depth interval where V(z) increases with depth is totally reflected
if the bottom point of the rays is in the interval.
Lutter W.J., R. D. Catchings and C. M. Jarchow, 1994,
An image of the Columbia Plateau from inversion of highresolution seismic data show the effect of a stratification basalt (5 km/s)  sediment (4 km/s)  basement (5.5 km/s) on first arrivals.
Hornby B., 1993,
Tomographic reconstruction of nearborehole slowness using refracted borehole sonic arrivals is an example of the use of head and refracted waves in seismic logging.
Russel D., Refraction of sound in the atmosphere
Garcés, M.A., Hansen, R. A. & Lindquist, K.G., 1998,
Traveltimes for infrasonic waves propagating in a stratified atmosphere
e = cosθ_{1}n+sinθ_{1}m  e_{r} = cosθ_{1}n+sinθ_{1}m  e_{r}e = 2cosθ_{1}n = 2(e.n)n  e_{r} = e2(e.n)n 
e/V_{1} = cosθ_{1}/V_{1}n+sinθ_{1}/V_{1}m  e_{t}/V_{2} = cosθ_{2}/V_{2}n+sinθ_{2}/V_{2}m  e_{t}/V_{2}e/V_{1} = (cosθ_{1}/V_{1}cosθ_{2}/V_{2})n  e_{t} = V_{2}/V_{1}(e(e.n+((V_{1}/V_{2})^{2}(1e.n^{2}))^{½})n) 
Here is an example for a medium including two dipping interfaces. The rays displayed are PP and PS reflections on the deepest interface. One can verify that the reflection points correspond to an extremum of the traveltimes.


refl3d.m
In the case of continuous variations of the propagation velocity with depth and lateral position, both the horizontal and vertical apparent velocities vary along a given ray.
The variation along a ray of the inverse of the apparent velocity in a given direction is equal to the component of the gradient of 1/V (called slowness) in this direction.
This relationship allows tracing rays by numerical integration with respect to the arclength along the ray and computing at each integration step the new propagation direction e.
In the case where V(x,z) varies linearly with x and z, the seismic rays are circles centered on the straigth line of equation V(x,z) = 0.
position x(s) along the ray as a function of arclength s
e = dx/ds  differential equation of the rays
e.e = 1 , e.∂e/∂x = 0 , gradT.∂e/∂x = 0 ∂(dT/ds)/∂x = ∂(e.gradT)/∂x = e.grad∂T/∂x = d(∂T/∂x)/ds  curvature of rays de/ds = e_{N}/R  torsion of rays
b = e∧e_{N} de_{N}/ds = e/Rb/τ 
gradT = e/V
T = ∫ds/V(x(s))  d(e/V)/ds = d(gradT)/ds = grad(dT/ds) = grad(1/V)  1/R = Ve_{N}.grad(1/V)  1/τ = Rb.d(Vgrad(1/V))/ds 
case V(x,z) = V_{0}+a_{x}x+a_{z}z = V_{0}+a(xsinα+zcosα)
by a rotation with angle α, one gets back to the case V(Z) = V_{0}+aZ seismic rays are circles centered on the straigth line of equation V(x,z) = 0 
d(sinθ/V)/ds = a_{x}/V^{2} = asinα/V^{2}
d(cosθ/V)/ds = a_{z}/V^{2} = acosα/V^{2}  1/R = dθ/ds = a_{x}cosθ/V+a_{z}sinθ/V = asin(θα)/V 
One can derive ray equations directly from the wave equation for a V(x) medium by looking for solutions expressed as quasiplane harmonic waves with wavefront T(x) and amplitude A(x). One gets terms depending on successive powers of the angular frequency ω. The term in ω^{2}, dominant at high frequency, gives the eikonal equation of traveltimes. It can be solved numerically in heterogeneous media by a finite difference method. The term in ω gives the transport equation that determines the variation of the amplitude A(x).
wave equation  quasiplane wave  term in ω^{2} → eikonal equation  term in ω → transport equation  ray tube with section dS(s) 
Δφ+(ω/V(x))^{2}φ = 0  φ = A(x)e^{iω(T(x)t)}  gradT.gradT=1/V^{2}  2gradA.gradT+AΔT = div(A^{2}gradT) = 0  ∫∫∫div(A^{2}gradT)dυ = ∫∫(A^{2}/V)e.dσ = 0
A^{2}dS/V = cte 
Langan, R.T., I. Lerche, and R. T. Cutler , 1985, Tracing of rays through heterogeneous media: An accurate and efficient procedure
Bernasconi G. and G. Drufuca, 2001, 3D traveltimes and amplitudes by gridded rays
give a simple method of 3D raytracing by a finitedifference method.
Sava P. and S. Fomel, M. 2001, 3D traveltime computation using Huygens wavefront tracing
ω and k_{x} = ωp = ωsin(θ)/V_{P} = ωsin(η)/V_{S} of the incident wave are the same for all reflected and transmitted waves on a plane horizontal interface  
downgoing waves :
k_{z}^{P} = +ωcos(θ)/V_{P} , k_{z}^{S} = +ωcos(η)/V_{S}
upgoing waves : k_{z}^{P} = ωcos(θ)/V_{P} , k_{z}^{S} = ωcos(η)/V_{S}  SH waves  PSV waves 
free surface  σ_{zy}^{Incident SH}+σ_{zy}^{Reflected SH} = 0 
σ_{zz}^{Incident P or SV}+σ_{zz}^{Reflected P}+σ_{zz}^{Reflected SV} = 0
σ_{zx}^{Incident P or SV}+σ_{zx}^{Reflected P}+σ_{zx}^{Reflected SV} = 0 
solid  solid interface 
u_{y }^{Incident SH}+u_{y }^{Reflected SH} = u_{y }^{Transmitted SH}
σ_{zy}^{Incident SH}+σ_{zy}^{Reflected SH} = σ_{zy}^{Transmitted SH} 
u_{x}^{Incident P or SV}+u_{x}^{Reflected P}+u_{x}^{Reflected SV} =
u_{x}^{Transmitted P}+u_{x}^{Transmitted SV}
u_{z}^{Incident P or SV}+u_{z}^{Reflected P}+u_{z}^{Reflected SV} = u_{z}^{Transmitted P}+u_{z}^{Transmitted SV} σ_{zz}^{Incident P or SV}+σ_{zz}^{Reflected P}+σ_{zz}^{Reflected SV} = σ_{zz}^{Transmitted P}+σ_{zz}^{Transmitted SV} σ_{zx}^{Incident P or SV}+σ_{zx}^{Reflected P}+σ_{zx}^{Reflected SV} = σ_{zx}^{Transmitted P}+σ_{zx}^{Transmitted SV} 
liquid  solid interface  . 
u_{z}^{Incident P}+u_{z}^{Reflected P} =
u_{z}^{Transmitted P}+u_{z}^{Transmitted SV}
σ_{zz}^{Incident P}+σ_{zz}^{Reflected P} = σ_{zz}^{Transmitted P}+σ_{zz}^{Transmitted SV} σ_{zx}^{Transmitted P}+σ_{zx}^{Transmitted SV} = 0 
Snell's law :
k_{x} = ωp = ωsinη_{1}/V_{S1} = ωsinη_{2}/V_{S2}  incident downgoing SH wave  reflected upgoing SH wave  transmitted downgoing SH wave  
particle motion SH wave
u_{y}(x,z,t) = U_{y}e^{ikzz}e^{i(kxxωt)} 
U_{y}^{I} = 1 k_{z}^{I} = ωcosη_{1}/V_{S1} 
U_{y}^{R} = R k_{z}^{R} = ωcosη_{1}/V_{S1} 
U_{y}^{T} = T k_{z}^{T} = ωcosη_{2}/V_{S2}  
stress on horizontal surface
σ_{zy}(x,z,t) = μ∂u_{y}/∂z = ρV_{S}^{2}ik_{z}u_{y} = iωS_{zy}e^{ikzz}e^{i(kxxωt)}  S_{zy}^{I} = ρ_{1}V_{S1}cosη_{1}  S_{zy}^{R} = ρ_{1}V_{S1}cosη_{1}R  S_{zy}^{T} = ρ_{2}V_{S2}cosη_{2}T  
continuity on interface z =0
(u_{y}^{I}+u_{y}^{R})(z=0) = u_{y}^{T}(z=0) (σ_{zy}^{I}+σ_{zy}^{R})(z=0) = σ_{zy}^{T}(z=0) 
1+R=T
(1R)ρ_{1}V_{S1}cosη_{1} = Tρ_{2}V_{S2}cosη_{2} 
R = (ρ_{1}V_{S1}cosη_{1}ρ_{2}V_{S2}cosη_{2})/(ρ_{1}V_{S1}cosη_{1}+ρ_{2}V_{S2}cosη_{2})
T = 2ρ_{1}V_{S1}cosη_{1}/(ρ_{1}V_{S1}cosη_{1}+ρ_{2}V_{S2}cosη_{2})  
weak contrast
ρ+δρ , V+δV δη = tanηδV/V (δp=0)  R ≈ δ(ρ_{}V_{}cosη_{})/2ρ_{}V_{}cosη_{} = (1/2)(δρ/ρ+(1tg^{2}η)δV/V)  T ≈ 1  
total reflection
sinη_{c}=V_{S1}/V_{S2} η_{1}>η_{c} , pV_{S2}>1 
k_{z}^{R} = k_{z}^{I} = ω(1p^{2}V_{S1}^{2})^{½}/V_{S1}
R = (ρ_{1}V_{S1}(1p^{2}V_{S1}^{2})^{½}iρ_{2}V_{S2}(p^{2}V_{S2}^{2}1)^{½})/(ρ_{1}V_{S1}(1p^{2}V_{S1}^{2})^{½}+iρ_{2}V_{S2}(p^{2}V_{S2}^{2}1)^{½}) R = e^{iχ(p)} , χ(p) = 2Arctan(ρ_{2}V_{S2}(p^{2}V_{S2}^{2}1)^{½}/ρ_{1}V_{S1}(1p^{2}V_{S1}^{2})^{½}) (u_{y}^{I}+u_{y}^{R})(x,z,t) = (e^{ikzIz}+e^{i(kzRz+χ(p))})e^{i(kxxωt)} = 2cos(k_{z}^{I}zχ(p)/2)e^{i(kxxωt+χ(p)/2)} 
transmitted evanescent wave : k_{z}^{T} = iω(p^{2}V_{S2}^{2}1)^{½}/V_{S2}
T = 2ρ_{1}V_{S1}(1p^{2}V_{S1}^{2})^{½}/(ρ_{1}V_{S1}(1p^{2}V_{S1}^{2})^{½}+iρ_{2}V_{S2}(p^{2}V_{S2}^{2}1)^{½}) T = Te^{iχ(p)/2} u_{y}^{T}(x,z,t) = Te^{ωz(p2VS221)½/VS2}e^{i(kxxωt+χ(p)/2)} 
Here are the reflection and transmission coefficients and the particle motion computed as a function of incidence angle in the case of partial and total reflection :

reflsh.m
Here is an animation in the case of a reflection on a free surface (R=1) that produces standing waves along depth with nodes (zero) of stress for z = nπ/k_{z} and antinodes (maxima) at z = (n+1/2)π/k_{z}.
Ocean Acoustics Library : total reflection of a plane wave on a plane interface
Snell's law : k_{x} = ωp = ωsinθ_{1}/V_{P}_{1} = ωsinη_{1}/V_{S}_{1} = ωsinθ_{2}/V_{P}_{2} = ωsinη_{2}/V_{S}_{2}  
P incident downgoing Φ^{I} = 1 k_{z}^{I} = ωcosθ_{1}/V_{P1}  P reflected upgoing Φ^{R} = R_{φ} k_{z}^{P} = ωcosθ_{1}/V_{P1}  SV reflected upgoing Ψ^{R} = R_{ψ} k_{z}^{S} = ωcosη_{1}/V_{S1}  P transmitted downgoing Φ^{T} = T_{φ} k_{z}^{P} = ωcosθ_{2}/V_{P2}  SV transmitted downgoing Ψ^{T} = T_{ψ} k_{z}^{S} = ωcosη_{2}/V_{S2} 
U_{x} = p  + pR_{φ}  + (cosη_{1}/V_{S1})R_{ψ}  = pT_{φ}  (cosη_{2}/V_{S2})T_{ψ} 
U_{z} = cosθ_{1}/V_{P1}  (cosθ_{1}/V_{P1})R_{φ}  + pR_{ψ}  = (cosθ_{2}/V_{P2})T_{φ}  + pT_{ψ} 
S_{zz} = ρ_{1}(12p^{2}V_{S1}^{2})  + ρ_{1}(12p^{2}V_{S1}^{2})R_{φ}   (2ρ_{1}V_{S1}^{2}pcosη_{1}/V_{S1})R_{ψ}  = ρ_{2}(12p^{2}V_{S2}^{2})T_{φ}  + (2ρ_{2}V_{S2}^{2}pcosη_{2}/V_{S2})T_{ψ} 
S_{zx} = 2ρ_{1}V_{S1}^{2}pcosθ_{1}/V_{P1}  + (2ρ_{1}V_{S1}^{2}pcosθ_{1}/V_{P1})R_{φ}  + ρ_{1}(12p^{2}V_{S1}^{2})R_{ψ}  = (2ρ_{2}V_{S2}^{2}pcosθ_{2}/V_{P2})T_{φ}  + ρ_{2}(12p^{2}V_{S2}^{2})T_{ψ} 
matrix notation case liquidsolid  AR = B ; R = [R_{φ} , R_{ψ} , T_{φ} , T_{ψ}]^{t} ;
B = [p , cosθ_{1}/V_{P1} , ρ_{1}(12p^{2}V_{S1}^{2}) , 2ρ_{1}V_{S1}^{2}pcosθ_{1}/V_{P1}]^{t}
A_{L}R_{L} = B_{L} ; R_{L} = [R_{φ} , T_{φ} , T_{ψ}]^{t} ; B_{L} = [cosθ_{1}/V_{P1} , ρ_{1} , 0]^{t}  
A[4,4]
A_{L}[3,3]  p  cosη_{1}/V_{S1}  p  cosη_{2}/V_{S2} 
cosθ_{1}/V_{P1}  p  cosθ_{2}/V_{P2}  p  
ρ_{1}(12p^{2}V_{S1}^{2}) (ρ_{1})  2ρ_{1}V_{S1}^{2}pcosη_{1}/V_{S1}  ρ_{2}(12p^{2}V_{S2}^{2})  2ρ_{2}V_{S2}^{2}pcosη_{2}/V_{S2}  
2ρ_{1}V_{S1}^{2}pcosθ_{1}/V_{P1} (0)  ρ_{1}(12p^{2}V_{S1}^{2})  2ρ_{2}V_{S2}^{2}pcosθ_{2}/V_{P2}  ρ_{2}(12p^{2}V_{S2}^{2})  
P incident downgoing or upgoing
reflpsv.m 
R_{φ} = det([B,A[:,2:4]])/det(A)
t_{φ} = det([b,A[:,2:4]])/det(A) 
R_{ψ} = det([A[:,1],B,A[:,3:4]])/det(A)
t_{ψ} = det([A[:,1],b,A[:,3:4]])/det(A) 
T_{φ} = det([A[:,1:2],B,A[:,4]])/det(A)
r_{φ} = det([A[:,1:2],b,A[:,4]])/det(A) 
T_{ψ} = det([A[:,1:3],B])/det(A)
r_{ψ} = det([A[:,1:3],b])/det(A) 
P incident downgoing
D = det(A) = PQ R_{φ} = (PQ)/D R_{ψ} =2S/D 
AR = B ; R = [R_{φ} , R_{ψ} , T_{φ} , T_{ψ}]^{t} ;
B = [p , cosθ_{1}/V_{P1} , ρ_{1}(12p^{2}V_{S1}^{2}) , 2ρ_{1}V_{S1}^{2}pcosθ_{1}/V_{P1}]^{t} = [a_{11} , a_{21} , a_{31} , a_{41}]^{t} ; D = (a_{11}A_{11}+a_{31}A_{31})(a_{21}A_{21}+a_{41}A_{41}) ; DR_{φ} = (a_{11}A_{11}+a_{31}A_{31})(a_{21}A_{21}+a_{41}A_{41}) DR_{ψ}/2 = a_{11}a_{21}(a_{33}a_{44}a_{43}a_{34}) +a_{11}a_{41}(a_{23}a_{34}a_{33}a_{24}) +a_{31}a_{41}(a_{13}a_{24}a_{23}a_{14}) DT_{φ}/2 = (a_{11}a_{21}(a_{32}a_{44}a_{42}a_{34}) +a_{11}a_{41}(a_{22}a_{34}a_{32}a_{24}) +a_{31}a_{41}(a_{12}a_{24}a_{22}a_{14})) DT_{ψ}/2 = (a_{11}a_{21}(a_{32}a_{43}a_{42}a_{33}) +a_{11}a_{41}(a_{22}a_{33}a_{32}a_{23}) +a_{31}a_{41}(a_{12}a_{23}a_{22}a_{13})  
P incident upgoing
k_{z}^{I} = ωcosθ_{2}/V_{P2} 
Ar = b ; r = [t^{φ},t^{ψ},r^{φ},r^{ψ}]^{t}
b = [p , cosθ_{2}/V_{P2} , ρ_{2}(12p^{2}V_{S2}^{2}) , 2ρ_{2}V_{S2}^{2}pcosθ_{2}/V_{P2}]^{t} = [a_{13} , a_{23} , a_{33} , a_{43}]^{t} D = (a_{13}A_{13}+a_{33}A_{33})(a_{23}A_{23}+a_{43}A_{43}) ; Dr_{φ} = (a_{13}A_{13}+a_{33}A_{33})(a_{23}A_{23}+a_{43}A_{43}) Dr_{ψ}/2 = a_{13}a_{23}(a_{31}a_{42}a_{41}a_{32}) +a_{13}a_{43}(a_{21}a_{32}a_{31}a_{22}) +a_{33}a_{43}(a_{11}a_{22}a_{21}a_{12}) Dt_{φ}/2 = (a_{13}a_{23}(a_{34}a_{42}a_{44}a_{32}) +a_{13}a_{43}(a_{24}a_{32}a_{34}a_{22}) +a_{33}a_{43}(a_{14}a_{22}a_{24}a_{12})) Dt_{ψ}/2 = a_{13}a_{23}(a_{34}a_{41}a_{44}a_{31}) +a_{13}a_{43}(a_{24}a_{31}a_{34}a_{21}) +a_{33}a_{43}(a_{14}a_{21}a_{24}a_{11})  
case liquidsolid
reflpsl.m 
A_{L}R_{L} = B_{L} ; R_{L} = [R_{φL} , T_{φL} , T_{ψL}]^{t} ;
B_{L} = [cosθ_{1}/V_{P1} , ρ_{1} , 0]^{t}
A_{L}r_{L} = b_{L} ; r_{L} = [t_{φL} , r_{φL} , r_{ψL}]^{t} ; b_{L} = [cosθ_{2}/V_{P2} , ρ_{2}(12p^{2}V_{S2}^{2}) , 2ρ_{2}V_{S2}^{2}pcosθ_{2}/V_{P2}]^{t} D_{L} = ρ_{2}^{2}cosθ_{1}/V_{P1}((12p^{2}V_{S2}^{2})^{2}+4V_{S2}^{4}p^{2}cosθ_{2}/V_{P2}cosη_{2}/V_{S2}) ρ_{1}ρ_{2}cosθ_{2}/V_{P2} D_{L}R_{φL} = ρ_{2}^{2}cosθ_{1}/V_{P1}((12p^{2}V_{S2}^{2})^{2}+4V_{S2}^{4}p^{2}cosθ_{2}/V_{P2}cosη_{2}/V_{S2}) +ρ_{1}ρ_{2}cosθ_{2}/V_{P2} D_{L}T_{φL} = 2ρ_{1}ρ_{2}cosθ_{1}/V_{P1}(12p^{2}V_{S2}^{2}) ; D_{L}t_{φL} = 2ρ_{2}^{2}cosθ_{2}/V_{P2}(12p^{2}V_{S2}^{2}) D_{L}T_{ψL} = 4ρ_{1}ρ_{2}cosθ_{1}/V_{P1}V_{S2}^{2}pcosθ_{2}/V_{P2}  
SV incident downgoing or upgoing
reflpsv.m 
R'_{φ} = det([B',A[:,2:4]])/D
t'_{φ} = det([b',A[:,2:4]])/D 
R'_{ψ} = det([A[:,1],B',A[:,3:4]])/D
t'_{ψ} = det([A[:,1],b',A[:,3:4]])/D 
T'_{φ} = det([A[:,1:2],B',A[:,4]])/D
r'_{φ} = det([A[:,1:2],b',A[:,4]])/D 
T'_{ψ} = det([A[:,1:3],B'])/D
r'_{ψ} = det([A[:,1:3],b'])/D 
SV incident downgoing
Ψ^{I} = 1 ; k_{z}^{I} = ωcosη_{1}/V_{S1} D = det(A) = P'+Q' R'_{ψ} = (P'Q')/D 
AR'=B' ; R' = [R'_{φ} , R'_{ψ} , T'_{φ} , T'_{ψ}]^{t} ;
B' = [cosη_{1}/V_{S1} , p , 2ρ_{1}V_{S1}^{2}pcosη_{1}/V_{S1} , ρ_{1}(12p^{2}V_{S1}^{2})]^{t} = [a_{12} , a_{22} , a_{32} , a_{42}]^{t} ; D = (a_{12}A_{12}+a_{32}A_{32})+(a_{22}A_{22}+a_{42}A_{42}) ; DR'_{ψ} = (a_{12}A_{12}+a_{32}A_{32})(a_{22}A_{22}+a_{42}A_{42})  
SV incident upgoing
k_{z}^{I} = ωcosη_{2}/V_{S2} 
Ar'=b' ; r' = [t'_{φ} , t'_{ψ} , r'_{φ} , r'_{ψ}]^{t}
b' = [cosη_{2}/V_{S2} , p , 2ρ_{2}V_{S2}^{2}pcosη_{2}/V_{S2} , ρ_{2}(12p^{2}V_{S2}^{2})]^{t} = [a_{14} , a_{24} , a_{34} , a_{44}]^{t} ; D = (a_{14}A_{14}+a_{34}A_{34})+(a_{24}A_{24}+a_{44}A_{44}) Dr'_{ψ} = (a_{14}A_{14}+a_{34}A_{34})(a_{24}A_{24}+a_{44}A_{44})  
total reflection of downgoing SV wave
p = sinη_{1}/V_{S}_{1} 1/V_{S1} > p > 1/V_{P1} p > 1/V_{S2} > 1/V_{P2} R'_{ψ} = e^{iχ'(p)} 
imaginary terms of A :
a_{14} = i(p^{2}V_{S2}^{2}1)^{½}/V_{S2} ;
a_{21} = i(p^{2}V_{P1}^{2}1)^{½}/V_{P1} ;
a_{23} = i(p^{2}V_{P2}^{2}1)^{½}/V_{P2} ;
a_{34} = 2iρ_{2}V_{S2}^{2}p(p^{2}V_{S2}^{2}1)^{½}/V_{S2} ; a_{41} = 2iρ_{1}V_{S1}^{2}p(p^{2}V_{P1}^{2}1)^{½}/V_{P1} ; a_{43} = 2iρ_{2}V_{S2}^{2}p(p^{2}V_{P2}^{2}1)^{½}/V_{P2} ; ⇒ A_{12}, A_{32} and P' are imaginary, A_{12}, A_{32} and Q' are real R'_{ψ} = (P'Q')/(P'+Q') = e^{iχ'(p)} χ'(p) = 2Arctan(P'/iQ') = 2Arctan((1/i)(a_{12}A_{12}+a_{32}A_{32})/(a_{22}A_{22}+a_{42}A_{42})) 
Here are the moduli of the reflectiontransmission coefficients PSV and SVP computed for incidence angles of 0°, 15° and 30° in water for a sedimentary stratified medium :

weak contrasts 
ρ_{2} = ρ+δρ ;
V_{P}_{2} = V_{P}+δV_{P} ;
θ_{2} = θ+δθ ; δθ = tanθδV_{P}/V_{P} ;
V_{S}_{2} = V_{S}+δV_{S} ;
η_{2} = η + δη ; δη = tanηδV_{S}/V_{S}
a_{11} = a_{13} = a_{22} = a_{24} = p ; a_{12} = cosη/V_{S} ; a_{21} = cosθ/V_{P} ; a_{31} = a_{42} = ρ(12p^{2}V_{S}^{2}) ; a_{32} = 2ρV_{S}^{2}pcosη/V_{S} ; a_{41} = 2ρV_{S}^{2}pcosθ/V_{P} a_{14} = (a_{12}+δa_{12}) ; a_{23} = (a_{21}+δa_{21}) ; a_{33} = a_{44} = (a_{31}+δa_{31}) ; a_{34} = (a_{32}+δa_{32}) ; a_{43} = (a_{41}+δa_{41}) 
D = det(A) 
D =
(p^{2}a_{21}a_{12})(a_{31}^{2}a_{41}a_{32})
 (2pa_{21})(2a_{31}a_{32})
+ (p^{2}a_{21}a_{12})(a_{32}a_{41}+a_{31}^{2})
+ (a_{12}a_{21}+p^{2})(a_{31}^{2}a_{41}a_{32})
 (2pa_{12})(2a_{31}a_{41})
+ (p^{2}a_{21}a_{12})(a_{31}^{2}a_{41}a_{32})
D = 2(p^{2}a_{21}a_{12})(a_{31}^{2}a_{41}a_{32}) 2(p^{2}+a_{21}a_{12})(a_{32}a_{41}+a_{31}^{2}) +4pa_{31}(a_{21}a_{32}+a_{12}a_{41}) = 4p^{2}a_{41}a_{32} 4a_{21}a_{12}a_{31}^{2} +4pa_{31}(a_{21}a_{32}+a_{12}a_{41}) D = 4cosθ/V_{P}cosη/V_{S} [p^{2}(2ρV_{S}^{2}p)^{2} +(ρ(12p^{2}V_{S}^{2}))^{2} +4ρ^{2}p^{2}V_{S}^{2}(12p^{2}V_{S}^{2})] = 4ρ^{2}cosθ/V_{P}cosη/V_{S} 
R_{φ} 
DR_{φ} =
(p^{2}a_{21}a_{12})((a_{31}+δa_{31})^{2}(a_{41}+δa_{41})(a_{32}+δa_{32}))
 (p(a_{21}+δa_{21})+pa_{21})(a_{32}(a_{31}+δa_{31})a_{31}(a_{32}+δa_{32}))
+ (p^{2}a_{21}(a_{12}+δa_{12}))(a_{32}(a_{41}+δa_{41})+a_{31}(a_{31}+δa_{31})) +
(a_{12}(a_{21}+δa_{21})+p^{2})(a_{31}(a_{31}+δa_{31})a_{41}(a_{32}+δa_{32}))
 (pa_{12}p(a_{12}+δa_{12}))(a_{31}(a_{41}+δa_{41})+a_{41}(a_{31}+δa_{31}))
+ (p^{2}(a_{21}+δa_{21})(a_{12}+δa_{12}))(a_{31}^{2}a_{41}a_{32}))
DR_{φ} = (p^{2}a_{21}a_{12})(2a_{31}δa_{31}a_{41}δa_{32}a_{32}δa_{41})  2pa_{32}a_{31}δa_{21} + (p^{2}a_{21}a_{12})(a_{32}δa_{41}+a_{31}δa_{31})a_{21}δa_{12}(a_{32}a_{41}+a_{31}^{2}) + (p^{2}+a_{12}a_{21})(a_{31}δa_{31}a_{41}δa_{32})+a_{12}δa_{21}(a_{31}^{2}a_{41}a_{32}) + 2pa_{12}(a_{31}δa_{41}+a_{41}δa_{31}) + (a_{21}δa_{12}+a_{12}δa_{21})(a_{31}^{2}+a_{41}a_{32}) DR_{φ} = 2(p(pa_{32}a_{12}a_{31})δa_{41} + a_{12}(pa_{41}a_{21}a_{31})δa_{31} + a_{31}(a_{31}a_{12}pa_{32})δa_{21})
DR_{φ} =
2ρ(2p^{2}cosη/V_{S}δ(ρV_{S}^{2}cosθ/V_{P})
+cosη/V_{S}cosθ/V_{P}δ(ρ(12p^{2}V_{S}^{2})
+ ρ(12p^{2}V_{S}^{2})cosη/V_{S}δ(cosθ/V_{P}))
R_{φ} =
(1/2)((14p^{2}V_{S}^{2})δρ/ρ
 8p^{2}V_{S}^{2}δV_{S}/V_{S}
 δ(cosθ/V_{P})/(cosθ/V_{P}))

R_{ψ} 
DR_{ψ} =
(2pa_{21})((a_{31}+δa_{31})^{2}(a_{41}+δa_{41})(a_{32}+δa_{32}))
 p((a_{21}+δa_{21})+a_{21})(a_{31}(a_{31}+δa_{31})a_{41}(a_{32}+δa_{32}))
+ (p^{2}a_{21}(a_{12}+δa_{12}))(a_{31}(a_{41}+δa_{41})+a_{41}(a_{31}+δa_{31})) +
p((a_{21}+δa_{21})+a_{21})(a_{31}(a_{31}+δa_{31})a_{41}(a_{32}+δa_{32}))
 (p^{2}a_{21}(a_{12}+δa_{12}))(a_{31}(a_{41}+δa_{41})+a_{41}(a_{31}+δa_{31}))
+ (p^{2}(a_{21}+δa_{21})(a_{12}+δa_{12}))(2a_{31}a_{41})
DR_{ψ} = 2pa_{21}(2a_{31}δa_{31}a_{41}δa_{32}a_{32}δa_{41})  2pa_{21}(a_{31}δa_{31}a_{41}δa_{32})pδa_{21}(a_{31}^{2}a_{41}a_{32}) + (p^{2}a_{21}a_{12})(a_{31}δa_{41})+a_{41}δa_{31})) + (pδa_{21})(a_{31}^{2}a_{41}a_{32})  (p^{2}a_{21}a_{12})(a_{31}δa_{41}+a_{41}δa_{31}) +a_{21}δa_{12}2a_{31}a_{41} + (2a_{31}a_{41}(a_{21}δa_{12}+a_{12}δa_{21})) DR_{ψ} = 2(p(a_{21}a_{31}pa_{41})δa_{31} + a_{21}(pa_{32}+a_{12}a_{31})δa_{41} + a_{41}(pa_{32}a_{31}a_{12})δa_{21})
DR_{ψ} =
2ρpcosθ/V_{P}(δ(ρ(12p^{2}V_{S}^{2}))
+ cosη/V_{S}δ(2ρV_{S}^{2}cosθ/V_{P})
+ 2ρV_{S}^{2}cosη/V_{S}δ(cosθ/V_{P}))
R_{ψ} =
(1/2)pV_{S}/cosη
((12p^{2}V_{S}^{2}+2V_{S}^{2}cosθ/V_{P}cosη/V_{S})δρ/ρ
+ 4V_{S}^{2}(p^{2}+cosη/V_{S}cosθ/V_{P})δV_{S}/V_{S})

In a similar way, an interface (z = 0) between two elastic halfspaces can guide a Stoneley wave (only for some values of velocities and densities). It is made of two coupled Rayleigh waves (two P_{1}, SV_{1} evanescent waves for z ≤ 0 coupled to two P_{2} , SV_{2} evanescent waves for z ≥ 0) propagating at the same horizontal velocity. V_{ST} is determined by looking for roots of the determinant of the matrix A obtained in the computation of the reflectiontransmission coefficients for PSV waves. There is not always a solution , except in the liquidsolid case.
Stoneley equation liquidsolid interf.  D_{L} = ρ_{2}^{2}cosθ_{1}/V_{P1}((12p^{2}V_{S2}^{2})^{2}+4V_{S2}^{4}p^{2}cosθ_{2}/V_{P2}cosη_{2}/V_{S2})
ρ_{1}ρ_{2}cosθ_{2}/V_{P2} = 0
(12p^{2}V_{S2}^{2})^{2}4V_{S2}^{4}p^{2}(p^{2}V_{P2}^{2}1)^{½}/V_{P2}(p^{2}V_{S2}^{2}1)^{½}/V_{S2} + (ρ_{1}/ρ_{2})((p^{2}V_{P2}^{2}1)^{½}/V_{P2})/((p^{2}V_{P1}^{2}1)^{½}/V_{P1}) = 0 
Stephen, R.A., F. CardoCasas, and C. H. Cheng , 1985, Finitedifference synthetic acoustic logs show the propagation of a Stoneley wave in a well for sonic logging.
standing wave in z
k_{z} for normal modes 
travelling wave in x , k_{x} = (ω^{2}/V_{S}^{2}k_{z}^{2})^{½}
phase velocity V(ω) = V_{ax} = ω/k_{x} 
ω ≥ highfrequency cutoff ω_{c}
k_{x} real for travelling wave in x  
normal modes 1D
C, C, G, C, E, G... 
e^{ikz}+e^{ikz} = 2cos(kz)
Hk_{n} = nπ , nλ_{n}/2 = H , n = 1,2,3...  
slab with free surfaces
R_{sup} = R_{inf} = 1 
e^{ikzz}+e^{ikzz} = 2cos(k_{z}z)
Hk_{z n} = nπ , nλ_{z n}/2 = H , n = 1,2,3... 
k_{x n} = (ω^{2}/V_{S}^{2}(nπ/H)^{2})^{½}
V_{n}(ω) = ω/k_{x n}  ω_{c n} = nπV_{S}/H 
layer (H ≤ z ≤ 0)
over halfspace (z≥ 0) R_{surf} = 1 ; R_{interf} = e^{iχ(p)} 
e^{ik1zz}+e^{iχ(p)}e^{ik1zz} = 2e^{iχ(p)/2}cos(k_{1z}zχ(p)/2)
Hk_{1z n}+χ(p)/2 = nπ , n = 0,1,2,3... 
k_{x n} = (ω^{2}/V_{S1}^{2}(nπχ(p)/2)^{2}/H^{2})^{½}
V_{L n}(ω) = ω/k_{x n} = 1/p_{n} 
k_{x nc} = ω_{c n}/V_{S2} ; χ(1/V_{S2}) = 0
ω_{c n} = nπ/(H(1/V_{S1}^{2}1/V_{S2}^{2})^{½}) 
displov.m 
k_{1z} = (ω/V_{S1})(1p^{2}V_{S1}^{2})^{½}
χ(p) = 2Arctan(ρ_{2}V_{S2}(p^{2}V_{S2}^{2}1)^{½}/ρ_{1}V_{S1}(1p^{2}V_{S1}^{2})^{½}) = 2Arctan(I_{2}/I_{1})
dispersion relation of Love wave :

An SV wave totally reflected in a layer such that V_{S1}<V_{S2}<V_{P1}<V_{P2} produces a guided Rayleigh wave. It propagates at the velocity V_{R} that is the apparent horizontal velocity of the SV waves in the layer, such that V_{S1}≤V_{R}≤V_{S2}. The total reflection on the free surface produces a phaseshift χ'_{1}(p) computed for PSV reflection coefficients on a free surface. The total reflection on the interface with the Halfspace produces a phaseshift χ'_{2}(p) computed for PSV reflection coefficients at the interface. The dispersion relation is written by taking into account the sum χ'_{12}(p) of these two phaseshifts for the standing wave condition in the layer. One obtains in this way the dispersion relation for the harmonics (n≥1 , V_{S1}≤V_{Rn}≤V_{S2}, V_{Rn}=V_{S2} at the cutoff frequency f_{cn}) of the Rayleigh wave guided in the layer.
boundary conditions  normal modes in z  travelling wave in x  cutoff frequency 
R_{surf} = e^{iχ'1(p)}
R_{interf} = e^{iχ'2(p)} χ'_{12}(p) = χ'_{1}(p)+χ'_{2}(p)  (nχ'_{12}(p)/2π)λ_{1z n}/2 = H
Hk_{1z n}+χ'_{12}(p)/2 = nπ  k_{x n} = (ω^{2}/V_{S1}^{2}(nπχ'_{12}(p)/2)^{2}/H^{2})^{½}
V_{R n}(ω) = ω/k_{x n} = 1/p_{n}  ω_{c n} = (nπχ'_{12}(1/V_{S2}))/2)/(H(1/V_{S1}^{2}1/V_{S2}^{2})^{½}) 
dispersion rel. of Rayleigh wave overtones n≥1 V_{S1}≤V_{R n}≤V_{S2} 
χ'_{1}(p) = 2Arctan((4V_{S1}^{4}p^{2}(p^{2}V_{P1}^{2}1)^{½}/V_{P1}(1p^{2}V_{S1}^{2})^{½}/V_{S1})/(12p^{2}V_{S1}^{2})^{2}) =
2Arctan(S)
χ'_{2}(p) = 2Arctan(Im(D)/Re(D)) = 2Arctan(I) tan(Hω(1p^{2}V_{S1}^{2})^{½}/V_{S1}) = tan((χ'_{1}(p)+χ'_{2}(p))/2) = (I+S)/(1IS) ω(p_{n}) = ((χ'_{1}(p_{n})+χ'_{2}(p_{n}))/2+nπ)V_{S1}/(H(1p_{n}^{2}V_{S1}^{2})^{½}) = (Arctan((I+S)/(1IS))+nπ)V_{S1}/(H(1p_{n}^{2}V_{S1}^{2})^{½}) 
The dispersion relation of the fundamental mode of the Rayleigh wave is obtained by considering four upgoings and downgoings P and SV waves in the layer and two evanescent P and SV waves in the halfspace. The propagation velocity V_{R0} correspond to the root of the determinant of the matrix expressing the boundary conditions on the free surface in z = H and on the interface on z = 0. It is such that V_{RS1}≤V_{R0}≤V_{RS2} where V_{RS1}is the Rayleigh wave velocity guided at the free surface of a halfspace of velocity V_{S1} and V_{RS2} that for a halfspace of velocity V_{S2}. At low frequency, V_{R0} = V_{RS2} (the wave with wavelength >> H does not see the layer).
boundary cond.  downgoing P1 , k_{zP1} = ω(1p^{2}V_{P1}^{2})^{½}/V_{P1}
Z_{P1} = e^{ikzP1H}  downgoing S1 , k_{zS1} = ω(1p^{2}V_{S1}^{2})^{½}/V_{S1}
Z_{S1} = e^{ikzS1H}  upgoing P1 , k_{zP1}  upgoing S1 , k_{zS1}  evanescent P2 , k_{zP2} = iω(p^{2}V_{P2}^{2}1)^{½}/V_{P2}  evanescent S2 , k_{zS2} = iω(p^{2}V_{S2}^{2}1)^{½}/V_{S2} 
z = H  σ_{zz}(k_{zP1}^{2})Z_{P1}  2μ_{1}k_{x}k_{zS1}Z_{S1}  σ_{zz}(k_{zP1}^{2})/Z_{P1}  2μ_{1}k_{x}k_{zS1}/Z_{S1}  0  0 
2μ_{1}k_{x}k_{zP1}Z_{P1}  σ_{zx}(k_{zS1}^{2})Z_{S1}  2μ_{1}k_{x}k_{zP1}/Z_{P1}  σ_{zx}(k_{zS1}^{2})/Z_{S1}  0  0  
z = 0  k_{x}  k_{zS1}  k_{x}  k_{zS1}  k_{x}  k_{zS2} 
k_{zP1}  k_{x}  k_{zP1}  k_{x}  k_{zP2}  k_{x}  
σ_{zz}(k_{zP1}^{2})  2μ_{1}k_{x}k_{zS1}  σ_{zz}(k_{zP1}^{2})  2μ_{1}k_{x}k_{zS1}  σ_{zz}(k_{zP2}^{2})  2μ_{2}k_{x}k_{zS2}  
2μ_{1}k_{x}k_{zP1}  σ_{zx}(k_{zS1}^{2})  2μ_{1}k_{x}k_{zP1}  σ_{zx}(k_{zS1}^{2})  2μ_{2}k_{x}k_{zP2}  σ_{zx}(k_{zS2}^{2}) 
The highfrequency part of the dispersion relation of the fundamental mode (V_{RS1}≤V_{R0}≤V_{S1}) corresponds to the Rayleigh wave guided by the free surface of the layer coupled to a Stoneley wave guided by the interface (or two P and SV evanescent waves in the layer from the surface coupled to four evanescent waves from the interface). At highfrequency, V_{R0} = V_{RS1} (the wave with wavelength << H sees only the layer).
boundary cond.  waves guided by free surface in z = H and prolongated in z = 0 by
Z_{P1} = e^{ω(p2VP121)½/VP1H} (0) ; Z_{S1} = e^{ω(p2VS121)½/VS1H} (0)  waves guided by interf. in z = 0 and prolongated in z = H by Z_{P1} , Z_{S1}  
z = H  Rayleigh free surf.  (12p^{2}V_{S1}^{2})Z_{P1} (0)  i2pV_{S1}(p^{2}V_{S1}^{2}1)^{½}Z_{S1} (0)  0  0  
+i2V_{S1}^{2}p(p^{2}V_{P1}^{2}1)^{½}/V_{P1}Z_{P1} (0)  +(12p^{2}V_{S1}^{2})Z_{S1} (0)  0  0  
z = 0  pZ_{P1} (0)  i(p^{2}V_{S1}^{2}1)^{½}/V_{S1}Z_{S1} (0)  Stoneley at interface  
i(p^{2}V_{P1}^{2}1)^{½}/V_{P1}Z_{P1} (0)  pZ_{S1} (0)  
ρ_{1}(12p^{2}V_{S1}^{2})Z_{P1} (0)  i2ρ_{1}V_{S1}^{2}p(p^{2}V_{S1}^{2}1)^{½}/V_{S1}Z_{S1} (0)  
i2ρ_{1}V_{S1}^{2}p(p^{2}V_{P1}^{2}1)^{½}/V_{P1}Z_{P1} (0)  ρ_{1}(12p^{2}V_{S1}^{2})Z_{S1} (0) 
Russel D., Waveguide
superposition of harmonic waves of phase velocities V(ω)=ω/k  constructive interference in (x,t) ⇔ stationnary phase for harmonic waves (the U variation is negligible) ⇒ group velocity V_{g} 
u(x,t) = ∫U(ω)e^{iω(x/V(ω)t)}dω = ∫U(ω)e^{i(k(ω)xωt)}dω = ∫U(k)e^{i(kxω(k)t)}dk  d(k(ω)xωt)/dω = (dk/dω)xt = 0 ⇒ V_{g} = dω/dk = V+kdV/dk = 1/(dk/dω) = 1/(p+ωdp/dω) 

displov.m 
dispersion relation of Love wave : tan(Hω(p_{n})(1p_{n}^{2}V_{S1}^{2})^{½}/V_{S1}) = (ρ_{2}V_{S2}(p^{2}V_{S2}^{2}1)^{½})/(ρ_{1}V_{S1}(1p^{2}V_{S1}^{2})^{½})
differentiation of dispersion relation :
dp_{n}/dω = (H/V_{S1})(1p_{n}^{2}V_{S1}^{2})/

Kelly K.R. , 1983, Numerical study of Love wave propagation
Roth M. and K. Holliger, 1999, Inversion of sourcegenerated noise in highresolution seismic data
z
I_{i} = ρ_{i}V_{Si}^{2}k_{iz} 
downgoing wave e^{kizz}
k_{iz} = ω(1/V_{Si}^{2}p^{2})^{½}  upgoing wave e^{ikizz} 
continuity of u_{y} and σ_{yz} at interface
total reflection if free surface (upgoing = downgoing) : dispersion relation 
(h_{3}+h_{2}+h_{1})
layer 3  D_{3}e^{ik3zh3} = D_{3}(a_{3}ib_{3})  M_{3}e^{ik3zh3} = M_{3}(a_{3}+ib_{3}) 
A_{3}b_{3}+a_{3}B_{3} = 0
tan(k_{3z}h_{3})+(I_{2}/I_{3})tan(k_{2z}h_{2})+((I_{1}/I_{3})(I_{1}/I_{2})tan(k_{3z}h_{3})tan(k_{2z}h_{2}))tan(k_{1z}h_{1}+χ/2) = 0 
(h_{2}+h_{1})
layer 3  D_{3} = A_{3}iB_{3}  M_{3} = A_{3}+iB_{3} 
A_{3} = A_{2}a_{2}B_{2}b_{2} = a_{1}a_{2}(I_{1}/I_{2})b_{1}b_{2}
B_{3} = (I_{2}/I_{3})(A_{2}b_{2}+B_{2}a_{2}) = (I_{2}/I_{3})a_{1}b_{2}+(I_{1}/I_{3})b_{1}a_{2} 
(h_{2}+h_{1})
layer 2  D_{2}e^{ik2zh2} = D_{2}(a_{2}ib_{2})  M_{2}e^{ik2zh2} = M_{2}(a_{2}+ib_{2}) 
A_{2}b_{2}+a_{2}B_{2} = 0
tan(k_{2z}h_{2})+(I_{1}/I_{2})tan(k_{1z}h_{1}+χ/2) = 0 
h_{1}
layer 2  D_{2} = A_{2}iB_{2}  M_{2} = A_{2}+iB_{2} 
A_{2} = a_{1}
B_{2} = (I_{1}/I_{2})b_{1} 
h_{1}
layer 1  D_{1}e^{ik1zh1} = a_{1}ib_{1}  M_{1}e^{ik1zh1} = a_{1}+ib_{1} 
b_{1} = 0
sin(k_{1z}h_{1}+χ/2) = 0 
0
layer 1  D_{1} = e^{iχ/2}  M_{1} = e^{iχ/2}  χ(p) = 2Arctan(ρV_{S}(p^{2}V_{S}^{2}1)^{½}/ρ_{1}V_{S1}(1p^{2}V_{S1}^{2})^{½}) = 2Arctan(I/I_{1}) 
> 0
1/2 space 
evanescent wave e^{kzz}
k_{z} = ω(p^{2}1/V^{2})^{½} 
Bessel equation  solution  derivative 
r^{2}d^{2}f/dr^{2}+rdf/dr+ (k_{r}^{2}r^{2}n^{2})f(r) = 0 
f(r) = H_{n}(k_{r}r) = J_{n}(k_{r}r)+iY_{n}(k_{r}r)
≈ (2/πk_{r}r)^{½}e^{i(krrπ/4nπ/2)}
H_{n}(ik_{r}r) ≈ i^{(n+1)}(2/πk_{r}r)^{½}e^{krr}  dH_{n}(k_{r}r)/dr = k_{r}H_{n1}(k_{r}r)(n/r)H_{n}(k_{r}r) = (n/r)H_{n}(k_{r}r)k_{r}H_{n+1}(k_{r}r) 
In the absence of azimutal variation, the displacement potentials for P et S waves are :
Δφ(r,z) = ∂^{2}φ/∂r^{2}+(1/r)∂φ/∂r+∂^{2}φ/∂z^{2} = (1/V_{P}^{2})∂^{2}φ/∂t^{2}  φ(r,z) = f(r)e^{i(kzzωt)}  k_{p}^{2} = ω^{2}/V_{P}^{2}k_{z}^{2}  f(r) = AH_{0}(k_{p}r) 
Δψ_{θ}(r,z) = (1/V_{S}^{2})∂^{2}ψ_{θ}/∂t^{2}  ψ_{θ}(r,z) = g(r)e^{i(kzzωt)}  k_{s}^{2} = ω^{2}/V_{S}^{2}k_{z}^{2}  g(r) = BH_{1}(k_{s}r) 
If there is an angular variation with θ, one has :
Δφ(r,θ,z) = Δφ(r,z) + (1/r)^{2}∂^{2}φ/∂θ^{2} = (1/V_{P}^{2})∂^{2}φ/∂t^{2}  φ(r,θ,z) =f(r)cos(nθ)e^{i(kzzωt)}  k_{p}^{2} = ω^{2}/V_{P}^{2}k_{z}^{2}  f(r) = AH_{n}(k_{p}r) 
Δψ_{r}(r,θ,z) = (1/V_{S}^{2})∂^{2}ψ_{r}/∂t^{2}  ψ_{r}(r,z) = g(r)sin(nθ)e^{i(kzzωt)}  k_{s}^{2} = ω^{2}/V_{S}^{2}k_{z}^{2}  g(r) = BH_{n+1}(k_{s}r) 
Δψ_{θ}(r,θ,z) = (1/V_{S}^{2})∂^{2}ψ_{θ}/∂t^{2}  ψ_{θ}(r,z) = g(r)cos(nθ)e^{i(kzzωt)}  
Δψ_{z}(r,θ,z) = (1/V_{S}^{2})∂^{2}ψ_{z}/∂t^{2}  ψ_{z}(r,z) = h(r)sin(nθ)e^{i(kzzωt)}  h(r) = CH_{n}(k_{s}r) 
Displacements and strains on the wall of a cylindrical borehole
u_{r} = cos(nθ)  u_{θ} = sin(nθ)  u_{z} = cos(nθ)  ε_{rr} = cos(nθ)  ε_{rθ} = sin(nθ)  ε_{rz} = cos(nθ)  
Ae^{i(kzzωt)}  n/rH_{n}(k_{p}r)k_{p}H_{n+1}(k_{p}r)  n/rH_{n}(k_{p}r)  ik_{z}H_{n}(k_{p}r)  (n(n1)/r^{2}k_{p}^{2})H_{n}(k_{p}r)+k_{p}/rH_{n+1}(k_{p}r)  n(n0.5)/r^{2}H_{n}(k_{p}r)+nk_{p}/rH_{n+1}(k_{p}r)  ik_{z}(n/rH_{n}(k_{p}r)k_{p}H_{n+1}(k_{p}r)) 
Be^{i(kzzωt)}  ik_{z}H_{n+1}(k_{s}r)  ik_{z}H_{n+1}(k_{s}r)  k_{s}H_{n}(k_{s}r)  ik_{z}(k_{s}H_{n}(k_{s}r)(n+1)/rH_{n+1}(k_{s}r))  ik_{z}(k_{s}/2H_{n}(k_{s}r)(n+1)/rH_{n+1}(k_{s}r))  (nk_{s}/rH_{n}(k_{s}r)+(k_{s}^{2}k_{z}^{2})H_{n+1}(k_{s}r))/2 
Ce^{i(kzzωt)}  n/rH_{n}(k_{s}r)  n/rH_{n}(k_{s}r)+k_{s}H_{n+1}(k_{s}r)  0  n(n1)/r^{2}H_{n}(k_{s}r)nk_{s}/rH_{n+1}(k_{s}r)  (n(n1)/r^{2}+k_{s}^{2}/2)H_{n}(k_{s}r)k_{s}/rH_{n+1}(k_{s}r)  (ik_{z}n/2r)H_{n}(k_{s}r) 
The P and S solutions in cylindrical coordinates are developped in : Sinha, B.K., E. Simsek, and S. Asvadurov, 2009 Influence of a pipe tool on borehole modes
Without body forces, the 3 equilibrium equations expressed in terms of the components of the particle displacement are :
ρ∂^{2}u/∂t^{2} = (λ+μ)graddivu+μΔu where
Δu = (Δu_{x}, Δu_{y}, Δu_{z}) = grad(divu)curl(curl(u))
The P wave equation, ∂^{2}φ/∂t^{2} = V_{P}^{2}Δφ,
is obtained from these 3 equations by writing u = gradφ.
Indeed, Δgradφ = gradΔφ as curl(gradφ) = 0.
Therefore, ρgrad∂^{2}φ/∂t^{2} = (λ+2μ)grad(Δφ).
The S wave equation, ∂^{2}ψ/∂t^{2} = V_{S}^{2}Δψ,
is obtained from these 3 equations by writing u = curlψ with divψ = 0.
Indeed, Δcurlψ = curlcurlcurlψ = curlΔψ as divcurlψ = 0 and divψ = 0.
Therefore, ρcurl∂^{2}ψ/∂t^{2} = μcurlΔψ.
For a source due to a point force with a module F(t) (N) and a direction e, the force per unit volume (N/m^{3}) is :
f = F(t)eδ(x) = F(t)/(4π)eΔ(1/r) = F(t)/(4π)Δ(e/r) =
(F(t)/(4π))(graddiv(e/r)curlcurl(e/r)).
Indeed, ∫δ(x)dx = 1 and
∫Δ(1/r)dx = ∫grad(1/r).dS =
(1/r^{2})e_{r}.4πr^{2}e_{r} = 4π
imply δ(x) = 1/(4π)Δ(1/r).
As e is constant, eΔ(1/r) = Δ(e/r).
In the harmonic case where the module of the force varies as a sinusoid of pulsation ω : F(t) = Fe^{iωt}.
θ being the angle between the direction e of the force and the radial direction e_{r} of propagation of the wave
from the point where the force applies, e = cosθe_{r}sinθe_{θ}.
One needs the expression of the differential operators
in spherical coordinates (r,&theta,ξ) :
gradφ(r,θ) = &partφ/∂re_{r}+1/r&partφ/∂θe_{θ} and
curl(ψ(r,θ)e_{ξ}) =
(1/r)(1/sinθ∂(sinθψ)/∂θe_{r}∂(rψ)/∂re_{θ}).
The equilibrium equations with the part of the source term producing P waves are :
ρgrad∂^{2}φ/∂t^{2} = (λ+2μ)grad(Δφ)(F(t)/(4π)(graddiv(e/r)
The P wave equation becomes :
∂^{2}φ/∂t^{2} = V_{P}^{2}Δφ(F(t)/4πρ)div(e/r).
It includes a source term proportional to
(dive/r) = grad(1/r).e = 1/r^{2}e_{r}.e = cosθ/r^{2}
The amplitude of the displacement potential and of the displacement in the far field is thus proprtional to cosθ.
By writing the particle displacement potential as the divergence of a vector Φ(r)e, one separates the dependency in r et θ.
u = gradφ(r,θ)  ∂^{2}φ/∂t^{2} = V_{P}^{2}Δφ(F(t)/4πρ)div(e/r)  φ(r,θ) = div(Φ(r)e)  ∂^{2}Φ/∂t^{2} = V_{P}^{2}ΔΦF(t)/(4πρr)  ∂^{2}(rΦ)/∂t^{2} = V_{P}^{2}∂^{2}(rΦ)/∂r^{2}F(t)/(4πρ) 
F(t) = Fe^{iωt}, k = ω/V_{P}  ∂^{2}(rΦ)/∂r^{2}+k^{2}(rΦ) = F/(4πρV_{P}^{2})  Φ = F/(4πρV_{P}^{2}k^{2})(1e^{ikr})/r  φ = gradΦ(r).e = cosθ∂Φ/∂r = Fcosθ/(4πρV_{P}^{2}k^{2})((ik1/r)e^{ikr}/r+1/r^{2})  
u =
cosθ∂^{2}Φ/∂r^{2}e_{r}1/r∂Φ/∂rsinθe_{θ} =
cosθ∂^{2}Φ/∂r^{2}e_{r}+1/r∂Φ/∂r(ecosθe_{r})
u = ( F/(4πρV_{P}^{2}k^{2})) ((cosθ(2/r^{2}2ik/rk^{2})e_{r}+1/r(ik1/r)(ecosθe_{r}))e^{ikr}/r) (1/r^{3})(3cosθe_{r}e))  u_{cl} = Fcosθe_{r}/(4πρV_{P}^{2})e^{ikr}/r  u_{cp} = F/(4πρV_{P}^{2})(e3cosθe_{r})((i/kr1/(kr)^{2})e^{ikr}/r+1/r^{3}) 
The equilibrium equations with the part of the source term producing S waves are :
ρcurl∂^{2}ψ/∂t^{2} = μcurlΔψ+(F(t)/(4π))curlcurl(e/r).
The S wave equation becomes :
∂^{2}ψ/∂t^{2} =V_{S}^{2}Δψ+(F(t)/(4πρ))curl(e/r).
It includes a source term proportional to
curl(e/r) = grad(1/r)∧e = (1/r^{2})e_{r}∧e =
sinθe_{ξ}/r^{2}
The amplitude of the displacement potential and of the displacement in the far field is thus proportional to sinθ.
By writing the particle displacement potential as the curl of a vector Ψ(r)e, one separates the dependency in r et θ.
u = curlψ(r,θ)  ∂^{2}ψ/∂t^{2} = V_{S}^{2}Δψ+(F(t)/4πρ)curl(e/r)  ψ(r,θ) = curl(Ψ(r)e)  ∂^{2}Ψ/∂t^{2} = V_{S}^{2}ΔΨF(t)/(4πρr)  ∂^{2}(rΨ)/∂t^{2} = V_{S}^{2}∂^{2}(rΨ)/∂r^{2}F(t)/(4πρ) 
F(t) = Fe^{iωt}, k = ω/V_{S}  ∂^{2}(rΨ)/∂r^{2}+k^{2}(rΨ) = F/(4πρV_{S}^{2})  Ψ = F/(4πρV_{S}^{2}k^{2})(1e^{ikr})/r  ψ(r,θ) = gradΨ(r)∧e = sinθ∂Ψ/∂re_{ξ} = Fsinθe_{ξ}/(4πρV_{S}^{2}k^{2})((ik1/r)e^{ikr}/r+1/r^{2})  
u = (1/r)(2cosθe_{r}∂Ψ/∂rsinθe_{θ}∂(r∂Ψ/∂r)/∂r) =
(1/r)(2(e+sinθe_{θ})∂Ψ/∂rsinθe_{θ}(∂Ψ/∂r+r∂^{2}Ψ/∂r^{2}))
u = F/(4πρV_{S}^{2}k^{2}) (1/r)((2(e+sinθe_{θ})(ik1/r)sinθe_{θ}(rk^{2}ik+1/r))e^{ikr}/r +(1/r^{2})(2e+3sinθe_{θ}))  u_{cl} = Fsinθe_{θ}/(4πρV_{S}^{2})e^{ikr}/r  u_{cp} = F/(4πρV_{S}^{2})(e3cosθe_{r}) ((i/kr1/(kr)^{2})e^{ikr}/r+1/r^{3}) 
The contribution of the terms of displacement in the near and far field to the displacement and polarization of P and S waves at distances of the source of the order
of the wavelength are illustrated in these figures :



force.m
isotropic case : 2 elastic coef. λ and μ  transverse isotropy (symmetry axis Oz) : 5 elastic coef.  example for clays (x10^{9}Pa) 
σ_{xx} = c_{11}ε_{xx}+c_{12}ε_{yy}+c_{12}ε_{zz}
σ_{xy} = (c_{11}c_{12})ε_{xy} c_{11} = λ+2μ , c_{12} = λ 
σ_{xx} = c_{11}ε_{xx}+c_{12}ε_{yy}+c_{13}ε_{zz} ,
σ_{yy} = c_{12}ε_{xx}+c_{11}ε_{yy}+c_{13}ε_{zz}
σ_{zz} = c_{13}ε_{xx}+c_{13}ε_{yy}+c_{33}ε_{zz} σ_{yz} = 2c_{44}ε_{yz} , σ_{zx} = 2c_{44}ε_{zx} σ_{xy} = 2c_{66}ε_{xy} , c_{66} = (c_{11}c_{12})/2 
c_{11} = 45 , c_{12} = 10 , c_{13} = 8
c_{33} = 28 c_{44} = 11 c_{66} = (c_{11}c_{12})/2 = 17.5 
σ , ε applied on thinly stratified medium of thickness H  horizontal layers of thicknesses h_{i}, p_{i} = h_{i}/H  anisotropic equivalent moduli 
layers compressed in series vertically : σ_{zz} ≠ 0
without lateral strain : ε_{xx} = ε_{yy} = 0 
σ_{zzi} = σ_{zz} , ε_{zzi} = σ_{zz}/(λ_{i}+2μ_{i})
ε_{zz} = ∑p_{i}ε_{zzi} = σ_{zz}∑p_{i}/(λ_{i}+2μ_{i}) σ_{xxi} = σ_{yyi} = λ_{i}ε_{zzi} 
σ_{zz} = c_{33}ε_{zz}
c_{33} = 1/∑p_{i}/(λ_{i}+2μ_{i}) 
layers compressed in parallel horizontaly along x : σ_{xx} ≠ 0
without strain along y : ε_{yy} = 0 without average vertical strain : ε_{zz} = ∑p_{i}ε_{zzi} = 0 uniform strain along x : ε_{xxi} = ε_{xx} 
σ_{xxi} = (λ_{i}+2μ_{i})ε_{xx}+λ_{i}ε_{zzi}
σ_{xx} = ∑p_{i}σ_{xxi} = ε_{xx}∑p_{i}(λ_{i}+2μ_{i})+∑p_{i}λ_{i}ε_{zzi} [1] σ_{zzi} = σ_{zz} , σ_{zz} = λ_{i}ε_{xx}+(λ_{i}+2μ_{i})ε_{zzi} [2] [2] ⇒ σ_{zz}∑p_{i}/(λ_{i}+2μ_{i}) = ε_{xx}∑p_{i}λ_{i}/(λ_{i}+2μ_{i}) + ∑p_{i}ε_{zzi} = ε_{xx}∑p_{i}λ_{i}/(λ_{i}+2μ_{i}) [3] [2,3] ⇒ ∑p_{i}λ_{i}ε_{zzi} = ε_{xx}((∑p_{i}λ_{i}/(λ_{i}+2μ_{i}))^{2}/∑p_{i}/(λ_{i}+2μ_{i})∑p_{i}λ_{i}^{2}/(λ_{i}+2μ_{i})) [4] 
σ_{xx} = c_{11}ε_{xx}
[1,4] ⇒ c_{11} = c_{33}(∑p_{i}λ_{i}/(λ_{i}+2μ_{i}))^{2}+4∑p_{i}μ_{i}(λ_{i}+μ_{i})/(λ_{i}+2μ_{i}) 
layers sheared in series horizontally : σ_{yz} ≠ 0 
σ_{yz} = 2μ_{i}ε_{yzi}
ε_{yz} = ∑p_{i}ε_{yzi} = (σ_{yz}/2)∑p_{i}/μ_{i} 
σ_{yz} = 2c_{44}ε_{yz}
c_{44} = 1/∑p_{i}/μ_{i} 
layers sheared in parallel horizontally : σ_{xy} ≠ 0
uniform strain : ε_{xyi} = ε_{xy} 
σ_{xyi} = 2μ_{i}ε_{xy}
σ_{xy} = ∑p_{i}σ_{xyi} = 2∑p_{i}μ_{i}ε_{xy} 
σ_{xy} = 2c_{66}ε_{xy}
c_{66} = ∑p_{i}μ_{i} , c_{12} = c_{11}2c_{66} 
layers compressed in parallel horizontally along x : σ_{xx} ≠ 0
without vertical stress : σ_{zz} = 0 without strain along y : ε_{yy} = 0 uniform strain along x : ε_{xxi} = ε_{xx} 
σ_{xxi} = (λ_{i}+2μ_{i})ε_{xx}+λ_{i}ε_{zzi}
σ_{zzi} = σ_{zz} = 0 , λ_{i}ε_{xx}+(λ_{i}+2μ_{i})ε_{zzi} = 0 ε_{zz} = ∑p_{i}ε_{zzi} = ε_{xx}∑p_{i}λ_{i}/(λ_{i}+2μ_{i}) 
σ_{zz} = c_{13}ε_{xx}+c_{33}ε_{zz} = 0
c_{13} = c_{33}(∑p_{i}λ_{i}/(λ_{i}+2μ_{i})) 
Here is a plot of these different strainstress states and the equivalent elastic moduli obtained for a stratification with an equal proportion of slow and fast layers :
anisomod.m
plane SH waves (horizontal polarization)  plane quasiP waves (quasilongitudinal polarization α : α ≈ θ)
plane quasiSV waves (quasitransverse polarization β : β ≈ ηπ/2) 
u_{y}(x,z,t) = U_{y}e^{i(kxx+kzzωt)}
ρω^{2} = (c_{66}k_{x}^{2}+c_{44}k_{z}^{2}) the dispersion is an ellipse with semiaxes (c_{66}/ρ)^{½} and (c_{44}/ρ)^{½}  u(x,z,t) = (U_{x},0,U_{z})e^{i(kxx+kzzωt)}
ρω^{2}U_{x} = (c_{11}k_{x}^{2}+c_{44}k_{z}^{2})U_{x}+(c_{13}+c_{44})k_{x}k_{z}U_{z} ρω^{2}U_{z} = (c_{13}+c_{44})k_{x}k_{z}U_{x}+(c_{44}k_{x}^{2}+c_{33}k_{z}^{2})U_{z} (ρω^{2}c_{11}k_{x}^{2}c_{44}k_{z}^{2})(ρω^{2}c_{44}k_{x}^{2}c_{33}k_{z}^{2})((c_{13}+c_{44})k_{x}k_{z})^{2} = 0 
V_{SH}^{2} = (1/ρ)(c_{66}sin^{2}η+c_{44}cos^{2}η)  V_{QP}^{2} = (1/2ρ)(c_{44}+c_{11}sin^{2}θ+c_{33}cos^{2}θ+D)
V_{QS}^{2} = (1/2ρ)(c_{44}+c_{11}sin^{2}θ+c_{33}cos^{2}θD) D^{2} = ((c_{11}c_{44})sin^{2}θ(c_{33}c_{44})cos^{2}θ)^{2}+(c_{13}+c_{44})^{2}sin^{2}2θ tanα = (c_{13}+c_{44})sin2θ/2(ρV_{QP}^{2}c_{11}sin^{2}θc_{44}cos^{2}θ) = (c_{13}+c_{44})sin2θ/((c_{44}c_{11})sin^{2}θ+(c_{33}c_{44})cos^{2}θ+D) tanβ = (c_{13}+c_{44})sin2θ/2(ρV_{QS}^{2}c_{11}sin^{2}θc_{44}cos^{2}θ) = (c_{13}+c_{44})sin2θ/((c_{44}c_{11})sin^{2}θ+(c_{33}c_{44})cos^{2}θD) tanα.tanβ = ((c_{13}+c_{44})sin2θ)^{2}/(((c_{44}c_{11})sin^{2}θ+(c_{33}c_{44})cos^{2}θ)^{2}D^{2}) = 1 ⇒ β = απ/2 
group velocity ⇔ stationnary phase in ∫∫U(k_{x},k_{z})e^{i(kxx+kzzω(kx,kz)t)}dk_{x}dk_{z}
in order to get constructive interferences of harmonic waves in (x,z,t)
∂(k_{x}x+k_{z}zω(k_{x},k_{z})t)/∂k_{x} = 0 , x(∂ω/∂k_{x})t = 0 , V_{Gx} = ∂ω/∂k_{x} ∂(k_{x}x+k_{z}zω(k_{x},k_{z})t)/∂k_{z} = 0 , z(∂ω/∂k_{z})t = 0 , V_{Gz} = ∂ω/∂k_{z} V_{G} = grad_{k}ω(k_{x},k_{z}) ⇒ e_{k}.V_{G} = e_{k}.grad_{k}ω(k_{x},k_{z}) = dω/dk = V_{phase} ⇒ V_{G} ≥ V_{phase}  
V_{Gx} = (c_{66}/ρω)k_{x} , V_{Gz} = (c_{44}/ρω)k_{z}
V_{G}^{2} = (1/ρ^{2}ω^{2})((c_{66}k_{x})^{2}+(c_{44}k_{z})^{2}) sinη_{g} = V_{Gx}/V_{G} = (c_{66}k_{x})/((c_{66}k_{x})^{2}+(c_{44}k_{z})^{2})^{½} cosη_{g} = V_{Gz}/V_{G} = (c_{44}k_{z})/((c_{66}k_{x})^{2}+(c_{44}k_{z})^{2})^{½} V_{G}^{2} = 1/(ρ(sin^{2}η_{g}/c_{66}+cos^{2}η_{g}/c_{44}))  V_{Gx} = ∂(kV)/∂k_{x} = V∂k/∂k_{x}+k∂V/∂k_{x}
k = (k_{x}^{2}+k_{z}^{2})^{½} , ∂k/∂k_{x} = k_{x}/k = sinθ k_{x} = ksinθ , 1 = ∂k/∂k_{x}sinθ+kcosθ∂θ/∂k_{x} , ∂θ/∂k_{x} = cosθ/k V_{Gx} = Vsinθ+cosθdV/dθ , V_{Gz} = VcosθsinθdV/dθ V_{G}^{2} = V^{2}+(dV/dθ)^{2} tanθ_{g} = V_{Gx}/V_{Gz} = (tanθ+(1/V)dV/dθ)/(1(tanθ/V)dV/dθ) 
Thomsen anisotropy parameters used in the processing of seismic reflection data : α_{T} , β_{T} , γ_{T} , δ_{T} , ε_{T}  
η = 0 , V_{SH}^{2} = c_{44}/ρ = β_{T}^{2} , η = 90° , V_{SH}^{2} = c_{66}/ρ
γ_{T} = (V_{SH}^{2}(90)V_{SH}^{2}(0))/2V_{SH}^{2}(0) = (c_{66}c_{44})/2c_{44} V_{SH}^{2} = β_{T}^{2}((c_{66}/c_{44})sin^{2}η+cos^{2}η) = β_{T}^{2}(1+2γ_{T}sin^{2}η) V_{G}^{2} = β_{T}^{2}((1+2γ_{T})/(1+2γ_{T}cos^{2}η_{g}))  θ = 0 , D = (c_{33}c_{44}) , V_{QP}^{2} = c_{33}/ρ = α_{T}^{2} , V_{QS}^{2} = c_{44}/ρ = β_{T}^{2}
θ = 90° , D = (c_{11}c_{44}) , V_{QP}^{2} = c_{11}/ρ , V_{QS}^{2} = c_{44}/ρ = β_{T}^{2} ε_{T} = (V_{QP}^{2}(90)V_{QP}^{2}(0))/2V_{QP}^{2}(0) = (c_{11}c_{33})/2c_{33} δ_{T} = (1/2α_{T})(d^{2}V_{QP}/dθ^{2})(0) = ((c_{13}+c_{44})^{2}(c_{33}c_{44})^{2})/(2c_{33}(c_{33}c_{44})) 
example for clays (ρ=2380 kg/m^{3}) : α_{T} = 3.43 km/s, β_{T} = 2.15 km/s , γ_{T} = 0.3 , δ_{T} = 0.07 , ε_{T} = 0.3  
weak anisotropy
γ_{T}<<1 , γ_{T} ≈ (V_{SH}(90)V_{SH}(0))/V_{SH}(0) V_{SH} ≈ β_{T}(1+γ_{T}sin^{2}η)  ε_{T}<<1 , δ_{T}<<1 , ε_{T} ≈ (V_{QP}(90)V_{QP}(0))/V_{QP}(0)
V_{QP} ≈ α_{T}(1+δ_{T}sin^{2}θcos^{2}θ+ε_{T}sin^{4}θ) = α_{T}(1+δ_{T}sin^{2}θ+(ε_{T}δ_{T})sin^{4}θ) V_{QS} ≈ β_{T}(1+(α_{T}^{2}/β_{T}^{2})(ε_{T}δ_{T})sin^{2}θcos^{2}θ) 
Here are the phase and group velocities and the polarizations obtained with the values of the above table, and also
groups of SH et QP waves for different ranges of phase angles :


anisopro.m
To compute the QPQSV reflectiontransmission coefficients, the components of particle motion are written as a function of the polarization angle α for the QP waves and of the angle β = QSV polarization angle + π/2 for the QS waves (same as in the isotropic case). Snell's law still applies but it does not give explicitely the different incidence angles since velocities depend on these angles. The dispersion relation provides the angles θ(p) and η(p), or more generally k_{z}^{P} = ωq^{P}(p) and k_{z}^{S} = ωq^{S}(p), wich are imaginary for evanescent waves.
Snell's law : p = sinθ_{1}/V_{QP1}(θ_{1}) = sinθ_{2}/V_{QP2}(θ_{2}) = sinη_{1}/V_{QS1}(η_{1}) = sinη_{2}/V_{QS2}(η_{2}) ;
q^{P1}(p) = (cosθ_{1}/V_{QP1})(p) ; q^{P2}(p) = (cosθ_{2}/V_{QP2})(p) ; q^{S1}(p) = (cosη_{1}/V_{QS1})(p) ; q^{S2}(p) = (cosη_{1}/V_{QS2})(p) ; Dispersion relation : (ρc_{11}p^{2}c_{44}q^{2})(ρc_{44}p^{2}c_{33}q^{2})((c_{13}+c_{44})pq)^{2} = 0 Aq^{4}+Bq^{2}+C = 0 ; A = c_{33}c_{44} ; B = c_{33}(ρc_{11}p^{2})c_{44}(ρc_{44}p^{2})(c_{13}+c_{44})^{2}p^{2} ; C = (ρc_{11}p^{2})(ρc_{44}p^{2}) ; Δ = (B^{2}4AC)^{½} q^{P} = ((BΔ)/2A)^{½} ; q^{S} = ((B+Δ)/2A)^{½} ; Polarization : α = Arctan((c_{13}+c_{44})pq^{P}/(ρc_{11}p^{2}c_{44}(q^{P})^{2})) ; β = π/2+Arctan((c_{13}+c_{44})pq^{S}/(ρc_{11}p^{2}c_{44}(q^{S})^{2}))  
QP incident downgoing + reflected upgoing U^{I} = (sinα_{1},cosα_{1}) , U^{RP} = R_{P}(sinα_{1},cosα_{1}) k_{z}^{I} = ωq^{P1}(p) , k_{z}^{RP} = ωq^{P1}(p)  QSV reflected upgoing U^{RS} = R_{S}(cosβ_{1},sinβ_{1}) k_{z}^{RS} = ωq^{S1}(p)  QP transmitted downgoing U^{TP} = T_{P}(sinα_{2},cosα_{2}) k_{z}^{TP} = ωq^{P2}(p)  QSV transmitted downgoing U^{TS} = T_{S}(cosβ_{2},sinβ_{2}) k_{z}^{TS} = ωq^{S2}(p)  
U_{x} =  sinα_{1}(1+R_{P})  + cosβ_{1}R_{S}  = sinα_{2}T_{P}  cosβ_{2}T_{S} 
U_{z} =  cosα_{1}(1R_{P})  + sinβ_{1}R_{S}  = cosα_{2}T_{P}  + sinβ_{2}T_{S} 
Σ_{zz} =  (c_{13}sinα_{1}p+c_{33}cosα_{1}q^{P1})(1+R_{P})  + (c_{13}cosβ_{1}pc_{33}sinβ_{1}q^{S1})R_{S}  = (d_{13}sinα_{2}p+d_{33}cosα_{2}q^{P2})T_{P}  + (d_{13}cosβ_{2}p+d_{33}sinβ_{2}q^{S2})T_{S} 
Σ_{zx} =  c_{44}(sinα_{1}q^{P1}+cosα_{1}p)(1R_{P})  + c_{44}(cosβ_{1}q^{S1}+sinβ_{1}p)R_{S}  = d_{44}(sinα_{2}q^{P2}+cosα_{2}p)T_{P}  + d_{44}(cosβ_{2}q^{S2}+sinβ_{2}p)T_{S} 
écriture matricielle : AR = B ; R = [R_{P} , R_{S} , T_{P} , T_{S}]^{t} ; B = [sinα_{1} , cosα_{1} , (c_{13}sinα_{1}p+c_{33}cosα_{1}q^{P1}) , c_{44}(sinα_{1}q^{P1}+cosα_{1}p)]^{t}  
A  sinα_{1}  cosβ_{1}  sinα_{2}  cosβ_{2} 
cosα_{1}  sinβ_{1}  cosα_{2}  sinβ_{2}  
(c_{13}sinα_{1}p+c_{33}cosα_{1}q^{P1})  (c_{13}cosβ_{1}pc_{33}sinβ_{1}q^{S1})   (d_{13}sinα_{2}p+d_{33}cosα_{2}q^{P2})  (d_{13}cosβ_{2}pd_{33}sinβ_{2}q^{S2})  
c_{44}(sinα_{1}q^{P1}+cosα_{1}p)  c_{44}(cosβ_{1}q^{S1}sinβ_{1}p)  d_{44}(sinα_{2}q^{P2}+cosα_{2}p)  d_{44}(cosβ_{2}q^{S2}+sinβ_{2}p) 
Here are the coefficients computed for an interface between a slow isotropic medium and a fast transverse isotropic medium in the case where a P wave is incident
in the isotropic medium :
anisopro.m
Winterstein D.F. and G. S. De, 2001, VTI documented
document the effect of transverse isotropy on S waves recorded in a 9component vertical seismic profile.
Johnston, J. E. and Christensen, N. I. 1995, Seismic anisotropy of shales
measure phase and group velocities in shales.
Godfrey, N.J. ; Christensen, N. I. and Okaya, D.A. 2000,
Anisotropy of schists: Contribution of crustal anisotropy to active source seismic experiments and shear wave splitting observations.
displacement
equilibrium eq.  cinetic energy density (Jm^{3})  elastic energy density (Jm^{3})  power density P_{t} (Wm^{2})  intensity I (Wm^{2}) 
u(x,t)
ρ∂^{2}u/∂t^{2} = ∂σ/∂x 
w_{c} = (ρ/2)(∂u/∂t)^{2}
∂w_{c}/∂t = (∂σ/∂x)(∂u/∂t) 
w_{e} = ∫σdε = (E/2)ε^{2}
∂w_{e}/∂t = σ∂/∂x(∂u/∂t) 
∂(w_{c}+w_{e})/∂t = ∂/∂x(σ∂u/∂t)
P_{t} = σ∂u/∂t  (1/T)∫_{0}^{T}P_{t}dt 
u(x,t) = Ucos(kxωt)  w_{c} = (ρ/2)ω^{2}U^{2}sin^{2}(kxωt)  w_{e} = (ρ/2)V^{2}k^{2}U^{2}sin^{2}(kxωt)  P_{t} = ρVω^{2}U^{2}sin^{2}(kxωt) = V(w_{c}+w_{e})  (ρ/2)Vω^{2}U^{2} 
u = gradφ(x,t)
ρ∂^{2}u/∂t^{2} = gradp 
w_{c} = (ρ/2)(∂u/∂t)^{2}
∂w_{c}/∂t = gradp.∂u/∂t 
w_{e} = (K/2)(divu)^{2}
∂w_{e}/∂t = pdiv∂u/∂t 
∫dx∂(w_{c}+w_{e})/∂t = ∫dxdiv(p∂u/∂t) = ∫(p∂u/∂t).dS
P_{t} = p∂u/∂t  (1/2)ℜ(p∂u*/∂t) 
2π/quality factor = energy dissipated per cycle / maximum elastic energy  cycle = λ  amplitude  displacement  complex velocity V_{Q} = V(1i/2Q) 
2π/Q = (δw)_{cycle}/w_{max}  π/Q = (dU/dx)λ/U  U(x) = Ue^{αx}
α = π/Qλ = k/2Q = ω/2QV  u(x,t) = Ue^{αx}cos(kxωt)  u(x,t) = Ue^{αx}e^{iω(x/Vt)} = Ue^{iω(x(1+i/2Q)/Vt)}
u(x,t) = Ue^{iω(x/VQt)} 
linear viscoelasticity  energy dissipated per cycle  maximum elastic energy  qualty factor 
σ = σ_{m}cos(ωt+χ)
ε = ε_{m}cosωt  (δw)_{cycle} = ∫_{0}^{T}σdε = πσ_{m}ε_{m}sinχ  w_{max} = ∫_{0}^{T/4}σdε ≈ (1/2)σ_{m}ε_{m}cosχ  1/Q = tanχ ≈ χ 
viscoelastic model , elastic modulus E , viscosity η  complex wavenumber k+iα  weak attenuation , 1/Q << 1  1/Q < 1 , dispersion V(ω) = ω/k 
E and η in parallel, relaxation time τ = η/E
σ = Eε+η∂ε/∂t = E(1iωτ)ε = Mε M = E(1+ω^{2}τ^{2})^{½}e^{iχ} , 1/Q = tanχ = ωτ  ρ∂^{2}u/∂t^{2} = E(1i/Q)∂^{2}u/∂x^{2}
ρω^{2} = E(1i/Q)(k+iα)^{2} , k_{0}^{2} = ρω^{2}/E k_{0}^{2} = (1i/Q)(k+iα)^{2}  k+iα = k_{0}(1+i/2Q)
k = k_{0} , V = V_{0} = ω/k_{0} α = k_{0}/(2Q) = ω^{2}τ/(2V_{0})  k^{2}+α^{2} = k_{0}^{2}/(1+1/Q^{2})^{½} ≈ k_{0}^{2}(11/2Q^{2})
k^{2}α^{2} = k_{0}^{2}/(1+1/Q^{2}) ≈ k_{0}^{2}(11/Q^{2}) V(ω) ≈ V_{0}(1+3/(8Q^{2})) = V_{0}(1+3ω^{2}τ^{2}/8) , α = k_{0}/(2Q) 
(E and η in parallel) in series with E' , σ = M''ε
1/M'' = 1/M+1/E' = 1/E(1iωτ)+1/E' = (1/E+1/E'iωτ/E')/(1iωτ) 1/E'' = 1/E+1/E' , τ'' = τE''/E' , M'' = E''(1iωτ)/(1iωτ'') = M''e^{i(χχ'')} 1/Q = tan(χχ'') = ω(ττ'')/(1+ω^{2}ττ'') d(1/Q)/dω = 0 pour ω_{r} = (ττ'')^{½} and (1/Q)_{max} = (ττ'')/2(ττ'')^{½}  ρω^{2} = M''(k+iα)^{2} , k''_{0}^{2} = ρω^{2}/E''
k^{2}+α^{2} = ρω^{2}/M'' = k''_{0}^{2}(1+ω^{2}τ''^{2})^{½}/(1+ω^{2}τ^{2})^{½} k^{2}α^{2} = ρω^{2}ℜ(1/M'') = k''_{0}^{2}(1+ω^{2}ττ'')/(1+ω^{2}τ^{2}) k^{2} = ω^{2}/V(ω)^{2} = (k''_{0}^{2}/2)((1+ω^{2}τ''^{2})^{½}/(1+ω^{2}τ^{2})^{½}+(1+ω^{2}ττ'')/(1+ω^{2}τ^{2})) 
constant Q model (Kjartansson, 1979, JGR, 84, 4737)
M = M_{0}(iω/ω_{0})^{2γ} = M_{0}(ω/ω_{0})^{2γ}e^{iπγsign(ω)} 1/Q = tan(πγ)  ρω^{2} = M_{0}(ω/ω_{0})^{2γ}e^{iπγsign(ω)}(k+iα)^{2}
k^{2}+α^{2} = k_{0}^{2}(ω/ω_{0})^{2γ} k^{2}α^{2} = k_{0}^{2}(ω/ω_{0})^{2γ}cos(πγ)  k = k_{0}(ω/ω_{0})^{γ}cos(πγ/2)
V(ω) = V_{0}(ω/ω_{0})^{γ}/cos(πγ/2) α = k_{0}(ω/ω_{0})^{γ}sin(πγ/2) = ktan(πγ/2)  dispersion pour 1/Q << 1
V(ω)/V_{0} ≈ (ω/ω_{0})^{1/πQ} = e^{Log(ω/ω0)/πQ} V(ω)/V_{0} ≈ 1+Log(ω/ω_{0})/πQ 
Φ = Ω_{Φ}/Ω = Ω_{f}/Ω  ε = δΩ/Ω  ζ = (δΩ_{Φ}δΩ_{f})/Ω = Φ(δΩ_{Φ}/Ω_{Φ}δΩ_{f}/Ω_{f}) = Φ(ε_{Φ}ε_{f}) 
La théorie de la poroélasticité linéaire isotrope exprime ε and ζ comme une combinaison linéaire de σ and p_{f} , la variation de pression de fluide dans les pores. La signification des coefficients de proportionnalité α (coefficient de BiotWillis) and B (coefficient de Skempton) apparaît en considérant les cas d'une roche drainée (pression de fluide constante, incompressibilité K), nondrainée (contenu en fluide constant, incompressibilité K_{sat}) and d'une pression différentielle nulle (pressions de confinement and de fluide égales). Dans ce dernier cas, la déformation homogène est celle de la roche sans pore d'incompressibilité K_{0} and la déformation volumique ε_{Φ} des pores est égale à la déformation volumique ε de la roche : on obtient ainsi la relation de Gassmann entre K_{sat}, K_{0} , K , K_{f} (l'incompressibilité du fluide) and Φ. Candte relation permet de déterminer la vitesse de propagation des waves de compression dans les roches poreuses saturées en fluide à basse frequency (f < 100Hz).
poroélasticité linéaire  ε = (1/K)(σ+αp_{f})
ζ = (α/K)(σ+p_{f}/B)  incompressibilité  coef. poroélastiques 
p_{f} = 0  ε = (1/K)σ  K  α = ζ/ε 
ζ = 0  ε = ((1αB)/K))σ  K = (1αB)K_{sat}  B = p_{f}/σ 
σ = p_{f}  ε = ((1α)/K))σ = σ/K_{0}  K = (1α)K_{0}  . 
ε_{Φ} = ζ/Φ+ε_{f} = ε
ε_{f} = p_{f}/K_{f}  (α/ΦK)(11/B)+1/K_{f} = 1/K_{0}  1/B = 1(ΦK/α)(1/K_{0}1/K_{f}) = α/(1K/K_{sat})  
relation de Gassmann 
1/K_{sat} = (1α)/K+α/K(11/(1(ΦK/α)(1/K_{0}1/K_{f})))
1/K_{sat} = 1/K_{0}+Φ/(ΦK/α+1/(1/K_{f}1/K_{0})) 1/(1/K_{sat}1/K_{0}) = K/α+1/(Φ(1/K_{f}1/K_{0})) = 1/(1/K1/K_{0})+1/(Φ(1/K_{f}1/K_{0}))  
densité and vitesses  ρ = (1Φ)ρ_{0}+Φρ_{f}  V_{P}^{2} = (K_{sat}+4μ/3)/ρ  V_{S}^{2} = μ/ρ 
Wang Z., 2001, Fundamentals of seismic rock physics
montre comment la relation de Gassmann est utilisée dans l'étude sismique des réservoirs.
Smith, T.M., C. H. Sondergeld, and C. S. Rai, 2003,
Gassmann fluid substitutions: A tutorial
Chapman, N.R., J. F. Gandtrust, R. Walia, D. Hannay, G. D. Spence, W. T. Wood, and R.D. Hyndman , 2002, Highresolution, deeptowed, multichannel seismic survey of deepsea gas hydrates off western Canada
Alternatively, one may consider one upgoing wave (sum of all upgoing waves in the layer) and one downgoing wave (sum of all downgoing waves in the layer). The upgoing and downgoing waves in 2 layers having a common interface are linked by continuity conditions. By using iteratively these conditions from the bottom to the top of the stratification, one gets the reflectiontransmission coefficients on the stratification. By computing these coefficients for different frequencies, one gets the seismic trace in time by inverse Fourier transform. To take into account the free surface (R = 1), the upgoing reflected wave is sent back in the stratification. This is done iteratively to take into account multiple reflections on the free surface.
layer
R_{interface} 
downgoing waves downgoing
D_{j} = T_{j+1}D_{j+1}e^{ikzjhj}R_{j+1}M_{j}e^{2ikzjhj} D_{j+1} = (D_{j}e^{ikzjhj}+R_{j+1}M_{j}e^{ikzjhj})/T_{j+1} 
upgoing waves
M_{j+1} = R_{j+1}D_{j+1}+(1R_{j+1})M_{j}e^{ikzjhj} M_{j+1} = (R_{j+1}D_{j}e^{ikzjhj}+M_{j}e^{ikzjhj})/T_{j+1} 
reflection coef. on stratification
R_{h1..j} = M_{j+1}/D_{j+1} R_{j+1} = (R_{h1..j}e^{ikzjhj}R_{h1..j1}e^{ikzjhj})/(e^{ikzjhj}R_{h1..j1}R_{h1..j}e^{ikzjhj}) 
transmission coef. on stratification
T_{h1..j} = D_{0}/D_{j+1} 
layer 3
R_{3}  D_{3} = (e^{ikz2h2}+R_{3}R_{h1}e^{ikz2h2})D_{2}/T_{3}  M_{3} = (R_{3}e^{ikz2h2}+R_{h1}e^{ikz2h2})D_{2}/T_{3} 
R_{h12} = M_{3}/D_{3}
R_{3} = (R_{h12}e^{ikz2h2}R_{h1}e^{ikz2h2})/(e^{ikz2h2}R_{h1}R_{h12}e^{ikz2h2})  T_{h12} = D_{0}/D_{3} = T_{h1}D_{2}/D_{3} 
layer 2
R_{2}  D_{2} = (e^{ikz1h1}+R_{2}R_{1}e^{ikz1h1})D_{1}/T_{2}  M_{2} = (R_{2}e^{ikz1h1}+R_{1}e^{ikz1h1})D_{1}/T_{2} 
R_{h1} = M_{2}/D_{2}
R_{2} = (R_{h1}e^{ikz1h1}R_{1}e^{ikz1h1})/(e^{ikz1h1}R_{1}R_{h1}e^{ikz1h1})  T_{h1} = D_{0}/D_{2} = T_{1}D_{1}/D_{2} 
layer 1
R_{1}  D_{1}  M_{1}  R_{1} = M_{1}/D_{1}  T_{1} = D_{0}/D_{1} 
1/2 space  D_{0}  M_{0} = 0 
stratif.m 
For PSV waves, there are two upgoing and two downgoing waves in each layer (φ and ψ_{y} potentials)
downgoing P and SV waves
D_{Pj} = (T_{PPj+1}D_{Pj+1}+T_{SPj+1}D_{Sj+1}+R^{PPj+1}M_{Pj}e^{ikzPjhj}+R^{SPj+1}M_{Sj}e^{ikzSjhj})e^{ikzPjhj} D_{Sj} = (T_{PSj+1}D_{Pj+1}+T_{SSj+1}D_{Sj+1}+R^{PSj+1}M_{Pj}e^{ikzPjhj}+R^{SSj+1}M_{Sj}e^{ikzSjhj})e^{ikzSjhj} 
upgoings P and SV waves
M_{Pj+1} = R_{PPj+1}D_{Pj+1}+R_{SPj+1}D_{Sj+1}+T^{PPj+1}M_{Pj}e^{ikzPjhj}+T^{SPj+1}M_{Sj}e^{ikzSjhj} M_{Sj+1} = R_{PSj+1}D_{Pj+1}+R_{SSj+1}D_{Sj+1}+T^{PSj+1}M_{Pj}e^{ikzPjhj}+T^{SSj+1}M_{Sj}e^{ikzSjhj}  
(T_{PPj+1}T_{SSj+1}T_{PSj+1}T_{SPj+1})D_{Pj+1}
(T_{SPj+1}T_{PSj+1}T_{SSj+1}T_{PPj+1})D_{Sj+1}  T_{SSj+1}e^{ikzPjhj}  T_{SPj+1}e^{ikzSjhj}  (T_{SPj+1}R^{PSj+1}T_{SSj+1}R^{PPj+1})e^{ikzPjhj}  (T_{SPj+1}R^{SSj+1}T_{SSj+1}R^{SPj+1})e^{ikzSjhj}) 
D_{Pj}
D_{Sj} M_{Pj} M_{Sj} 
T_{PSj+1}e^{ikzPjhj}  T_{PPj+1}e^{ikzSjhj}  (T_{PPj+1}R^{PSj+1}T_{PSj+1}R^{PPj+1})e^{ikzPjhj}  (T_{PPj+1}R^{SSj+1}T_{PSj+1}R^{SPj+1})e^{ikzSjhj})  
M_{Pj+1}
M_{Sj+1}  R_{PPj+1}  R_{SPj+1}  T^{PPj+1}e^{ikzPjhj}  T^{SPj+1}e^{ikzSjhj} 
D_{Pj+1}
D_{Sj+1} M_{Pj} M_{Sj} 
R_{PSj+1}  R_{SSj+1}  T^{PSj+1}e^{ikzPjhj}  T^{SSj+1}e^{ikzSjhj}  
1 (0)
0 (1) R_{PPh1.n} (R_{SPh1.n}) R_{PSh1.n} (R_{SSh1.n})  DM(4,4) = ∏ transfert matrices in layers 1 à n 
1
0 R_{PP1} R_{PS1} 
0
1 R_{SP1} R_{SS1} 
D_{P1}
D_{S1}  
liquidsolid case
D_{Pj} = (T_{PPj+1}D_{Pj+1}+R^{PPj+1}M_{Pj}e^{ikzPjhj}+R^{SPj+1}M_{Sj}e^{ikzSjhj})e^{ikzPjhj}  M_{Pj+1} = R_{PPj+1}D_{Pj+1}+T^{PPj+1}M_{Pj}e^{ikzPjhj}+T^{SPj+1}M_{Sj}e^{ikzSjhj}  
1
R_{PPeau.h1.n}  DM_{liqsol}(2,4)  DM(4,4) 
1
0 R_{PP1} R_{PS1} 
0
1 R_{SP1} R_{SS1} 
D_{P1}
D_{S1} 
Exemples of potentials computed with a densityvelocity model and incidence angles in water θ = 0°, 15° and 30° in cases corresponding to
land and marine reflection seismics (recording at the surface or at the bottom of the sea).




stratifpsv.m
Liu Y. and D. R. Schmitt, 2003, Amplitude and AVO responses of a single thin bed
Chapman, C. H., 2003, Yet another elastic planewave, layermatrix algorithm
give a more elaborate program, that can be used for incidence angles larger than the critical value.




monitorpsv.m
1 path with (T_{e}T_{s})^{m+1}  T_{S0} = T_{e}e^{ikzH}T_{s}∏(e^{ikzH}T_{e}e^{ikzH}T_{s}) = (1R_{e}^{2})^{m+1}e^{(2m+1)ikzH} 
transtrat.m 
N_{1}(2m+1) = 2m+1 paths with R^{2}  T_{S1} = T_{S0}N_{1}(2m+1)R_{e}^{2}e^{2ikzH}  
N_{2}(2m+1) = N_{1}(2m+1)+...+N_{1}(1) paths with R^{4}
N_{1}(2m) paths with R^{2},T_{e}T_{s} 
T_{S2} = T_{S0}(N_{2}(2m+1)R_{e}^{4}
N_{1}(2m)(1R_{e}^{2})R_{e}^{2})e^{4ikzH}  
N_{3}(2m+1) = N_{2}(2m+1)+...+N_{2}(1) paths with R^{6}
2N_{2}(2m+1)2 paths with R^{4},T_{e}T_{s} N_{1}(2m1) paths with R^{2},(T_{e}T_{s})^{2} 
T_{S3} = T_{S0}(N_{3}(2m+1)R_{e}^{6}
(2N_{2}(2m+1)2)(1R_{e}^{2})R_{e}^{4} +N_{1}(2m1)(1R_{e}^{2})^{2}R_{e}^{2})e^{6ikzH}  
N_{4}(2m+1) = N_{3}(2m+1)+...+N_{3}(1) paths with R^{8}
3N_{3}(2m+1)+N_{2}(2m)3 paths with R^{6},T_{e}T_{s} 3N_{2}(2m+1)N_{2}(2)2 paths with R^{4},(T_{e}T_{s})^{2} N_{1}(2m2) paths with R^{2},(T_{e}T_{s})^{3} 
T_{S4} = T_{S0}(N_{4}(2m+1)R_{e}^{8}
(3N_{3}(2m+1)+N_{2}(2m)3)(1R_{e}^{2})R_{e}^{6} +(3N_{2}(2m+1)N_{2}(2)2)(1R_{e}^{2})^{2}R_{e}^{4} +N_{1}(2m2)(1R_{e}^{2})^{3}R_{e}^{2})e^{8ikzH} 
Les termes successifs des coefficients de transmission and reflection sont calculables sans avoir à dénombrer tous les trajands possibles en utilisant les récurrences reliant les waves upgoings M_{j,n+1} and downgoings D_{j,n+1} de part and d'autre de chaque interface j après n+1 reflectiontransmissions à celles provenant des interfaces immédiatement au dessus D_{j1,n} and en dessous M_{j+1,n} après n reflectiontransmissions. On obtient ainsi les réponses impulsionelle and harmonique complètes de la stratification. Les réponses impulsionnelles en reflection and transmission sont caractérisées par des codas avec une amplitude décroissant exponentiellement. Les réponses harmoniques sont caractérisées par des pics de réflectivité (trous de transmissivité) pour H=(2n+1)λ_{z}/4 .
initialisation  reflectiontransmission on odd interface  reflectiontransmission on even interface 
refltranstrat.m 
M_{1,1} = R_{e} , D_{1,1} = T_{e}
M_{2,2} = R_{s}T_{e}e^{ikzH} , D_{2,2} = T_{s}T_{e}e^{ikzH} 
R_{2n+1} = M_{1,2n+1} = T_{s}M_{2,2n}e^{ikzH}
D_{1,2n+1} = R_{s}M_{2,2n}e^{ikzH} M_{2j+1,2n+1} = (R_{e}D_{2j,2n}+T_{s}M_{2j+2,2n})e^{ikzH} D_{2j+1,2n+1} = (T_{e}D_{2j,2n}+R_{s}M_{2j+2,2n})e^{ikzH} 
M_{2j,2n+2} = (R_{s}D_{2j1,2n+1}+T_{e}M_{2j+1,2n+1})e^{ikzH}
D_{2j,2n+2} = (T_{s}D_{2j1,2n+1}+R_{e}M_{2j+1,2n+1})e^{ikzH} M_{2m+2,2n+2} = R_{s}D_{2m+1,2n+1}e^{ikzH} T_{2n+2} = D_{2m+2,2n+2} = T_{s}D_{2m+1,2n+1}e^{ikzH} 
Anstey N.A. and R. F. O'Doherty, 2002, Cycles, layers, and reflections: Part 1
explain the origin on reflections in stratified sedimentary basin as an effect of cyclical stratification.
Stewart, R.R., P. D. Huddleston, and T. K. Kan, 1984,
Seismic versus sonic velocities: A vertical seismic profiling study
explain the difference between local velocity measurements in wells and velocities measured from the surface as an effect of cyclical stratification.
SH case  incident wave at top (a)  incident wave at bottom (b)  reversal (a) + propagation  reversal (b) + propagation  
propagation direction  ↓  ↑  ↓  ↑  ↓  ↑  ↓  ↑ 
top  1  R_{h}  0  t_{h}  R_{h}*  1 = R_{h}R_{h}*+t_{h}T_{h}*  t_{h}*  0 = R_{h}t_{h}*+t_{h}r_{h}* 
bottom  T_{h}  0  r_{h}  1  0 = T_{h}R_{h}*+r_{h}T_{h}*  T_{h}*  1 = T_{h}t_{h}*+r_{h}r_{h}*  r_{h}* 
Pour le cas d'une couche d'épaisseur H dans un milieu homogène, R_{e} étant le coefficient de reflection SH sur l'interface à l'entrée de la couche, on vérifie :
layer  incident wave at top (a)  incident wave at bottom (b)  reversal (a) + propagation  
top  1  R_{h} = R_{e}(1e^{2ikzH})/(1R_{e}^{2}e^{2ikzH})  0  t_{h} = T_{h}  R_{h}R_{h}*+t_{h}T_{h}* = (R_{e}^{2}(1e^{2ikzH})(1e^{2ikzH})+ (1R_{e}^{2})^{2}) /((1R_{e}^{2}e^{2ikzH})(1R_{e}^{2}e^{2ikzH})) = 1 
bottom  T_{h} = (1R_{e}^{2})e^{ikzH}/(1R_{e}^{2}e^{2ikzH})  0  r_{h} = R_{h}  1  T_{h}R_{h}*+r_{h}T_{h}* = R_{e}(1R_{e}^{2})(e^{ikzH}(1e^{2ikzH})+e^{ikzH}(1e^{2ikzH})) /((1R_{e}^{2}e^{2ikzH})(1R_{e}^{2}e^{2ikzH})) = 0 
Par ailleurs, les flux d'énergie entrant and sortant de la stratification doivent être égaux.
SH impedance x cos(incidence)  incident wave at top  incident wave at bottom  
energy flux  entering  outgoing  entering  outgoing 
top : I_{a} = (ρV_{S}cosη)_{haut}  I_{a}  I_{a}R_{h}R_{h}*  0  I_{a}t_{h}t_{h}* 
bottom : I_{b} = (ρV_{S}cosη)_{bas}  0  I_{b}T_{h}T_{h}*  I_{b}  I_{b}r_{h}r_{h}* 
flux conservation  I_{a}(1R_{h}R_{h}*) = I_{b}T_{h}T_{h}*  I_{b}(1r_{h}r_{h}*) = I_{a}t_{h}t_{h}*  
transmission/reflection relations  t_{h}T_{h}* = 1R_{h}R_{h}* = 1r_{h}r_{h}* = (I_{b}/I_{a})T_{h}T_{h}*  t_{h} = (I_{b}/I_{a})T_{h}  R_{h}R_{h}* = r_{h}r_{h}* 
Si la stratification est limitée par une surface libre, il y a reflection totale sur celleci. Pour une onde incidente downgoing, la réponse en reflection de la stratification est renvoyée vers le bas. L'autocorrélation de la réponse en transmission downgoing D dans le demiespace sous la stratification est proportionnelle à la somme de l'onde incidente, de la réponse en reflection and de sa randournée temporelle. De même, pour une onde incidente upgoing, l'autocorrélation de la réponse en transmission upgoing M dans la couche sous la surface libre est proportionnelle à la somme de l'onde incidente, de la réponse en reflection and de sa randournée temporelle.
incident wave  transmission + 1 reflection on free surface  autocorrelation  reflection + time reversal 
downgoing in layer  D = T_{h}(1+R_{h})  DD* = T_{h}T_{h}*(1+R_{h})(1+R_{h})* = (I_{a}/I_{b})T_{h}t_{h}*(1+R_{h})(1+R_{h})* = (I_{a}/I_{b})(1R_{h}R_{h}*)(1+R_{h})(1+R_{h})*  DD* ≈ (I_{a}/I_{b})(1+R_{h}+R_{h}*) 
upgoing in halfspace  M = t_{h}(1+R_{h})  MM* = t_{h}t_{h}*(1+R_{h})(1+R_{h})* = (I_{b}/I_{a})T_{h}t_{h}*(1+R_{h})(1+R_{h})* = (I_{b}/I_{a})(1R_{h}R_{h}*)(1+R_{h})(1+R_{h})*  MM* ≈ (I_{b}/I_{a})(1+R_{h}+R_{h}*) 
incident+reflected harmonic wave on free surface z=0  pressure = 0 on free surface
z = h(x) = He^{ikhx}  wave = p_{0} + scattered wave  P_{s}<<P_{0}
k_{z0}H<<1  harmonic wave scattered in θ_{s} direction 
p_{0}(x,z) = P_{0}e^{ikx0x}(e^{ikz0z}e^{ikz0z})
p_{0}(x,0) = 0  p(x,h(x)) = 0 p(x,0)+h(x)∂p/∂z(x,0) ≈ 0  p(x,z) = (p_{0}+p_{s})(x,z) p_{s}(x,0) = h(x)∂(p_{0}+p_{s})/∂z(x,0)  p_{s}(x,0) ≈ h(x)∂p_{0}/∂z(x,0) p_{s}(x,0) =P_{s}e^{ikxsx} P_{s} = 2ik_{z0}HP_{0} k_{xs} = (k_{x0}+k_{h})  p_{s}(x,z) = P_{s}e^{ikxsx}e^{ikzsz}
sinθ_{s} = sinθ_{0}+Vk_{h}/ω k_{zs} = (ω^{2}/V^{2}k_{xs}^{2})^{½} 