JMM ex EOST

Seismic : Waves - Imaging - Processing - Propagation - Applications

Elastic waves : Characteristics - Elasticity - Wave equations - Plane waves - Spherical waves
Geometrical seismics : Wavefronts and rays - Point source and head wave - V(z) - V(x,y,z)
Reflection-transmission : Continuity - SH - P/SV at free surface - P/SV at interface - Approximations
Guided waves : Rayleigh by free surface and Stoneley by interface - Love and Rayleigh in layer - Group velocity - Love in stratification
Elastic waves + : Cylindrical waves - P and S waves from point force - Anisotropy - Attenuation - Fluid saturation
Reflection-transmission + : Velocity gradient - Layer and stratification - Seismic monitoring - Cyclical stratification - Reflection-transmission relations - Wavy surface
Documents : Figures, programs, animations - Quizs - Books and links

Introduction

Seismic waves are elastic waves that propagate through the Earth. To get acquainted with the way they propagate, visit the following sites :

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.

Characteristics

Seismic waves propagate as undulations (λ = VT), move particles around their equilibrium positions, strain et stress the elastic medium of propagation. Physical quantities involved in the propagation and their magnitude range in geophysics are the following (vectors are in bold).

quantity unit notation components magnitude range
wavelength m λ . 10^-3 m (grains in rocks) - 10⁷ m (size of the Earth)
period
frequency
angular frequency s
Hz
rad/s T
f = 1/T
ω = 2πf . 10⁶ Hz (ultrasonic waves on rock samples) - 10^-3 Hz (free oscillations of the Earth)
propagation velocities
P and S waves m/s V_P
V_S 10² m/s (V_S in soils) - 10⁴ m/s (V_P in the mantle)
V_P/V_S = 1.7 in consolidated rocks
density
impedances kg/m³
kg/m²/s ρ
I = ρV 10³ kg/m³ (water) - 10⁴ kg/m³ (Earth's core)
particle motion m u u_x , u_y , u_z U : 1 m (motion on a fault) - 10^-8m (ambient noise)
particle velocity m/s ∂u/∂t ∂u_x/∂t , ∂u_y/∂t , ∂u_z/∂t U/T
particle acceleration m/s² ∂²u/∂t² ∂²u_x/∂t² , ∂²u_y/∂t² , ∂²u_z/∂t² U/T² : 10 m/s²≈1g (at the epicenters of large earthquakes) - 10^-11m/s² (free oscillations of the Earth)
strains . ε ε_xx , ε_yy , ε_zz
ε_xy , ε_yz , ε_zx U/λ : 10^-6
stresses Pa σ σ_xx , σ_yy , σ_zz
σ_xy , σ_yz , σ_zx ρV²U/λ = IU/T : 10⁴ - 1 Pa (pressure recorded in water in seismic reflection experiments, f = 10-100Hz)
direction of propagation . e e_x , e_y , e_z
sinθ , 0 , cosθ seismic sources emit in wide angular ranges
wave vector m^-1 k = (2π/λ)e k_x , k_y , k_z characterizes an harmonic plane wave wavelength and direction of propagation
quality factor . Q . ∞ - 10 (unconsolidated sediments)

Linear isotropic elasticity

Rocks are springs. Strains due seismic waves are small enough (except close to sources) for stresses to be proportional to strains. The coefficients of proportionality are the elastic moduli. An isotropic medium (behaviour independant of the orientation of applied stresses) is characterized by two elastic moduli. Uniaxial compression involves the Young modulus and the Poisson coefficient of a rock sample. Lamé coefficients λ and μ are used to characterize any states of stress and strain. The bulk modulus K is the ratio of mean pressure to volumetric strain.

. 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 = (1-2ν)ε_xx
stress-strain relation coef. Lamé λ , μ
bulk modulus K σ_xx = λdivu+2με_xx σ_xy = 2με_xy divu = (σ_xx+σ_yy+σ_zz)/3K
relationships among moduli E = 10¹⁰ - 10¹¹ Pa
ν = 0.2 - 0.4 (rocks) ; 0.5 (water) λ = Eν/((1+ν)(1-2ν)) μ = E/(2(1+ν)) K = λ+2μ/3 = E/(3(1-2ν))

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³) 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

with divu = ε_rr+ε_θθ+ε_zz = (1/r)(∂(ru_r)/∂r+∂u_θ/∂θ)+∂u_z/∂z

Elastic wave equations in an homogeneous isotropic medium

Two kinds of waves propagate in an homogeneous isotropic medium with distinct propagation velocities. Compressional waves, called also P waves, have a longitudinal polarization (the particle motion is colinear with the direction of propagation). They propagate at the highest velocity. They are equivalent for elastic media to acoustic waves (sounds) propagating in air and water. Shear waves, called also S waves, have a transverse polarization (the particle motion is in a plane orthogonal to the direction of propagation). They propagate only in elastic media that can be sheared, not in acoustic media. They propagate more slowly than P waves. For compacted rocks (λ ≈ μ or ν ≈ 0.25), V_P/V_S ≈ 1.7. For uncompacted rocks, V_P/V_S can be much larger (> 4). Mathematically, the correct polarization of the waves is automatically fulfilled by writing the particle motion as the gradient of a scalar potential φ for P waves, or as the curl of a vector potential ψ for S waves.

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 = ∂²φ/∂x² , ε_yy , ε_zz
ε_xy = ∂²φ/∂x∂y , ε_yz , ε_zx
stresses σ_xx = E∂u_x/∂x p = -Kdivu σ_xx = λΔφ+2μ∂²φ/∂x² , σ_yy , σ_zz
σ_xy = 2μ∂²φ/∂x∂y , σ_yz , σ_zx
equilibrium equations ρ∂²u_x/∂t² = ∂σ_xx/∂x ρ∂²u/∂t² = -gradp ρ∂²u_x/∂t² = ∂σ_xx/∂x+∂σ_xy/∂y+∂σ_xz/∂z
ρ∂/∂x(∂²φ/∂t²) = (λ+2μ)∂/∂x(Δφ)
ρgrad(∂²φ/∂t²) = (λ+2μ)grad(Δφ)
wave equations ∂²u_x/∂t² = V_L²∂²u_x/∂x² ∂²p/∂t² = V_A²Δp ∂²φ/∂t² = V_P²Δφ
propagation velocity V_L² = E/ρ V_A² = K/ρ V_P² = (λ+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 non-zero 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 = -∂²ψ_y/∂z∂x = -ε_zz
ε_xz = (-∂²ψ_y/∂z²+∂²ψ_y/∂x²)/2
stresses σ_xy = μ∂u_y/∂x σ_xy = μ∂u_y/∂x , σ_yz σ_xx = -2μ∂²ψ_y/∂z∂x = -σ_zz
σ_xz = μ(-∂²ψ_y/∂z²+∂²ψ_y/∂x²)
equilibrium equations ρ∂²u_y/∂t² = ∂σ_xy/∂x ρ∂²u_y/∂t² = ∂σ_xy/∂x+∂σ_yz/∂z ρ∂²u_x/∂t² = ∂σ_xx/∂x+∂σ_xz/∂z
ρ∂/∂z(∂²ψ_y/∂t²) = μ∂/∂z(∂²ψ_y/∂x²+∂²ψ_y/∂z²)
wave equations ∂²u_y/∂t² = V_S²∂²u_y/∂x² ∂²u_y/∂t² = V_S²Δu_y ∂²ψ_y/∂t² = V_S²Δψ_y
propagation velocity V_S² = μ/ρ V_S² = μ/ρ V_S² = μ/ρ

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²/V_S²-2)/(2V_P²/V_S²-2) 0.32 0.24 0.31 0.22 0.29 0.29 0.25
ρ (kg/m³) 2100 2400 2400 2650 2880 3000 3300

Plane waves

A plane wave propagates in a direction e with a velocity V_P or V_S. The sign of the components of e determines the way the wave propagates in the three directions x, y and z (i.e. towards positive or negative x, y and z). The wave equation does not say anything about the solution of the equation, except its phase. A plane harmonic wave corresponds to the choice of a sinusoidal function for the solution, a period T. It is easy to verify that the polarization of the P wave is longitudinal (u // e), and of the S wave transverse (u ⊥ e).

1D wave P waves S waves
u→ = u(x/V-t)
u← = u(-x/V-t) φ(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/V-t) = 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 : ω² = (k_x²+k_y²+k_z²)V² . 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

Spherical waves

Spherical P waves due to a point compressional source with spherical symmetry are solutions of a 1D wave equation in terms of (rφ(r)). The radial particle motion has two terms : one in 1/r, dominant far from the source, and one in 1/r², dominant near the source. For a harmonic wave, the boundary between the far and the near field is defined by the value of kr, i.e. by the ratio r/λ.

∂²(rφ)/∂t² = V²∂²(rφ)/∂r² particle motion near/far field
φ = (A/r)f(±r/V-t) u = (A/r)e_r(-f/r±f'/V) .
φ = (A/r)e^i(kr-ωt) u = (A/r)e_r(ik-1/r)e^i(kr-ωt) kr>>1 u = (ikA/r)e_re^i(kr-ωt)
kr<<1 u = (-A/r²)e_re^i(kr-ωt)
A = -U₀r₀² for a compressional source with radius r₀<<λ

Δφ(r) = (1/r)∂²(rφ)/∂r² , divr = ∂x/∂x+∂y/∂y+∂z/∂z = 3 , dive_r = div(r/r) = grad(1/r).r+divr/r = 2/r

Wavefronts and rays in an homogeneous isotropic medium

The propagation of a wave from a point M reached by the wave at time t is in a direction e at the propagation velocity V. A wavefront corresponds to the surface T(x) on which the perturbation is localized at time t. Rays represent the paths along wich the wave advances with the propagation velocity of the medium. In a homogeneous isotropic medium, the rays are straight lines orthogonal to the wavefronts. A wavefront sweeps along a direction e' different from the propagation direction e with an apparent velocity V_ae' larger than the propagation velocity : V_ae' = V/cos(angle(e,e')). Apparent velocities in three orthogonal directions are geometrically linked by the eikonal equation. The Snell-Descartes law of reflection and transmission expresses the equality of the apparent velocities of the incident, reflected and transmitted wavefronts along the interface between two homogeneous media, thus allowing continuity of displacements and stresses across the interface when the waves sweep the interface.

. 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)²+(∂T/∂y)²+(∂T/∂z)² = 1/V²
(1/V_ax)²+(1/V_ay)²+(1/V_az)² = 1/V²
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.

Point source at a distance d from a plane interface : direct, reflected and head waves

A spherical wavefont incident on a plane interface generates a reflected spherical wave and a transmitted (non-spherical) wave. If the velocity increases across the interface (V₂>V₁), it is only the part of the incident spherical wave that corresponds to incidence angles less than the critical incidence θ_c that is transmitted. The part corresponding to incidence angles larger than the critical incidence is totally reflected. The transmitted wavefront, initially dragged by the incident wavefront along the interface, propagates horizontally under the interface at the velocity V₂ as soon as θ₂ reaches 90°. Since V₂ is larger than the apparent velocity of the incident and reflected waves at incidences θ₁>θ_c, the transmitted wave is in advance relative to the incident wave. It produces a diffracted wave, called head wave, that travels at the constant incidence angle θ_c in the medium 1. The head wave is the "wake" in the medium 1 of the transmitted wave. Its wavefront has a conical shape.
For a non dipping interface parallel to the recording surface, the apparent velocity of the head wave at the surface is the propagation velocity V₂ of the transmitted wave along the interface.
For a dipping interface making an angle α with the recording surface, the apparent velocity of the head wave at the surface is different if the receivers are located up-dip (V_au = V₁/sin(θ_c-α)) or down-dip (V_ad = V₁/sin(θ_c+α)) relative to the point source. By measuring the up-dip and down-dip apparent velocities in reverse experiments, it is possible to determine V₂ and α.

Traveltime curves T(X) of direct, reflected and head waves for source-receiver offset SR=X and homogeneous media of velocity V₁ and V₂
dip α direct wave reflected wave
total reflection for incidence θ₁>θ_c head wave (only if V₂>V₁)
critical incidence θ_c : sinθ_c = V₁/V₂

direco.m T = X/V₁ T = (4d²+X²)^½/V₁ critical distance : X_c = 2dtgθ_c
apparent horizontal et vertical velocities : V_ax = V₁/sinθ_c = V₂ ; V_az = V₁/cosθ_c
traveltime : T = X/V_ax+2d/V_az = X/V₂+2dcosθ_c/V₁
T = X/V₁ virtual source symmetrical of the source with respect to interface
x_v = 2dsinα ; z_v = 2dcosα
T = (4(dcosα)²+(X±2dsinα)²)^½/V₁
+ for receivers located downdip of the source
- for receivers located updip of the source
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 :
V_aα = V₁/sinθ_c = V₂ ; V_a⊥α = V₁/cosθ_c
distances along and orthogonally to interface :
x_α = Xcosα ; d_⊥α = 2d±Xsinα
traveltime : T = x_α/V_aα+d_⊥α/V_a⊥α
T = Xcosα/V₂+(2d±Xsinα)cosθ_c/V₁ = Xsin(θ_c±α)/V₁+2dcosθ_c/V₁

If the source-receiver 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 half-space, 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₁ of water, the velocity V₂ and tickness d of a sedimentary layer, and the velocity V₃ 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.

Rays, distances and times of propagation for V(z) - V_RMS

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₀+az), the rays are circles centered at the depth z = -V₀/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_itanθ_i = p∑d_iV_i/(1-p²V_i²)^½ T(p) = ∑d_i/cosθ_iV_i = ∑d_i/(1-p²V_i²)^½V_i
τ(p) = T(p)-pX(p) = ∑d_icosθ_i/V_i = ∑d_i(1-p²V_i²)^½/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)/(1-p²V(z)²)^½ T(p) = ∫dz/(1-p²V(z)²)^½V(z)
τ(p) = ∫dz(1-p²V(z)²)^½/V(z)
case V(z) = V₀+az
the rays are circles of radius R(p)=1/(ap) centered at z = -V₀/a where V(z) = 0 sinθ(z) = (V₀/a+z)/R = pV(z) X(p) = R[cosθ] = R[(1-p²V(z)²)^½] T(p) = ∫Rdθ/V(z(θ)) = (1/a)∫dθ/sinθ = (1/a)[Log(tan(θ/2))] = (1/a)[Log(pV(z)/(1+(1-p²V(z)²)^½))]
case V(z) = V₀+az
the wavefronts at times T = -(1/a)[Log(tan(θ₀/2))] are the circles orthogonal to the rays
the equation of the wavefronts is x²+(z-Z_max(T))² = X(T)²

hodovz.m tan(θ₀/2) = e^-aT X(T) = V₀/atanθ₀ = (V₀/2a)(e^aT-e^-aT)
X(T) = (V₀/a)sinh(aT)
Z_max(T) = (V₀/a)(1/sinθ₀-1) = (V₀/2a)(e^aT+e^-aT-2)
Z_max(T) = (V₀/a)(cosh(aT)-1)
V_RMS approximation for small incidence angles

vrms.m θ ≈ 0
p ≈ θ_i/V_i X(p) ≈ p∑d_iV_i T(p) ≈ ∑d_i(1+p²V_i²/2)/V_i = ∑d_i/V_i+(p²/2)∑d_iV_i = T(X=0)+X²/2∑d_iV_i
T(X)² ≈ T(0)²+X²/V_RMS²
V_RMS² = (∑d_iV_i)/(∑d_i/V_i) = (∫V(z)dz)/(∫dz/V(z))
V_RMS(V₀+az) = (a(V₀z+az²/2)/Log(1+az/V₀))½

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 θ₀ 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.

- -
hodovzt.m - hodovzo.m

Lutter W.J., R. D. Catchings and C. M. Jarchow, 1994, An image of the Columbia Plateau from inversion of high-resolution 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 near-borehole 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

Rays, wavefronts, amplitudes for V(x,y,z)

For a stratified medium with interfaces dipping in different azimuths, there is no common normal to the interfaces. Snell's law is written for each interface as a function of the direction of the incident ray e and the normal to the interface n. The reflected and transmitted rays are within the plane (e , n). Let m be the unit vector normal to n in the plane (e , n), θ₁ the angle of incidence relative to n, e_r and e_t the directions of the reflected and transmitted rays, θ₂ the angle (e_t , n). The following relations allow tracing the rays in 3D :

e = -cosθ₁n+sinθ₁m e_r = cosθ₁n+sinθ₁m e_r-e = 2cosθ₁n = -2(e.n)n e_r = e-2(e.n)n
e/V₁ = -cosθ₁/V₁n+sinθ₁/V₁m e_t/V₂ = -cosθ₂/V₂n+sinθ₂/V₂m e_t/V₂-e/V₁ = (cosθ₁/V₁-cosθ₂/V₂)n e_t = V₂/V₁(e-(e.n+((V₁/V₂)²-(1-e.n²))^½)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/R-b/τ
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₀+a_xx+a_zz = V₀+a(xsinα+zcosα)
by a rotation with angle α, one gets back to the case V(Z) = V₀+aZ
seismic rays are circles centered on the straigth line of equation V(x,z) = 0 d(sinθ/V)/ds = -a_x/V² = -asinα/V²
d(cosθ/V)/ds = -a_z/V² = -acosα/V² 1/R = dθ/ds = -a_xcosθ/V+a_zsinθ/V = asin(θ-α)/V

One can derive ray equations directly from the wave equation for a V(x) medium by looking for solutions expressed as quasi-plane harmonic waves with wavefront T(x) and amplitude A(x). One gets terms depending on successive powers of the angular frequency ω. The term in ω², 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 quasi-plane wave term in ω² → eikonal equation term in ω → transport equation ray tube with section dS(s)
Δφ+(ω/V(x))²φ = 0 φ = A(x)e^iω(T(x)-t) gradT.gradT=1/V² 2gradA.gradT+AΔT = div(A²gradT) = 0 ∫∫∫div(A²gradT)dυ = ∫∫(A²/V)e.dσ = 0
A²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, 3-D traveltimes and amplitudes by gridded rays give a simple method of 3D ray-tracing by a finite-difference method.
Sava P. and S. Fomel, M. 2001, 3-D traveltime computation using Huygens wavefront tracing

Continuity of displacement and stresses across a plane interface

The reflection-transmission coefficients for a plane interface between two elastic media are obtained by writing the continuity of the components of particle motion and stresses that apply on both sides of the interface. Stresses that do not apply on the interface do not have to be taken into account, even at the depth of the interface. For a free surface, the stresses are zero outside of the elastic medium : they are also zero on the surface. In contrast, there is no condition on the free surface displacement, as proved by the surface damages due to earthquakes.
The coefficients are computed for plane harmonic waves for which the frequency and the apparent horizontal wavelength are the same for the incident, reflected and transmitted waves (Snell's law). k_z depends for each wave on its kind (P or S) and its incidence angle (angles are linked by Snell's law). The sign of k_z is determined by the vertical propagation direction (upgoing or downgoing) of the wave.
Here are the boundary conditions for different types of interfaces and waves, in the case of waves propagating in the vertical plane (x,z) and incident on a plane horizontal interface :

ω 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 P-SV 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

Reflection of a plane harmonic SH wave on a plane interface

A plane harmonic SH wave propagating in the vertical plane (x,z) and incident on a plane horizontal inteface between two elastic media produces on this interface a perturbation of displacement u_y and shear stress σ_zy. This perturbation sweeps along the interface at the apparent velocity of the incident wave V_S1/sinη₁, where η₁ is the angle of the incident ray with respect to the vertical. It is balanced by the perturbations due to the reflected and transmitted SH waves that sweep along the interface at the same apparent velocity (Snell's law). The displacement u_y and shear stress σ_zy resulting from the incident and reflected waves on one side of the interface are equal to u_y and σ_zy due to the transmitted SH wave on the other side of the interface. This is how one gets the reflection and transmission coefficients, expressed here for the particle motion u_y.
In the case V_S1>V_S2 and η₁>η_c, the transmitted wave becomes evanescent : its amplitude decreases in an exponential way with the distance to the interface. There is a phase shift χ for the totally reflected wave given by the argument of the complex reflection coefficient. The interference between the incident and totally reflected wave produces a wave that is standing vertically and propagating horizontally.

Snell's law :
k_x = ωp = ωsinη₁/V_S1 = ωsinη₂/V_S2 incident downgoing SH wave reflected upgoing SH wave transmitted downgoing SH wave
particle motion SH wave
u_y(x,z,t) = U_ye^ik_zze^i(k_xx-ωt) U_y^I = 1
k_z^I = ωcosη₁/V_S1 U_y^R = R
k_z^R = -ωcosη₁/V_S1 U_y^T = T
k_z^T = ωcosη₂/V_S2
stress on horizontal surface
σ_zy(x,z,t) = μ∂u_y/∂z = ρV_S²ik_zu_y = iωS_zye^ik_zze^i(k_xx-ωt) S_zy^I = ρ₁V_S1cosη₁ S_zy^R = -ρ₁V_S1cosη₁R S_zy^T = ρ₂V_S2cosη₂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
(1-R)ρ₁V_S1cosη₁ = Tρ₂V_S2cosη₂ R = (ρ₁V_S1cosη₁-ρ₂V_S2cosη₂)/(ρ₁V_S1cosη₁+ρ₂V_S2cosη₂)
T = 2ρ₁V_S1cosη₁/(ρ₁V_S1cosη₁+ρ₂V_S2cosη₂)
weak contrast
ρ+δρ , V+δV
δη = tanηδV/V (δp=0) R ≈ -δ(ρVcosη)/2ρVcosη = -(1/2)(δρ/ρ+(1-tg²η)δV/V) T ≈ 1
total reflection
sinη_c=V_S1/V_S2
η₁>η_c , pV_S2>1 k_z^R = -k_z^I = -ω(1-p²V_S1²)^½/V_S1
R = (ρ₁V_S1(1-p²V_S1²)^½-iρ₂V_S2(p²V_S2²-1)^½)/(ρ₁V_S1(1-p²V_S1²)^½+iρ₂V_S2(p²V_S2²-1)^½)
R = e^iχ(p) , χ(p) = -2Arctan(ρ₂V_S2(p²V_S2²-1)^½/ρ₁V_S1(1-p²V_S1²)^½)
(u_y^I+u_y^R)(x,z,t) = (e^{ik_z^Iz}+e^{i(k_z^Rz+χ(p))})e^i(k_xx-ωt) = 2cos(k_z^Iz-χ(p)/2)e^{i(k_xx-ωt+χ(p)/2)} transmitted evanescent wave : k_z^T = iω(p²V_S2²-1)^½/V_S2
T = 2ρ₁V_S1(1-p²V_S1²)^½/(ρ₁V_S1(1-p²V_S1²)^½+iρ₂V_S2(p²V_S2²-1)^½)
T = ||T||e^iχ(p)/2
u_y^T(x,z,t) = ||T||e^{-ωz(p²V_S2²-1)^½/V_S2}e^{i(k_xx-ω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

Reflection of a plane harmonic P or SV wave on a plane free surface

P and SV waves propagating in the vertical plane (x,z) produce a compression σ_zz and a shear stress σ_zx on the free forizontal surface of a half-space. There is therefore a coupling between the two kinds of waves at reflection : an incident P or SV wave is reflected as a P and a SV wave. The stresses σ_zz and σ_zx resulting from the incident wave must go to zero at the free surface. The two resulting equations determine the reflection coefficients, expressed here for the potential of the particle motion.
In the case where a SV wave is incident with an angle η larger than the critical incidence corresponding to the ration V_S/V_P, the reflected P wave becomes evanescent. The SV wave is totally reflected with a phase shift χ' given by the argument of the complex reflection coefficient.

Snell's law :
k_x = ωp = ωsinθ/V_P = ωsinη/V_S incident upgoing P wave reflected downgoing P wave reflected downgoing SV wave
displacement potentials P and SV
φ(x,z,t) = Φ e^{ik_z^Pz}e^i(k_xx-ωt)
ψ_y(x,z,t) = Ψ e^{ik_z^Sz}e^i(k_xx-ωt) Φ^I = 1
k_z^P = -ωcosθ/V_P Φ^R = R_φ
k_z^P = ωcosθ/V_P Ψ^R = R_ψ
k_z^S = ωcosη/V_S
stresses due to P wave on a horizontal surface
σ_zz^P = λΔφ+2μ∂²φ/∂z² = -ρω²(1-2p²V_S²)φ = ρω²S_zze^{ik_z^Pz}e^i(k_xx-ωt)
σ_zx^P = 2μ∂²φ/∂z∂x = -2ρV_S²ωpk_zφ = ρω²S_zxe^{ik_z^Pz}e^i(k_xx-ωt) S_zz = -(1-2p²V_S²)
S_zx = 2V_S²pcosθ/V_P S_zz = -(1-2p²V_S²)R_φ
S_zx = -(2V_S²pcosθ/V_P)R_φ .
stresses due to SV wave on a horizontal surface
σ_zz^S = 2μ∂²ψ_y/∂z∂x = -2ρV_S²ωpk_zψ = ρω²S_zze^{ik_z^Sz}e^i(k_xx-ωt)
σ_zx^S = μ(-∂²ψ_y/∂z²+∂²ψ_y/∂x²) = ρω²(1-2p²V_S²)ψ = ρω²S_zxe^{ik_z^Sz}e^i(k_xx-ωt) . . S_zz = -(2V_S²pcosη/V_S)R_ψ
S_zx = (1-2p²V_S²)R_ψ
incident P wave, stresses = 0 on free surface z = 0
(σ_zz^PI + σ_zz^PR + σ_zz^SR)(z=0) = 0
(σ_zx^PI + σ_zx^PR + σ_zx^SR)(z=0) = 0 (1-2p²V_S²)(1+R_φ)+(2V_S²pcosη/V_S)R_ψ = 0
(2V_S²pcosθ/V_P)(1-R_φ)+(1-2p²V_S²)R_ψ = 0
D = (1-2p²V_S²)²+4V_S⁴p²cosθ/V_Pcosη/V_S
R_φ = (-(1-2p²V_S²)²+4V_S⁴p²cosθ/V_Pcosη/V_S)/D
R_ψ = (-4(1-2p²V_S²)V_S²pcosθ/V_P)/D
incident SV wave, stresses = 0 on free surface z = 0
(σ_zz^SI + σ_zz^PR + σ_zz^SR)(z=0) = 0
(σ_zx^SI + σ_zx^PR + σ_zx^SR)(z=0) = 0 (1-2p²V_S²)R'_φ+(2V_S²pcosη/V_S)(-1+R'_ψ) = 0
-(2V_S²pcosθ/V_P)R'_φ+(1-2p²V_S²)(1+R'_ψ) = 0
R'_φ = (4(1-2p²V_S²)V_S²pcosη/V_S)/D
R'_ψ = (-(1-2p²V_S²)²+4V_S⁴p²cosθ/V_Pcosη/V_S))/D
total reflection of SV wave, sinη > V_S/V_P (η > 36°)
reflected evanescent P wave : k_z^P = iω(p²V_P²-1)^½/V_P R'_ψ = e^iχ'(p)
χ'(p) = -2Arctan((4V_S⁴p²(p²V_P²-1)^½/V_Pcosη/V_S)/(1-2p²V_S²)²)

Reflection of a plane harmonic P or SV wave on a plane interface

On a plane horizontal interface between two elastic half-spaces, the P and SV waves propagating in the vertical plane (x,z) are linked by Snell's law and by the continuity of the components of particle motion u_x and u_z and stresses σ_zz and σ_zx across the interface. An incident P or SV wave produces two P and SV reflected waves and two P and SV transmitted waves. The four continuity equations can be written as the product of a matrix by the vector of reflection-transmission coefficients, expressed here for the potentials of particle motion.
For a liquid-solid interface , there is continuity for u_z and σ_zz across the interface. The shear stress σ_zx in the elastic medium must go to zero at the interface, since there is no shear stress in the liquid.
As soon as a value of critical incidence is reached, the reflection-transmission coefficients become complex quantities. The phase-shift χ' at the total reflection of the SV wave is equal to the argument of the reflection coefficient.

Snell's law : k_x = ωp = ωsinθ₁/V_P₁ = ωsinη₁/V_S₁ = ωsinθ₂/V_P₂ = ωsinη₂/V_S₂
P incident downgoing
Φ^I = 1
k_z^I = ωcosθ₁/V_P1 P reflected upgoing
Φ^R = R_φ
k_z^P = -ωcosθ₁/V_P1 SV reflected upgoing
Ψ^R = R_ψ
k_z^S = -ωcosη₁/V_S1 P transmitted downgoing
Φ^T = T_φ
k_z^P = ωcosθ₂/V_P2 SV transmitted downgoing
Ψ^T = T_ψ
k_z^S = ωcosη₂/V_S2
U_x = p + pR_φ + (cosη₁/V_S1)R_ψ = pT_φ -(cosη₂/V_S2)T_ψ
U_z = cosθ₁/V_P1 -(cosθ₁/V_P1)R_φ + pR_ψ = (cosθ₂/V_P2)T_φ + pT_ψ
S_zz = ρ₁(1-2p²V_S1²) + ρ₁(1-2p²V_S1²)R_φ - (2ρ₁V_S1²pcosη₁/V_S1)R_ψ = ρ₂(1-2p²V_S2²)T_φ + (2ρ₂V_S2²pcosη₂/V_S2)T_ψ
S_zx = -2ρ₁V_S1²pcosθ₁/V_P1 + (2ρ₁V_S1²pcosθ₁/V_P1)R_φ + ρ₁(1-2p²V_S1²)R_ψ = -(2ρ₂V_S2²pcosθ₂/V_P2)T_φ + ρ₂(1-2p²V_S2²)T_ψ
matrix notation
case liquid-solid AR = B ; R = [R_φ , R_ψ , T_φ , T_ψ]^t ; B = [-p , -cosθ₁/V_P1 , -ρ₁(1-2p²V_S1²) , 2ρ₁V_S1²pcosθ₁/V_P1]^t
A_LR_L = B_L ; R_L = [R_φ , T_φ , T_ψ]^t ; B_L = [-cosθ₁/V_P1 , -ρ₁ , 0]^t
A[4,4]
A_L[3,3] p cosη₁/V_S1 -p cosη₂/V_S2
-cosθ₁/V_P1 p -cosθ₂/V_P2 -p
ρ₁(1-2p²V_S1²) (ρ₁) -2ρ₁V_S1²pcosη₁/V_S1 -ρ₂(1-2p²V_S2²) -2ρ₂V_S2²pcosη₂/V_S2
2ρ₁V_S1²pcosθ₁/V_P1 (0) ρ₁(1-2p²V_S1²) 2ρ₂V_S2²pcosθ₂/V_P2 -ρ₂(1-2p²V_S2²)
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) = P-Q
R_φ = (-P-Q)/D
R_ψ =2S/D AR = B ; R = [R_φ , R_ψ , T_φ , T_ψ]^t ;
B = [-p , -cosθ₁/V_P1 , -ρ₁(1-2p²V_S1²) , 2ρ₁V_S1²pcosθ₁/V_P1]^t = [-a₁₁ , a₂₁ , -a₃₁ , a₄₁]^t ;
D = (a₁₁A₁₁+a₃₁A₃₁)-(a₂₁A₂₁+a₄₁A₄₁) ;
DR_φ = -(a₁₁A₁₁+a₃₁A₃₁)-(a₂₁A₂₁+a₄₁A₄₁)
DR_ψ/2 = a₁₁a₂₁(a₃₃a₄₄-a₄₃a₃₄) +a₁₁a₄₁(a₂₃a₃₄-a₃₃a₂₄) +a₃₁a₄₁(a₁₃a₂₄-a₂₃a₁₄)
DT_φ/2 = -(a₁₁a₂₁(a₃₂a₄₄-a₄₂a₃₄) +a₁₁a₄₁(a₂₂a₃₄-a₃₂a₂₄) +a₃₁a₄₁(a₁₂a₂₄-a₂₂a₁₄))
DT_ψ/2 = (a₁₁a₂₁(a₃₂a₄₃-a₄₂a₃₃) +a₁₁a₄₁(a₂₂a₃₃-a₃₂a₂₃) +a₃₁a₄₁(a₁₂a₂₃-a₂₂a₁₃)
P incident upgoing
k_z^I = -ωcosθ₂/V_P2 Ar = b ; r = [t^φ,t^ψ,r^φ,r^ψ]^t
b = [p , -cosθ₂/V_P2 , ρ₂(1-2p²V_S2²) , 2ρ₂V_S2²pcosθ₂/V_P2]^t = [-a₁₃ , a₂₃ , -a₃₃ , a₄₃]^t
D = (a₁₃A₁₃+a₃₃A₃₃)-(a₂₃A₂₃+a₄₃A₄₃) ;
Dr_φ = -(a₁₃A₁₃+a₃₃A₃₃)-(a₂₃A₂₃+a₄₃A₄₃)
Dr_ψ/2 = a₁₃a₂₃(a₃₁a₄₂-a₄₁a₃₂) +a₁₃a₄₃(a₂₁a₃₂-a₃₁a₂₂) +a₃₃a₄₃(a₁₁a₂₂-a₂₁a₁₂)
Dt_φ/2 = -(a₁₃a₂₃(a₃₄a₄₂-a₄₄a₃₂) +a₁₃a₄₃(a₂₄a₃₂-a₃₄a₂₂) +a₃₃a₄₃(a₁₄a₂₂-a₂₄a₁₂))
Dt_ψ/2 = a₁₃a₂₃(a₃₄a₄₁-a₄₄a₃₁) +a₁₃a₄₃(a₂₄a₃₁-a₃₄a₂₁) +a₃₃a₄₃(a₁₄a₂₁-a₂₄a₁₁)
case liquid-solid

reflpsl.m A_LR_L = B_L ; R_L = [R_φL , T_φL , T_ψL]^t ; B_L = [-cosθ₁/V_P1 , -ρ₁ , 0]^t
A_Lr_L = b_L ; r_L = [t_φL , r_φL , r_ψL]^t ; b_L = [-cosθ₂/V_P2 , ρ₂(1-2p²V_S2²) , 2ρ₂V_S2²pcosθ₂/V_P2]^t
D_L = -ρ₂²cosθ₁/V_P1((1-2p²V_S2²)²+4V_S2⁴p²cosθ₂/V_P2cosη₂/V_S2) -ρ₁ρ₂cosθ₂/V_P2
D_LR_φL = -ρ₂²cosθ₁/V_P1((1-2p²V_S2²)²+4V_S2⁴p²cosθ₂/V_P2cosη₂/V_S2) +ρ₁ρ₂cosθ₂/V_P2
D_LT_φL = -2ρ₁ρ₂cosθ₁/V_P1(1-2p²V_S2²) ; D_Lt_φL = -2ρ₂²cosθ₂/V_P2(1-2p²V_S2²)
D_LT_ψL = -4ρ₁ρ₂cosθ₁/V_P1V_S2²pcosθ₂/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η₁/V_S1
D = det(A) = -P'+Q'
R'_ψ = (-P'-Q')/D AR'=B' ; R' = [R'_φ , R'_ψ , T'_φ , T'_ψ]^t ;
B' = [cosη₁/V_S1 , -p , -2ρ₁V_S1²pcosη₁/V_S1 , -ρ₁(1-2p²V_S1²)]^t = [a₁₂ , -a₂₂ , a₃₂ , -a₄₂]^t ;
D = -(a₁₂A₁₂+a₃₂A₃₂)+(a₂₂A₂₂+a₄₂A₄₂) ;
DR'_ψ = -(a₁₂A₁₂+a₃₂A₃₂)-(a₂₂A₂₂+a₄₂A₄₂)
SV incident upgoing
k_z^I = -ωcosη₂/V_S2 Ar'=b' ; r' = [t'_φ , t'_ψ , r'_φ , r'_ψ]^t
b' = [cosη₂/V_S2 , p , -2ρ₂V_S2²pcosη₂/V_S2 , ρ₂(1-2p²V_S2²)]^t = [a₁₄ , -a₂₄ , a₃₄ , -a₄₄]^t ;
D = -(a₁₄A₁₄+a₃₄A₃₄)+(a₂₄A₂₄+a₄₄A₄₄)
Dr'_ψ = -(a₁₄A₁₄+a₃₄A₃₄)-(a₂₄A₂₄+a₄₄A₄₄)
total reflection of downgoing SV wave
p = sinη₁/V_S₁
1/V_S1 > p > 1/V_P1
p > 1/V_S2 > 1/V_P2
R'_ψ = e^iχ'(p) imaginary terms of A : a₁₄ = i(p²V_S2²-1)^½/V_S2 ; a₂₁ = -i(p²V_P1²-1)^½/V_P1 ; a₂₃ = -i(p²V_P2²-1)^½/V_P2 ;
a₃₄ = -2iρ₂V_S2²p(p²V_S2²-1)^½/V_S2 ; a₄₁ = 2iρ₁V_S1²p(p²V_P1²-1)^½/V_P1 ; a₄₃ = 2iρ₂V_S2²p(p²V_P2²-1)^½/V_P2 ;
⇒ A₁₂, A₃₂ and P' are imaginary, A₁₂, A₃₂ and Q' are real
R'_ψ = (-P'-Q')/(-P'+Q') = e^iχ'(p)
χ'(p) = 2Arctan(P'/iQ') = 2Arctan((1/i)(a₁₂A₁₂+a₃₂A₃₂)/(a₂₂A₂₂+a₄₂A₄₂))

Here are the moduli of the reflection-transmission coefficients P-SV and SV-P computed for incidence angles of 0°, 15° and 30° in water for a sedimentary stratified medium :
-

Approximations of P-P and P-SV reflection coefficients for weak contrasts

When the contrast in elastic parameters and density across the interface is weak, the reflection coefficients can be approximated with linear relations in terms of these contrasts. These approximations are used in the analysis of reflected amplitudes in reflection seismic prospecting.

weak contrasts ρ₂ = ρ+δρ ; V_P₂ = V_P+δV_P ; θ₂ = θ+δθ ; δθ = tanθδV_P/V_P ; V_S₂ = V_S+δV_S ; η₂ = η + δη ; δη = tanηδV_S/V_S
a₁₁ = -a₁₃ = a₂₂ = -a₂₄ = p ; a₁₂ = cosη/V_S ; a₂₁ = -cosθ/V_P ; a₃₁ = a₄₂ = ρ(1-2p²V_S²) ; a₃₂ = -2ρV_S²pcosη/V_S ; a₄₁ = 2ρV_S²pcosθ/V_P
a₁₄ = (a₁₂+δa₁₂) ; a₂₃ = (a₂₁+δa₂₁) ; a₃₃ = a₄₄ = -(a₃₁+δa₃₁) ; a₃₄ = (a₃₂+δa₃₂) ; a₄₃ = (a₄₁+δa₄₁)
D = det(A) D = (p²-a₂₁a₁₂)(a₃₁²-a₄₁a₃₂) - (2pa₂₁)(-2a₃₁a₃₂) + (-p²-a₂₁a₁₂)(a₃₂a₄₁+a₃₁²) + (a₁₂a₂₁+p²)(-a₃₁²-a₄₁a₃₂) - (-2pa₁₂)(2a₃₁a₄₁) + (p²-a₂₁a₁₂)(a₃₁²-a₄₁a₃₂)
D = 2(p²-a₂₁a₁₂)(a₃₁²-a₄₁a₃₂) -2(p²+a₂₁a₁₂)(a₃₂a₄₁+a₃₁²) +4pa₃₁(a₂₁a₃₂+a₁₂a₄₁) = -4p²a₄₁a₃₂ -4a₂₁a₁₂a₃₁² +4pa₃₁(a₂₁a₃₂+a₁₂a₄₁)
D = 4cosθ/V_Pcosη/V_S [p²(2ρV_S²p)² +(ρ(1-2p²V_S²))² +4ρ²p²V_S²(1-2p²V_S²)] = 4ρ²cosθ/V_Pcosη/V_S
R_φ DR_φ = (-p²-a₂₁a₁₂)((a₃₁+δa₃₁)²-(a₄₁+δa₄₁)(a₃₂+δa₃₂)) - (-p(a₂₁+δa₂₁)+pa₂₁)(-a₃₂(a₃₁+δa₃₁)-a₃₁(a₃₂+δa₃₂)) + (p²-a₂₁(a₁₂+δa₁₂))(a₃₂(a₄₁+δa₄₁)+a₃₁(a₃₁+δa₃₁)) + (a₁₂(a₂₁+δa₂₁)+p²)(a₃₁(a₃₁+δa₃₁)-a₄₁(a₃₂+δa₃₂)) - (-pa₁₂-p(a₁₂+δa₁₂))(-a₃₁(a₄₁+δa₄₁)+a₄₁(a₃₁+δa₃₁)) + (p²-(a₂₁+δa₂₁)(a₁₂+δa₁₂))(-a₃₁²-a₄₁a₃₂))
DR_φ = (-p²-a₂₁a₁₂)(2a₃₁δa₃₁-a₄₁δa₃₂-a₃₂δa₄₁) - 2pa₃₂a₃₁δa₂₁ + (p²-a₂₁a₁₂)(a₃₂δa₄₁+a₃₁δa₃₁)-a₂₁δa₁₂(a₃₂a₄₁+a₃₁²) + (p²+a₁₂a₂₁)(a₃₁δa₃₁-a₄₁δa₃₂)+a₁₂δa₂₁(a₃₁²-a₄₁a₃₂) + 2pa₁₂(-a₃₁δa₄₁+a₄₁δa₃₁) + (a₂₁δa₁₂+a₁₂δa₂₁)(a₃₁²+a₄₁a₃₂)
DR_φ = 2(p(pa₃₂-a₁₂a₃₁)δa₄₁ + a₁₂(pa₄₁-a₂₁a₃₁)δa₃₁ + a₃₁(a₃₁a₁₂-pa₃₂)δa₂₁)
DR_φ = 2ρ(-2p²cosη/V_Sδ(ρV_S²cosθ/V_P) +cosη/V_Scosθ/V_Pδ(ρ(1-2p²V_S²) + ρ(1-2p²V_S²)cosη/V_Sδ(-cosθ/V_P))
DR_φ = 2ρcosη/V_S(cosθ/V_P((1-4p²V_S²)δρ - 4ρp²δ(V_S²)) - ρδ(cosθ/V_P))
R_φ = (1/2)((1-4p²V_S²)δρ/ρ - 8p²V_S²δV_S/V_S - δ(cosθ/V_P)/(cosθ/V_P))
R_φ = (1/2)((1-4p²V_S²)δρ/ρ - 8p²V_S²δV_S/V_S + (1/cosθ)²δV_P/V_P)
R_φ = (1/2)((δρ/ρ+δV_P/V_P) + (δV_P/V_P -4V_S²/V_P²(δρ/ρ+2δV_S/V_S)sin²θ + (tg²θ-sin²θ)δV_P/V_P) ≈ R₀+Gsin²θ
R_ψ DR_ψ = (2pa₂₁)((a₃₁+δa₃₁)²-(a₄₁+δa₄₁)(a₃₂+δa₃₂)) - p((a₂₁+δa₂₁)+a₂₁)(a₃₁(a₃₁+δa₃₁)-a₄₁(a₃₂+δa₃₂)) + (-p²-a₂₁(a₁₂+δa₁₂))(-a₃₁(a₄₁+δa₄₁)+a₄₁(a₃₁+δa₃₁)) + p(-(a₂₁+δa₂₁)+a₂₁)(-a₃₁(a₃₁+δa₃₁)-a₄₁(a₃₂+δa₃₂)) - (p²-a₂₁(a₁₂+δa₁₂))(a₃₁(a₄₁+δa₄₁)+a₄₁(a₃₁+δa₃₁)) + (p²-(a₂₁+δa₂₁)(a₁₂+δa₁₂))(2a₃₁a₄₁)
DR_ψ = 2pa₂₁(2a₃₁δa₃₁-a₄₁δa₃₂-a₃₂δa₄₁) - 2pa₂₁(a₃₁δa₃₁-a₄₁δa₃₂)-pδa₂₁(a₃₁²-a₄₁a₃₂) + (-p²-a₂₁a₁₂)(-a₃₁δa₄₁)+a₄₁δa₃₁)) + (-pδa₂₁)(-a₃₁²-a₄₁a₃₂) - (p²-a₂₁a₁₂)(a₃₁δa₄₁+a₄₁δa₃₁) +a₂₁δa₁₂2a₃₁a₄₁ + (-2a₃₁a₄₁(a₂₁δa₁₂+a₁₂δa₂₁))
DR_ψ = 2(p(a₂₁a₃₁-pa₄₁)δa₃₁ + a₂₁(-pa₃₂+a₁₂a₃₁)δa₄₁ + a₄₁(pa₃₂-a₃₁a₁₂)δa₂₁)
DR_ψ = -2ρpcosθ/V_P(δ(ρ(1-2p²V_S²)) + cosη/V_Sδ(2ρV_S²cosθ/V_P) + 2ρV_S²cosη/V_Sδ(-cosθ/V_P))
DR_ψ = -2ρpcosθ/V_P ((1-2p²V_S²+2V_S²cosθ/V_Pcosη/V_S)δρ + 2ρ(-p²+cosη/V_Scosθ/V_P)δ(V_S²)
R_ψ = -(1/2)pV_S/cosη ((1-2p²V_S²+2V_S²cosθ/V_Pcosη/V_S)δρ/ρ + 4V_S²(-p²+cosη/V_Scosθ/V_P)δV_S/V_S)
R_ψ ≈ -(1/2)sinη ((1+2V_S/V_P)δρ/ρ + 4V_S/V_PδV_S/V_S) pour incidence faible

Rayleigh wave guided by the free surface of a half-space and Stoneley wave guided by an interface

The Rayleigh wave corresponds to the propagation along the free surface of a half-space (z ≥ 0 ) of two evanescent P and SV waves having the same horizontal propagation velocity V_R < V_S. By writing that stresses are zero on the free surface, one gets the equation giving the value of V_R, as a function of the V_P/V_S ratio. V_R ≈ 0.9V_S does not depend on frequency. The particle motion at the surface is elliptical because there is a phase shift π/2 between u_x and u_z.

1/V_R = p = k_x/ω > 1/V_S > 1/V_P evanescent P wave
k_z^P = iω(p²V_P²-1)^½/V_P evanescent S wave
k_z^S = iω(p²V_S²-1)^½/V_S
σ_zz = σ_zx = 0 on free surface z = 0 (1-2p²V_S²)Φ+i2pV_S(p²V_S²-1)^½Ψ = 0
-i2V_S²p(p²V_P²-1)^½/V_PΦ+(1-2p²V_S²)Ψ = 0
déterminant = 0
(1-2p²V_S²)²-4p²V_S⁴(p²V_S²-1)^½/V_S(p²V_P²-1)^½/V_P = 0
(2-V_R²/V_S²)²-4(1-V_R²/V_P²)^½(1-V_R²/V_S²)^½ = 0

(u_x , u_z) = (U_x(z) , U_z(z))e^iω(x/V_R-t) U_x = iωpΦe^{-ωz(p²V_P²-1)^½/V_P}+ω(p²V_S²-1)^½/V_SΨe^{-ωz(p²V_S²-1)^½/V_S}
U_z = -ω(p²V_P²-1)^½/V_PΦe^{-ωz(p²V_P²-1)^½/V_P}+iωpΨe^{-ωz(p²V_S²-1)^½/V_S}
U_x = e^iπ/2ωpΦ(e^{-ωz(p²V_P²-1)^½/V_P}+(1/(2p²V_S²)-1)e^{-ωz(p²V_S²-1)^½/V_S})
U_z = e^iπω(p²V_P²-1)^½/V_PΦ(e^{-ωz(p²V_P²-1)^½/V_P}+(1/(2p²V_S²)-1)^-1e^{-ωz(p²V_S²-1)^½/V_S})

In a similar way, an interface (z = 0) between two elastic half-spaces can guide a Stoneley wave (only for some values of velocities and densities). It is made of two coupled Rayleigh waves (two P₁, SV₁ evanescent waves for z ≤ 0 coupled to two P₂ , SV₂ 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 reflection-transmission coefficients for P-SV waves. There is not always a solution , except in the liquid-solid case.

Stoneley equation
liquid-solid interf. D_L = -ρ₂²cosθ₁/V_P1((1-2p²V_S2²)²+4V_S2⁴p²cosθ₂/V_P2cosη₂/V_S2) -ρ₁ρ₂cosθ₂/V_P2 = 0
(1-2p²V_S2²)²-4V_S2⁴p²(p²V_P2²-1)^½/V_P2(p²V_S2²-1)^½/V_S2 + (ρ₁/ρ₂)((p²V_P2²-1)^½/V_P2)/((p²V_P1²-1)^½/V_P1) = 0

Stephen, R.A., F. Cardo-Casas, and C. H. Cheng , 1985, Finite-difference synthetic acoustic logs show the propagation of a Stoneley wave in a well for sonic logging.

Love and Rayleigh waves guided in a layer with thickness H

The Love wave is a guided wave resulting from constructive interferences of SH waves totally reflected in a layer having a propagation velocity V_S1 less than the propagation velocity V_S2 of the half-space beneath the layer. The total reflection on the free surface is without phase-shift. On the interface with the half-space, there is a phase-shift χ , computed for the reflection coefficient of a SH wave for an incidence angle in the slow medium larger than η_c. The propagation velocity V_L of the Love wave is the apparent horizontal velocity of the SH wave in the layer, such that V_S1≤V_L≤V_S2. It depends on frequency (dispersion). The different normal modes (fundamental and harmonics) in the layer correspond to a node of stress at the free surface for the stationary wave imposed by the total reflection at the interface.

standing wave in z
k_z for normal modes travelling wave in x , k_x = (ω²/V_S²-k_z²)^½
phase velocity V(ω) = V_ax = ω/k_x ω ≥ high-frequency cut-off ω_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^ik_zz+e^-ik_zz = 2cos(k_zz)
Hk_{z n} = nπ , nλ_{z n}/2 = H , n = 1,2,3... k_{x n} = (ω²/V_S²-(nπ/H)²)^½
V_n(ω) = ω/k_{x n} ω_{c n} = nπV_S/H
layer (-H ≤ z ≤ 0)
over half-space (z≥ 0)
R_surf = 1 ; R_interf = e^iχ(p) e^ik_1zz+e^iχ(p)e^-ik_1zz = 2e^iχ(p)/2cos(k_1zz-χ(p)/2)
Hk_{1z n}+χ(p)/2 = nπ , n = 0,1,2,3... k_{x n} = (ω²/V_S1²-(nπ-χ(p)/2)²/H²)^½
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²-1/V_S2²)^½)

displov.m k_1z = (ω/V_S1)(1-p²V_S1²)^½
χ(p) = -2Arctan(ρ₂V_S2(p²V_S2²-1)^½/ρ₁V_S1(1-p²V_S1²)^½) = -2Arctan(I₂/I₁)
dispersion relation of Love wave :
ω(p_n) = (-χ(p_n)/2+nπ)V_S1/(H(1-p_n²V_S1²)^½) = (Arctan(I₂/I₁)+nπ)V_S1/(H(1-p_n²V_S1²)^½)

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 phase-shift χ'₁(p) computed for P-SV reflection coefficients on a free surface. The total reflection on the interface with the Half-space produces a phase-shift χ'₂(p) computed for P-SV reflection coefficients at the interface. The dispersion relation is written by taking into account the sum χ'₁₂(p) of these two phase-shifts 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 cut-off frequency f_cn) of the Rayleigh wave guided in the layer.

boundary conditions normal modes in z travelling wave in x cut-off frequency
R_surf = e^iχ'₁(p)
R_interf = e^iχ'₂(p)
χ'₁₂(p) = χ'₁(p)+χ'₂(p) (n-χ'₁₂(p)/2π)λ_{1z n}/2 = H
Hk_{1z n}+χ'₁₂(p)/2 = nπ k_{x n} = (ω²/V_S1²-(nπ-χ'₁₂(p)/2)²/H²)^½
V_{R n}(ω) = ω/k_{x n} = 1/p_n ω_{c n} = (nπ-χ'₁₂(1/V_S2))/2)/(H(1/V_S1²-1/V_S2²)^½)
dispersion rel. of Rayleigh wave
overtones n≥1
V_S1≤V_{R n}≤V_S2 χ'₁(p) = -2Arctan((4V_S1⁴p²(p²V_P1²-1)^½/V_P1(1-p²V_S1²)^½/V_S1)/(1-2p²V_S1²)²) = -2Arctan(S)
χ'₂(p) = -2Arctan(Im(D)/Re(D)) = -2Arctan(I)
tan(Hω(1-p²V_S1²)^½/V_S1) = -tan((χ'₁(p)+χ'₂(p))/2) = (I+S)/(1-IS)
ω(p_n) = (-(χ'₁(p_n)+χ'₂(p_n))/2+nπ)V_S1/(H(1-p_n²V_S1²)^½) = (Arctan((I+S)/(1-IS))+nπ)V_S1/(H(1-p_n²V_S1²)^½)

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 half-space. 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_RS1is the Rayleigh wave velocity guided at the free surface of a half-space of velocity V_S1 and V_RS2 that for a half-space 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 = ω(1-p²V_P1²)^½/V_P1
Z_P1 = e^-ik_zP1H downgoing S1 , k_zS1 = ω(1-p²V_S1²)^½/V_S1
Z_S1 = e^-ik_zS1H upgoing P1 , -k_zP1 upgoing S1 , -k_zS1 evanescent P2 , k_zP2 = iω(p²V_P2²-1)^½/V_P2 evanescent S2 , k_zS2 = iω(p²V_S2²-1)^½/V_S2
z = -H σ_zz(k_zP1²)Z_P1 -2μ₁k_xk_zS1Z_S1 σ_zz(k_zP1²)/Z_P1 2μ₁k_xk_zS1/Z_S1 0 0
-2μ₁k_xk_zP1Z_P1 σ_zx(k_zS1²)Z_S1 2μ₁k_xk_zP1/Z_P1 σ_zx(k_zS1²)/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μ₁k_xk_zS1 σ_zz(k_zP1²) 2μ₁k_xk_zS1 -σ_zz(k_zP2²) 2μ₂k_xk_zS2
-2μ₁k_xk_zP1 σ_zx(k_zS1²) 2μ₁k_xk_zP1 σ_zx(k_zS1²) 2μ₂k_xk_zP2 -σ_zx(k_zS2²)

The high-frequency 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 high-frequency, 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^{-ω(p²V_P1²-1)^½/V_P1H} (0) ; Z_S1 = e^{-ω(p²V_S1²-1)^½/V_S1H} (0) waves guided by interf. in z = 0 and prolongated in z = -H by Z_P1 , Z_S1
z = -H Rayleigh free surf. (1-2p²V_S1²)Z_P1 (0) -i2pV_S1(p²V_S1²-1)^½Z_S1 (0) 0 0
+i2V_S1²p(p²V_P1²-1)^½/V_P1Z_P1 (0) +(1-2p²V_S1²)Z_S1 (0) 0 0
z = 0 pZ_P1 (0) -i(p²V_S1²-1)^½/V_S1Z_S1 (0) Stoneley at interface
i(p²V_P1²-1)^½/V_P1Z_P1 (0) pZ_S1 (0)
ρ₁(1-2p²V_S1²)Z_P1 (0) i2ρ₁V_S1²p(p²V_S1²-1)^½/V_S1Z_S1 (0)
-i2ρ₁V_S1²p(p²V_P1²-1)^½/V_P1Z_P1 (0) ρ₁(1-2p²V_S1²)Z_S1 (0)

Russel D., Waveguide

Group velocity

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ω)x-t = 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)(1-p_n²V_S1²)^½/V_S1) = (ρ₂V_S2(p²V_S2²-1)^½)/(ρ₁V_S1(1-p²V_S1²)^½)

differentiation of dispersion relation :
(H/V_S1)((1-p_n²V_S1²)^½-ωp_nV_S1²(1-p_n²V_S1²)^-½dp_n/dω)/cos²(Hω(1-p_n²V_S1²)^½/V_S1)) = (ρ₂V_S2)/(ρ₁V_S1)(1-p_n²V_S1²)^-½(p_nV_S2²(p_n²V_S2²-1)^-½+p_nV_S1²(p_n²V_S2²-1)^½/(1-p_n²V_S1²))dp_n/dω

dp_n/dω = (H/V_S1)(1-p_n²V_S1²)/
(Hωp_nV_S1+(ρ₂V_S2)/(ρ₁V_S1)cos²(Hω(1-p_n²V_S1²)^½/V_S1)(p_nV_S2²(p_n²V_S2²-1)^-½+p_nV_S1²(p_n²V_S2²-1)^½/(1-p_n²V_S1²)))

Kelly K.R. , 1983, Numerical study of Love wave propagation
Roth M. and K. Holliger, 1999, Inversion of source-generated noise in high-resolution seismic data

Love wave in stratified medium

The dispersion relations can be obtained by stacking layers of thickness h_i above the interface with the half-space z > 0 where there is total reflection with a phase-shift χ. In each layer, the upgoing and downgoing harmonic waves are prolongated upwards by a vertical phase-shift ±k_izh_i. If the upper surface of the layer is the free surface, the total reflection condition gives the Love wave dispersion. If the upper surface of the layer is an interface with another layer, the continuity conditions on displacement u_y and stress σ_yz give the amplitudes of the upgoing and downgoing harmonic waves in the upper layer (the table reads from bottom to top).

z
I_i = ρ_iV_Si²k_iz downgoing wave e^k_izz
k_iz = ω(1/V_Si²-p²)^½ upgoing wave e^-ik_izz continuity of u_y and σ_yz at interface
total reflection if free surface (upgoing = downgoing) : dispersion relation
-(h₃+h₂+h₁)
layer 3 D₃e^-ik_3zh₃ = D₃(a₃-ib₃) M₃e^ik_3zh₃ = M₃(a₃+ib₃) A₃b₃+a₃B₃ = 0
tan(k_3zh₃)+(I₂/I₃)tan(k_2zh₂)+((I₁/I₃)-(I₁/I₂)tan(k_3zh₃)tan(k_2zh₂))tan(k_1zh₁+χ/2) = 0
-(h₂+h₁)
layer 3 D₃ = A₃-iB₃ M₃ = A₃+iB₃ A₃ = A₂a₂-B₂b₂ = a₁a₂-(I₁/I₂)b₁b₂
B₃ = (I₂/I₃)(A₂b₂+B₂a₂) = (I₂/I₃)a₁b₂+(I₁/I₃)b₁a₂
-(h₂+h₁)
layer 2 D₂e^-ik_2zh₂ = D₂(a₂-ib₂) M₂e^ik_2zh₂ = M₂(a₂+ib₂) A₂b₂+a₂B₂ = 0
tan(k_2zh₂)+(I₁/I₂)tan(k_1zh₁+χ/2) = 0
-h₁
layer 2 D₂ = A₂-iB₂ M₂ = A₂+iB₂ A₂ = a₁
B₂ = (I₁/I₂)b₁
-h₁
layer 1 D₁e^-ik_1zh₁ = a₁-ib₁ M₁e^ik_1zh₁ = a₁+ib₁ b₁ = 0
sin(k_1zh₁+χ/2) = 0
0
layer 1 D₁ = e^-iχ/2 M₁ = e^iχ/2 χ(p) = -2Arctan(ρV_S(p²V_S²-1)^½/ρ₁V_S1(1-p²V_S1²)^½) = -2Arctan(I/I₁)
> 0
1/2 space evanescent wave e^-k_zz
k_z = |ω|(p²-1/V²)^½

Cylindrical waves

In the case of cylindrical symetry, harmonic solutions of the wave equation are expressed in terms of Bessel et Hankel functions of order n = 0,1,2... where n corresponds to the azimutal variation e^inθ.
For k_rr>>1, they are similar to sinusoids with amplitudes decreasing as (1/r)^½.

Bessel equation solution derivative
r²d²f/dr²+rdf/dr+ (k_r²r²-n²)f(r) = 0 f(r) = H_n(k_rr) = J_n(k_rr)+iY_n(k_rr) ≈ (2/πk_rr)^½e^{i(k_rr-π/4-nπ/2)}
H_n(ik_rr) ≈ i^-(n+1)(2/πk_rr)^½e^-k_rr
dH_n(k_rr)/dr = k_rH_n-1(k_rr)-(n/r)H_n(k_rr) = (n/r)H_n(k_rr)-k_rH_n+1(k_rr)

In the absence of azimutal variation, the displacement potentials for P et S waves are :

Δφ(r,z) = ∂²φ/∂r²+(1/r)∂φ/∂r+∂²φ/∂z² = (1/V_P²)∂²φ/∂t² φ(r,z) = f(r)e^i(k_zz-ωt) k_p² = ω²/V_P²-k_z² f(r) = AH₀(k_pr)
Δψ_θ(r,z) = (1/V_S²)∂²ψ_θ/∂t² ψ_θ(r,z) = g(r)e^i(k_zz-ωt) k_s² = ω²/V_S²-k_z² g(r) = BH₁(k_sr)

If there is an angular variation with θ, one has :

Δφ(r,θ,z) = Δφ(r,z) + (1/r)²∂²φ/∂θ² = (1/V_P²)∂²φ/∂t² φ(r,θ,z) =f(r)cos(nθ)e^i(k_zz-ωt) k_p² = ω²/V_P²-k_z² f(r) = AH_n(k_pr)
Δψ_r(r,θ,z) = (1/V_S²)∂²ψ_r/∂t² ψ_r(r,z) = g(r)sin(nθ)e^i(k_zz-ωt) k_s² = ω²/V_S²-k_z² g(r) = BH_n+1(k_sr)
Δψ_θ(r,θ,z) = (1/V_S²)∂²ψ_θ/∂t² ψ_θ(r,z) = -g(r)cos(nθ)e^i(k_zz-ωt)
Δψ_z(r,θ,z) = (1/V_S²)∂²ψ_z/∂t² ψ_z(r,z) = h(r)sin(nθ)e^i(k_zz-ωt) h(r) = CH_n(k_sr)

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(k_zz-ωt) n/rH_n(k_pr)-k_pH_n+1(k_pr) -n/rH_n(k_pr) ik_zH_n(k_pr) (n(n-1)/r²-k_p²)H_n(k_pr)+k_p/rH_n+1(k_pr) -n(n-0.5)/r²H_n(k_pr)+nk_p/rH_n+1(k_pr) ik_z(n/rH_n(k_pr)-k_pH_n+1(k_pr))
Be^i(k_zz-ωt) ik_zH_n+1(k_sr) ik_zH_n+1(k_sr) -k_sH_n(k_sr) ik_z(k_sH_n(k_sr)-(n+1)/rH_n+1(k_sr)) ik_z(k_s/2H_n(k_sr)-(n+1)/rH_n+1(k_sr)) (-nk_s/rH_n(k_sr)+(k_s²-k_z²)H_n+1(k_sr))/2
Ce^i(k_zz-ωt) n/rH_n(k_sr) -n/rH_n(k_sr)+k_sH_n+1(k_sr) 0 n(n-1)/r²H_n(k_sr)-nk_s/rH_n+1(k_sr) (-n(n-1)/r²+k_s²/2)H_n(k_sr)-k_s/rH_n+1(k_sr) (ik_zn/2r)H_n(k_sr)

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

P et S waves from a point force

Without body forces, the 3 equilibrium equations expressed in terms of the components of the particle displacement are :
ρ∂²u/∂t² = (λ+μ)graddivu+μΔu where Δu = (Δu_x, Δu_y, Δu_z) = grad(divu)-curl(curl(u))

The P wave equation, ∂²φ/∂t² = V_P²Δφ, is obtained from these 3 equations by writing u = gradφ. Indeed, Δgradφ = gradΔφ as curl(gradφ) = 0. Therefore, ρgrad∂²φ/∂t² = (λ+2μ)grad(Δφ).
The S wave equation, ∂²ψ/∂t² = V_S²Δψ, 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∂²ψ/∂t² = μ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³) 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²)e_r.4πr²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∂²φ/∂t² = (λ+2μ)grad(Δφ)-(F(t)/(4π)(graddiv(e/r)
The P wave equation becomes : ∂²φ/∂t² = V_P²Δφ-(F(t)/4πρ)div(e/r). It includes a source term proportional to (dive/r) = grad(1/r).e = -1/r²e_r.e = -cosθ/r²
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,θ) ∂²φ/∂t² = V_P²Δφ-(F(t)/4πρ)div(e/r) φ(r,θ) = div(Φ(r)e) ∂²Φ/∂t² = V_P²ΔΦ-F(t)/(4πρr) ∂²(rΦ)/∂t² = V_P²∂²(rΦ)/∂r²-F(t)/(4πρ)
F(t) = Fe^-iωt, k = ω/V_P ∂²(rΦ)/∂r²+k²(rΦ) = F/(4πρV_P²) Φ = F/(4πρV_P²k²)(1-e^ikr)/r φ = gradΦ(r).e = cosθ∂Φ/∂r = -Fcosθ/(4πρV_P²k²)((ik-1/r)e^ikr/r+1/r²)
u = cosθ∂²Φ/∂r²e_r-1/r∂Φ/∂rsinθe_θ = cosθ∂²Φ/∂r²e_r+1/r∂Φ/∂r(e-cosθe_r)
u = ( -F/(4πρV_P²k²)) ((cosθ(2/r²-2ik/r-k²)e_r+1/r(ik-1/r)(e-cosθe_r))e^ikr/r) -(1/r³)(3cosθe_r-e)) u_cl = Fcosθe_r/(4πρV_P²)e^ikr/r u_cp = -F/(4πρV_P²)(e-3cosθe_r)((i/kr-1/(kr)²)e^ikr/r+1/r³)

The equilibrium equations with the part of the source term producing S waves are : ρcurl∂²ψ/∂t² = μcurlΔψ+(F(t)/(4π))curlcurl(e/r).
The S wave equation becomes : ∂²ψ/∂t² =V_S²Δψ+(F(t)/(4πρ))curl(e/r). It includes a source term proportional to curl(e/r) = grad(1/r)∧e = -(1/r²)e_r∧e = sinθe_ξ/r²
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,θ) ∂²ψ/∂t² = V_S²Δψ+(F(t)/4πρ)curl(e/r) ψ(r,θ) = -curl(Ψ(r)e) ∂²Ψ/∂t² = V_S²ΔΨ-F(t)/(4πρr) ∂²(rΨ)/∂t² = V_S²∂²(rΨ)/∂r²-F(t)/(4πρ)
F(t) = Fe^-iωt, k = ω/V_S ∂²(rΨ)/∂r²+k²(rΨ) = F/(4πρV_S²) Ψ = F/(4πρV_S²k²)(1-e^ikr)/r ψ(r,θ) = -gradΨ(r)∧e = sinθ∂Ψ/∂re_ξ = -Fsinθe_ξ/(4πρV_S²k²)((ik-1/r)e^ikr/r+1/r²)
u = (1/r)(2cosθe_r∂Ψ/∂r-sinθe_θ∂(r∂Ψ/∂r)/∂r) = (1/r)(2(e+sinθe_θ)∂Ψ/∂r-sinθe_θ(∂Ψ/∂r+r∂²Ψ/∂r²))
u = -F/(4πρV_S²k²) (1/r)((2(e+sinθe_θ)(ik-1/r)-sinθe_θ(-rk²-ik+1/r))e^ikr/r +(1/r²)(2e+3sinθe_θ)) u_cl = -Fsinθe_θ/(4πρV_S²)e^ikr/r u_cp = F/(4πρV_S²)(e-3cosθe_r) ((i/kr-1/(kr)²)e^ikr/r+1/r³)

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

Transverse isotropy

In an anisotropic medium, the characteristics of the seismic waves are more complex : the propagation velocities depend on the direction of propagation, the polarizations of the waves are no longer exactly longitudinal or transverse, the velocities of SH and SV waves are no longer equal (splitting).
A stratified medium made of thin (with respect to the wavelength) isotropic horizontal layers has an anisotropic behaviour similar to a millefeuille pastry : compressed vertically, the layers are in series ; compressed horizontally, the layers are in parallel. Such a composite medium can be replaced by an equivalent homogeneous medium with an anisotropic behaviour characterized by transverse isotropy. The stress-strain relationships involve five elastic coefficients instead of two for the isotropic case.
The coefficients of the equivalent medium can be obtained by considering various stress-strain states allowing the determination of the different series/parallel averages.

isotropic case : 2 elastic coef. λ and μ transverse isotropy (symmetry axis Oz) : 5 elastic coef. example for clays (x10⁹Pa)
σ_xx = c₁₁ε_xx+c₁₂ε_yy+c₁₂ε_zz
σ_xy = (c₁₁-c₁₂)ε_xy
c₁₁ = λ+2μ , c₁₂ = λ σ_xx = c₁₁ε_xx+c₁₂ε_yy+c₁₃ε_zz , σ_yy = c₁₂ε_xx+c₁₁ε_yy+c₁₃ε_zz
σ_zz = c₁₃ε_xx+c₁₃ε_yy+c₃₃ε_zz
σ_yz = 2c₄₄ε_yz , σ_zx = 2c₄₄ε_zx
σ_xy = 2c₆₆ε_xy , c₆₆ = (c₁₁-c₁₂)/2 c₁₁ = 45 , c₁₂ = 10 , c₁₃ = 8
c₃₃ = 28
c₄₄ = 11
c₆₆ = (c₁₁-c₁₂)/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₃₃ε_zz
c₃₃ = 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))²/∑p_i/(λ_i+2μ_i)-∑p_iλ_i²/(λ_i+2μ_i)) [4] σ_xx = c₁₁ε_xx
[1,4] ⇒ c₁₁ = c₃₃(∑p_iλ_i/(λ_i+2μ_i))²+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₄₄ε_yz
c₄₄ = 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₆₆ε_xy
c₆₆ = ∑p_iμ_i , c₁₂ = c₁₁-2c₆₆
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₁₃ε_xx+c₃₃ε_zz = 0
c₁₃ = c₃₃(∑p_iλ_i/(λ_i+2μ_i))

Here is a plot of these different strain-stress 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 quasi-P waves (quasi-longitudinal polarization α : α ≈ θ)
plane quasi-SV waves (quasi-transverse polarization β : β ≈ η-π/2)
u_y(x,z,t) = U_ye^{i(k_xx+k_zz-ωt)}
ρω² = (c₆₆k_x²+c₄₄k_z²)
the dispersion is an ellipse with semi-axes (c₆₆/ρ)^-½ and (c₄₄/ρ)^-½ u(x,z,t) = (U_x,0,U_z)e^{i(k_xx+k_zz-ωt)}
ρω²U_x = (c₁₁k_x²+c₄₄k_z²)U_x+(c₁₃+c₄₄)k_xk_zU_z
ρω²U_z = (c₁₃+c₄₄)k_xk_zU_x+(c₄₄k_x²+c₃₃k_z²)U_z
(ρω²-c₁₁k_x²-c₄₄k_z²)(ρω²-c₄₄k_x²-c₃₃k_z²)-((c₁₃+c₄₄)k_xk_z)² = 0
V_SH² = (1/ρ)(c₆₆sin²η+c₄₄cos²η) V_QP² = (1/2ρ)(c₄₄+c₁₁sin²θ+c₃₃cos²θ+D)
V_QS² = (1/2ρ)(c₄₄+c₁₁sin²θ+c₃₃cos²θ-D)
D² = ((c₁₁-c₄₄)sin²θ-(c₃₃-c₄₄)cos²θ)²+(c₁₃+c₄₄)²sin²2θ
tanα = (c₁₃+c₄₄)sin2θ/2(ρV_QP²-c₁₁sin²θ-c₄₄cos²θ) = (c₁₃+c₄₄)sin2θ/((c₄₄-c₁₁)sin²θ+(c₃₃-c₄₄)cos²θ+D)
tanβ = (c₁₃+c₄₄)sin2θ/2(ρV_QS²-c₁₁sin²θ-c₄₄cos²θ) = (c₁₃+c₄₄)sin2θ/((c₄₄-c₁₁)sin²θ+(c₃₃-c₄₄)cos²θ-D)
tanα.tanβ = ((c₁₃+c₄₄)sin2θ)²/(((c₄₄-c₁₁)sin²θ+(c₃₃-c₄₄)cos²θ)²-D²) = -1 ⇒ β = α-π/2
group velocity ⇔ stationnary phase in ∫∫U(k_x,k_z)e^{i(k_xx+k_zz-ω(k_x,k_z)t)}dk_xdk_z in order to get constructive interferences of harmonic waves in (x,z,t)
∂(k_xx+k_zz-ω(k_x,k_z)t)/∂k_x = 0 , x-(∂ω/∂k_x)t = 0 , V_Gx = ∂ω/∂k_x
∂(k_xx+k_zz-ω(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₆₆/ρω)k_x , V_Gz = (c₄₄/ρω)k_z
V_G² = (1/ρ²ω²)((c₆₆k_x)²+(c₄₄k_z)²)
sinη_g = V_Gx/V_G = (c₆₆k_x)/((c₆₆k_x)²+(c₄₄k_z)²)^½
cosη_g = V_Gz/V_G = (c₄₄k_z)/((c₆₆k_x)²+(c₄₄k_z)²)^½
V_G² = 1/(ρ(sin²η_g/c₆₆+cos²η_g/c₄₄)) V_Gx = ∂(kV)/∂k_x = V∂k/∂k_x+k∂V/∂k_x
k = (k_x²+k_z²)^½ , ∂k/∂k_x = k_x/k = sinθ
k_x = ksinθ , 1 = ∂k/∂k_xsinθ+kcosθ∂θ/∂k_x , ∂θ/∂k_x = cosθ/k
V_Gx = Vsinθ+cosθdV/dθ , V_Gz = Vcosθ-sinθdV/dθ
V_G² = V²+(dV/dθ)²
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² = c₄₄/ρ = β_T² , η = 90° , V_SH² = c₆₆/ρ
γ_T = (V_SH²(90)-V_SH²(0))/2V_SH²(0) = (c₆₆-c₄₄)/2c₄₄
V_SH² = β_T²((c₆₆/c₄₄)sin²η+cos²η) = β_T²(1+2γ_Tsin²η)
V_G² = β_T²((1+2γ_T)/(1+2γ_Tcos²η_g)) θ = 0 , D = (c₃₃-c₄₄) , V_QP² = c₃₃/ρ = α_T² , V_QS² = c₄₄/ρ = β_T²
θ = 90° , D = (c₁₁-c₄₄) , V_QP² = c₁₁/ρ , V_QS² = c₄₄/ρ = β_T²
ε_T = (V_QP²(90)-V_QP²(0))/2V_QP²(0) = (c₁₁-c₃₃)/2c₃₃
δ_T = (1/2α_T)(d²V_QP/dθ²)(0) = ((c₁₃+c₄₄)²-(c₃₃-c₄₄)²)/(2c₃₃(c₃₃-c₄₄))
example for clays (ρ=2380 kg/m³) : α_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+γ_Tsin²η) ε_T<<1 , δ_T<<1 , ε_T ≈ (V_QP(90)-V_QP(0))/V_QP(0)
V_QP ≈ α_T(1+δ_Tsin²θcos²θ+ε_Tsin⁴θ) = α_T(1+δ_Tsin²θ+(ε_T-δ_T)sin⁴θ)
V_QS ≈ β_T(1+(α_T²/β_T²)(ε_T-δ_T)sin²θcos²θ)

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 QP-QSV reflection-transmission 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θ₁/V_QP1(θ₁) = sinθ₂/V_QP2(θ₂) = sinη₁/V_QS1(η₁) = sinη₂/V_QS2(η₂) ;
q^P1(p) = (cosθ₁/V_QP1)(p) ; q^P2(p) = (cosθ₂/V_QP2)(p) ; q^S1(p) = (cosη₁/V_QS1)(p) ; q^S2(p) = (cosη₁/V_QS2)(p) ;
Dispersion relation : (ρ-c₁₁p²-c₄₄q²)(ρ-c₄₄p²-c₃₃q²)-((c₁₃+c₄₄)pq)² = 0
Aq⁴+Bq²+C = 0 ; A = c₃₃c₄₄ ; B = -c₃₃(ρ-c₁₁p²)-c₄₄(ρ-c₄₄p²)-(c₁₃+c₄₄)²p² ; C = (ρ-c₁₁p²)(ρ-c₄₄p²) ; Δ = (B²-4AC)^½
q^P = ((-B-Δ)/2A)^½ ; q^S = ((-B+Δ)/2A)^½ ;
Polarization : α = Arctan((c₁₃+c₄₄)pq^P/(ρ-c₁₁p²-c₄₄(q^P)²)) ; β = π/2+Arctan((c₁₃+c₄₄)pq^S/(ρ-c₁₁p²-c₄₄(q^S)²))
QP incident downgoing + reflected upgoing
U^I = (sinα₁,cosα₁) , U^RP = R_P(sinα₁,-cosα₁)
k_z^I = ωq^P1(p) , k_z^RP = -ωq^P1(p) QSV reflected upgoing
U^RS = R_S(cosβ₁,sinβ₁)
k_z^RS = -ωq^S1(p) QP transmitted downgoing
U^TP = T_P(sinα₂,cosα₂)
k_z^TP = ωq^P2(p) QSV transmitted downgoing
U^TS = T_S(-cosβ₂,sinβ₂)
k_z^TS = ωq^S2(p)
U_x = sinα₁(1+R_P) + cosβ₁R_S = sinα₂T_P -cosβ₂T_S
U_z = cosα₁(1-R_P) + sinβ₁R_S = cosα₂T_P + sinβ₂T_S
Σ_zz = (c₁₃sinα₁p+c₃₃cosα₁q^P1)(1+R_P) + (c₁₃cosβ₁p-c₃₃sinβ₁q^S1)R_S = (d₁₃sinα₂p+d₃₃cosα₂q^P2)T_P + (-d₁₃cosβ₂p+d₃₃sinβ₂q^S2)T_S
Σ_zx = c₄₄(sinα₁q^P1+cosα₁p)(1-R_P) + c₄₄(-cosβ₁q^S1+sinβ₁p)R_S = d₄₄(sinα₂q^P2+cosα₂p)T_P + d₄₄(-cosβ₂q^S2+sinβ₂p)T_S
écriture matricielle : AR = B ; R = [R_P , R_S , T_P , T_S]^t ; B = [-sinα₁ , -cosα₁ , -(c₁₃sinα₁p+c₃₃cosα₁q^P1) , c₄₄(sinα₁q^P1+cosα₁p)]^t
A sinα₁ cosβ₁ -sinα₂ cosβ₂
-cosα₁ sinβ₁ -cosα₂ -sinβ₂
(c₁₃sinα₁p+c₃₃cosα₁q^P1) (c₁₃cosβ₁p-c₃₃sinβ₁q^S1) - (d₁₃sinα₂p+d₃₃cosα₂q^P2) (d₁₃cosβ₂p-d₃₃sinβ₂q^S2)
c₄₄(sinα₁q^P1+cosα₁p) c₄₄(cosβ₁q^S1-sinβ₁p) d₄₄(sinα₂q^P2+cosα₂p) d₄₄(-cosβ₂q^S2+sinβ₂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 9-component 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.

Attenuation

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) ρ∂²u/∂t² = ∂σ/∂x	w_c = (ρ/2)(∂u/∂t)² ∂w_c/∂t = (∂σ/∂x)(∂u/∂t)	w_e = ∫σdε = (E/2)ε² ∂w_e/∂t = σ∂/∂x(∂u/∂t)	∂(w_c+w_e)/∂t = ∂/∂x(σ∂u/∂t) P_t = -σ∂u/∂t	(1/T)∫₀^TP_tdt
u(x,t) = Ucos(kx-ωt)	w_c = (ρ/2)ω²U²sin²(kx-ωt)	w_e = (ρ/2)V²k²U²sin²(kx-ωt)	P_t = ρVω²U²sin²(kx-ωt) = V(w_c+w_e)	(ρ/2)Vω²U²
u = gradφ(x,t) ρ∂²u/∂t² = -gradp	w_c = (ρ/2)(∂u/∂t)² ∂w_c/∂t = -gradp.∂u/∂t	w_e = (K/2)(divu)² ∂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(1-i/2Q)
2π/Q = -(δw)_cycle/w_max π/Q = -(dU/dx)λ/U U(x) = Ue^-αx
α = π/Qλ = k/2Q = ω/2QV u(x,t) = Ue^-αxcos(kx-ωt) u(x,t) = Ue^-αxe^iω(x/V-t) = Ue^{iω(x(1+i/2Q)/V-t)}
u(x,t) = Ue^iω(x/V_Q-t)

linear viscoelasticity energy dissipated per cycle maximum elastic energy qualty factor
σ = σ_mcos(ωt+χ)
ε = ε_mcosωt (δw)_cycle = ∫₀^Tσdε = πσ_mε_msinχ w_max = ∫₀^T/4σdε ≈ -(1/2)σ_mε_mcosχ 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(1-iωτ)ε = Mε
M = E(1+ω²τ²)^½e^-iχ , 1/Q = tanχ = ωτ ρ∂²u/∂t² = E(1-i/Q)∂²u/∂x²
ρω² = E(1-i/Q)(k+iα)² , k₀² = ρω²/E
k₀² = (1-i/Q)(k+iα)² k+iα = k₀(1+i/2Q)
k = k₀ , V = V₀ = ω/k₀
α = k₀/(2Q) = ω²τ/(2V₀) k²+α² = k₀²/(1+1/Q²)^½ ≈ k₀²(1-1/2Q²)
k²-α² = k₀²/(1+1/Q²) ≈ k₀²(1-1/Q²)
V(ω) ≈ V₀(1+3/(8Q²)) = V₀(1+3ω²τ²/8) , α = k₀/(2Q)
(E and η in parallel) in series with E' , σ = M''ε
1/M'' = 1/M+1/E' = 1/E(1-iωτ)+1/E' = (1/E+1/E'-iωτ/E')/(1-iωτ)
1/E'' = 1/E+1/E' , τ'' = τE''/E' , M'' = E''(1-iωτ)/(1-iωτ'') = ||M''||e^-i(χ-χ'')
1/Q = tan(χ-χ'') = ω(τ-τ'')/(1+ω²ττ'')
d(1/Q)/dω = 0 pour ω_r = (ττ'')^-½ and (1/Q)_max = (τ-τ'')/2(ττ'')^½
ρω² = M''(k+iα)² , k''₀² = ρω²/E''
k²+α² = ρω²/||M''|| = k''₀²(1+ω²τ''²)^½/(1+ω²τ²)^½
k²-α² = ρω²ℜ(1/M'') = k''₀²(1+ω²ττ'')/(1+ω²τ²)
k² = ω²/V(ω)² = (k''₀²/2)((1+ω²τ''²)^½/(1+ω²τ²)^½+(1+ω²ττ'')/(1+ω²τ²))

constant Q model (Kjartansson, 1979, JGR, 84, 4737)
M = M₀(iω/ω₀)^2γ = M₀(|ω/ω₀|)^2γe^{iπγsign(ω)}
1/Q = tan(πγ) ρω² = M₀(|ω/ω₀|)^2γe^{iπγsign(ω)}(k+iα)²
k²+α² = k₀²(|ω/ω₀|)^-2γ
k²-α² = k₀²(|ω/ω₀|)^-2γcos(πγ) k = k₀(|ω/ω₀|)^-γcos(πγ/2)
V(ω) = V₀(|ω/ω₀|)^γ/cos(πγ/2)
α = k₀(|ω/ω₀|)^-γsin(πγ/2) = ktan(πγ/2) dispersion pour 1/Q << 1
V(ω)/V₀ ≈ (|ω/ω₀|)^1/πQ = e^{Log(|ω/ω₀|)/πQ}
V(ω)/V₀ ≈ 1+Log(|ω/ω₀|)/πQ

Poroelasticity and velocities in fluid saturated rocks

La porosité Φ est la proportion volumique d'espace poreux Ω_Φ dans un volume élémentaire Ω. Seuls les pores interconnectés dans lesquels un fluide peut s'écouler interviennent dans Φ, qui atteint 30% pour certains grès. Pour une porosité saturée par un fluide, le volume de fluide Ω_f est égal au volume des pores. On considère une déformation volumique ε de la roche due à une variation de la pression de confinement -σ. L'augmentation ζ du contenu volumique en fluide est due à l'augmentation du volume des pores and à la compression volumique du fluide.

Φ = Ω_Φ/Ω = Ω_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 Biot-Willis) and B (coefficient de Skempton) apparaît en considérant les cas d'une roche drainée (pression de fluide constante, incompressibilité K), non-drainé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₀ 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₀ , 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₀ K = (1-α)K₀ .
ε_Φ = ζ/Φ+ε_f = ε
ε_f = -p_f/K_f (α/ΦK)(1-1/B)+1/K_f = 1/K₀ 1/B = 1-(ΦK/α)(1/K₀-1/K_f) = α/(1-K/K_sat)
relation de Gassmann 1/K_sat = (1-α)/K+α/K(1-1/(1-(ΦK/α)(1/K₀-1/K_f)))
1/K_sat = 1/K₀+Φ/(ΦK/α+1/(1/K_f-1/K₀))
1/(1/K_sat-1/K₀) = K/α+1/(Φ(1/K_f-1/K₀)) = 1/(1/K-1/K₀)+1/(Φ(1/K_f-1/K₀))
densité and vitesses ρ = (1-Φ)ρ₀+Φρ_f V_P² = (K_sat+4μ/3)/ρ V_S² = μ/ρ

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

1D reflection on a velocity gradient (ρ=cte)

layers	particle motions	continuity of u and σ on z=0 and z=H
z<0 V=V₀	u(z) = e^iωz/V₀+Re^-iωz/V₀	1+R = AV₀ⁿ+BV₀^m (iω/V₀)(1-R) = a(AnV₀^n-1+BmV₀^m-1) Te^iωH/V₁ = AV₁ⁿ+BV₁^m (iω/V₁)Te^iωH/V₁ = a(AnV₁^n-1+BmV₁^m-1) ⇒ B/A = -V₁^n-m(n-iω/a)/(m-iω/a)
0≤z≤H V(z)=V₀+az	-ρω²u(z) = ∂/∂z(ρV(z)²∂u/∂z) u(z) = A(V₀+az)ⁿ+B(V₀+az)^m n,m = -1/2±(1/4-ω²/a²)^½ n-m = 2(1/4-ω²/a²)^½
z>H V=V₁	u(z) = Te^iωz/V₁
	R = -(V₀^n-m-V₁^n-m)/((n+iω/a)/(n-iω/a)V₀^n-m-(m+iω/a)/(m-iω/a)V₁^n-m) = 1/(2iω/a+(n-m)(V₀^n-m+V₁^n-m)/(V₀^n-m-V₁^n-m))

Chapman, N.R., J. F. Gandtrust, R. Walia, D. Hannay, G. D. Spence, W. T. Wood, and R.D. Hyndman , 2002, High-resolution, deep-towed, multichannel seismic survey of deep-sea gas hydrates off western Canada

Reflection-transmission of a monochromatic wave on a layer or a stratified medium

The reflection (transmission) of a monochromatic plane SH wave on (through) a layer of thickness H within a homogeneous space is made of interferences between the primary reflection (transmission) when the incident wave enters the layer and the multiple reflections that propagate out of the layer . If R_e and T_e=(1+R_e) are the reflection and transmission coefficients when entering the layer, R_s=(-R_e) and T_s=(1-R_e) when leaving the layer, and k_z the vertical wavenumber in the layer, R_H and T_H are the sum of the series of successive multiples. If R_e<<1, only the first term of the series is significant in amplitude :

R_H1 = R_e+1 reflection in layer , R_e<<1 , T_eT_s = 1-R_e² ≈ 1 R_Hn = R_e+n reflections in layer , R_H = R_Hn (n→∞) R_H1 , R_H , T_H1, T_H
R_H1 = R_e+T_ee^ik_zHR_se^ik_zHT_s = R_e(1-(1-R_e²)e^2ik_zH) ≈ R_e(1-e^2ik_zH)
R_H1 = 0 pour k_zH=nπ or λ_z = 2H/n
R_H1 = 2R_e pour k_zH=(2n+1)π/2 or λ_z = 4H/(2n+1) R_Hn = R_e+T_eR_sT_se^2ik_zH(1+∑(R_se^ik_zH)²ⁿ)
R_H = R_e(1-(1-R_e²)e^2ik_zH/(1-R_e²e^2ik_zH)) = R_e(1-e^2ik_zH)/(1-R_e²e^2ik_zH)
T_H1 = T_ee^ik_zHT_s = (1-R_e²)e^ik_zH ≈ e^ik_zH T_Hn = T_ee^ik_zHT_s(1+∑(R_se^ik_zH)²ⁿ)
T_H = (1-R_e²)e^ik_zH/(1-R_e²e^2ik_zH)

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 reflection-transmission 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+1D_j+1e^ik_zjh_j-R_j+1M_je^2ik_zjh_j
D_j+1 = (D_je^-ik_zjh_j+R_j+1M_je^ik_zjh_j)/T_j+1 upgoing waves
M_j+1 = R_j+1D_j+1+(1-R_j+1)M_je^ik_zjh_j
M_j+1 = (R_j+1D_je^-ik_zjh_j+M_je^ik_zjh_j)/T_j+1 reflection coef. on stratification
R_h1..j = M_j+1/D_j+1
R_j+1 = (R_h1..je^-ik_zjh_j-R_h1..j-1e^ik_zjh_j)/(e^-ik_zjh_j-R_h1..j-1R_h1..je^ik_zjh_j) transmission coef. on stratification
T_h1..j = D₀/D_j+1
layer 3
R₃ D₃ = (e^-ik_z2h₂+R₃R_h1e^ik_z2h₂)D₂/T₃ M₃ = (R₃e^-ik_z2h₂+R_h1e^ik_z2h₂)D₂/T₃ R_h12 = M₃/D₃
R₃ = (R_h12e^-ik_z2h₂-R_h1e^ik_z2h₂)/(e^-ik_z2h₂-R_h1R_h12e^ik_z2h₂) T_h12 = D₀/D₃ = T_h1D₂/D₃
layer 2
R₂ D₂ = (e^-ik_z1h₁+R₂R₁e^ik_z1h₁)D₁/T₂ M₂ = (R₂e^-ik_z1h₁+R₁e^ik_z1h₁)D₁/T₂ R_h1 = M₂/D₂
R₂ = (R_h1e^-ik_z1h₁-R₁e^ik_z1h₁)/(e^-ik_z1h₁-R₁R_h1e^ik_z1h₁) T_h1 = D₀/D₂ = T₁D₁/D₂
layer 1
R₁ D₁ M₁ R₁ = M₁/D₁ T₁ = D₀/D₁
1/2 space D₀ M₀ = 0
stratif.m

For P-SV waves, there are two upgoing and two downgoing waves in each layer (φ and ψ_y potentials)

downgoing P and SV waves
D_Pj = (T_PPj+1D_Pj+1+T_SPj+1D_Sj+1+R^PPj+1M_Pje^ik_zPjh_j+R^SPj+1M_Sje^ik_zSjh_j)e^ik_zPjh_j
D_Sj = (T_PSj+1D_Pj+1+T_SSj+1D_Sj+1+R^PSj+1M_Pje^ik_zPjh_j+R^SSj+1M_Sje^ik_zSjh_j)e^ik_zSjh_j upgoings P and SV waves
M_Pj+1 = R_PPj+1D_Pj+1+R_SPj+1D_Sj+1+T^PPj+1M_Pje^ik_zPjh_j+T^SPj+1M_Sje^ik_zSjh_j
M_Sj+1 = R_PSj+1D_Pj+1+R_SSj+1D_Sj+1+T^PSj+1M_Pje^ik_zPjh_j+T^SSj+1M_Sje^ik_zSjh_j
(T_PPj+1T_SSj+1-T_PSj+1T_SPj+1)D_Pj+1
(T_SPj+1T_PSj+1-T_SSj+1T_PPj+1)D_Sj+1
T_SSj+1e^-ik_zPjh_j -T_SPj+1e^-ik_zSjh_j (T_SPj+1R^PSj+1-T_SSj+1R^PPj+1)e^ik_zPjh_j (T_SPj+1R^SSj+1-T_SSj+1R^SPj+1)e^ik_zSjh_j) D_Pj
D_Sj
M_Pj
M_Sj
T_PSj+1e^-ik_zPjh_j -T_PPj+1e^-ik_zSjh_j (T_PPj+1R^PSj+1-T_PSj+1R^PPj+1)e^ik_zPjh_j (T_PPj+1R^SSj+1-T_PSj+1R^SPj+1)e^ik_zSjh_j)
M_Pj+1
M_Sj+1
R_PPj+1 R_SPj+1 T^PPj+1e^ik_zPjh_j T^SPj+1e^ik_zSjh_j D_Pj+1
D_Sj+1
M_Pj
M_Sj
R_PSj+1 R_SSj+1 T^PSj+1e^ik_zPjh_j T^SSj+1e^ik_zSjh_j
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
liquid-solid case
D_Pj = (T_PPj+1D_Pj+1+R^PPj+1M_Pje^ik_zPjh_j+R^SPj+1M_Sje^ik_zSjh_j)e^ik_zPjh_j M_Pj+1 = R_PPj+1D_Pj+1+T^PPj+1M_Pje^ik_zPjh_j+T^SPj+1M_Sje^ik_zSjh_j
1
R_PPeau.h1.n DM_liq-sol(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 density-velocity 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 plane-wave, layer-matrix algorithm give a more elaborate program, that can be used for incidence angles larger than the critical value.

Seismic monitoring in a stratified elastic medium

Sesimic monitoring of reservoirs involves measurements repeated in time to detect the changes on traveltimes and amplitudes due to changes in elastic properties. The following example shows the perturbations on P and SV waves traveltimes and amplitudes for a model with a 40 m thick layer à 1000m depth in which the velocity V_P decreases from 3200 to 2800 m/s, the other parameters being unchanged.

- - - -
monitorpsv.m

Reflection-transmission of a SH wave on a cyclical stratification

Si l'onde est incidente sur une stratification cyclique constituée d'une alternance régulière de 2 couches telles que les coefficients de reflection-transmission à chaque interface soient alternativement (R_e,T_e) and (R_s,T_s) , l'effand des reflections multiples de même ordre sur les différentes interfaces est cumulatif. Pour un grand nombre d'interfaces, les modules des termes du coefficient de transmission tenant compte des premiers multiples dans chaque couche sont supérieurs à celui de la transmission directe. En supposant que k_zH est le même dans les 2 couches and que la stratification comporte 2m+1 couches (2m+2 interfaces), les premiers termes du coefficient de transmission sont :

1 path with (T_eT_s)^m+1 T_S0 = T_ee^ik_zHT_s∏(e^ik_zHT_ee^ik_zHT_s) = (1-R_e²)^m+1e^(2m+1)ik_zH
transtrat.m
N₁(2m+1) = 2m+1 paths with R² T_S1 = T_S0N₁(2m+1)R_e²e^2ik_zH
N₂(2m+1) = N₁(2m+1)+...+N₁(1) paths with R⁴
N₁(2m) paths with -R²,T_eT_s T_S2 = T_S0(N₂(2m+1)R_e⁴
-N₁(2m)(1-R_e²)R_e²)e^4ik_zH
N₃(2m+1) = N₂(2m+1)+...+N₂(1) paths with R⁶
2N₂(2m+1)-2 paths with -R⁴,T_eT_s
N₁(2m-1) paths with R²,(T_eT_s)² T_S3 = T_S0(N₃(2m+1)R_e⁶
-(2N₂(2m+1)-2)(1-R_e²)R_e⁴
+N₁(2m-1)(1-R_e²)²R_e²)e^6ik_zH
N₄(2m+1) = N₃(2m+1)+...+N₃(1) paths with R⁸
3N₃(2m+1)+N₂(2m)-3 paths with -R⁶,T_eT_s
3N₂(2m+1)-N₂(2)-2 paths with R⁴,(T_eT_s)²
N₁(2m-2) paths with -R²,(T_eT_s)³ T_S4 = T_S0(N₄(2m+1)R_e⁸
-(3N₃(2m+1)+N₂(2m)-3)(1-R_e²)R_e⁶
+(3N₂(2m+1)-N₂(2)-2)(1-R_e²)²R_e⁴
+N₁(2m-2)(1-R_e²)³R_e²)e^8ik_zH

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 reflection-transmissions à celles provenant des interfaces immédiatement au dessus D_j-1,n and en dessous M_j+1,n après n reflection-transmissions. 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 reflection-transmission on odd interface reflection-transmission on even interface
refltranstrat.m
M_1,1 = R_e , D_1,1 = T_e
M_2,2 = R_sT_ee^ik_zH , D_2,2 = T_sT_ee^ik_zH R_2n+1 = M_1,2n+1 = T_sM_2,2ne^ik_zH
D_1,2n+1 = R_sM_2,2ne^ik_zH
M_2j+1,2n+1 = (R_eD_2j,2n+T_sM_2j+2,2n)e^ik_zH
D_2j+1,2n+1 = (T_eD_2j,2n+R_sM_2j+2,2n)e^ik_zH M_2j,2n+2 = (R_sD_2j-1,2n+1+T_eM_2j+1,2n+1)e^ik_zH
D_2j,2n+2 = (T_sD_2j-1,2n+1+R_eM_2j+1,2n+1)e^ik_zH
M_2m+2,2n+2 = R_sD_2m+1,2n+1e^ik_zH
T_2n+2 = D_2m+2,2n+2 = T_sD_2m+1,2n+1e^ik_zH

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.

Reflection-transmission relationships for a stratified medium

On obtient des relations entre les coefficients de reflection-transmission pour des waves incidentes par en haut or par en bas sur la stratification en remarquant qu'une inversion du signe du temps inverse le sens de propagation de toutes les waves harmoniques and remplace leurs amplitudes complexes par leurs conjuguées (les randards de phase correspondant au temps de propagation vertical dans la stratification deviennent des avances). La propagation dans la stratification des waves reflecteds and transmitteds auxquelles ce randournement temporel a été préalablement appliqué reconstitue l'onde incidente initiale. En effand, les waves refont le même chemin en sens inverse and arrivent en phase sur les plans limitant la 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_hR_h*+t_hT_h* t_h* 0 = R_ht_h*+t_hr_h*
bottom T_h 0 r_h 1 0 = T_hR_h*+r_hT_h* T_h* 1 = T_ht_h*+r_hr_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(1-e^2ik_zH)/(1-R_e²e^2ik_zH) 0 t_h = T_h R_hR_h*+t_hT_h* = (R_e²(1-e^2ik_zH)(1-e^-2ik_zH)+ (1-R_e²)²) /((1-R_e²e^2ik_zH)(1-R_e²e^-2ik_zH)) = 1
bottom T_h = (1-R_e²)e^ik_zH/(1-R_e²e^2ik_zH) 0 r_h = R_h 1 T_hR_h*+r_hT_h* = R_e(1-R_e²)(e^ik_zH(1-e^-2ik_zH)+e^-ik_zH(1-e^2ik_zH)) /((1-R_e²e^2ik_zH)(1-R_e²e^-2ik_zH)) = 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_Scosη)_haut I_a I_aR_hR_h* 0 I_at_ht_h*
bottom : I_b = (ρV_Scosη)_bas 0 I_bT_hT_h* I_b I_br_hr_h*
flux conservation I_a(1-R_hR_h*) = I_bT_hT_h* I_b(1-r_hr_h*) = I_at_ht_h*
transmission/reflection relations t_hT_h* = 1-R_hR_h* = 1-r_hr_h* = (I_b/I_a)T_hT_h* t_h = (I_b/I_a)T_h R_hR_h* = r_hr_h*

Si la stratification est limitée par une surface libre, il y a reflection totale sur celle-ci. 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 demi-espace 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_hT_h*(1+R_h)(1+R_h)* = (I_a/I_b)T_ht_h*(1+R_h)(1+R_h)* = (I_a/I_b)(1-R_hR_h*)(1+R_h)(1+R_h)* DD* ≈ (I_a/I_b)(1+R_h+R_h*)
upgoing in half-space M = t_h(1+R_h) MM* = t_ht_h*(1+R_h)(1+R_h)* = (I_b/I_a)T_ht_h*(1+R_h)(1+R_h)* = (I_b/I_a)(1-R_hR_h*)(1+R_h)(1+R_h)* MM* ≈ (I_b/I_a)(1+R_h+R_h*)

Scattering of a plane acoustic wave on a free surface with a weak undulation

A free surface with an undulating topography z=h(x) around the plane z = 0 causes a perturbation of the reflection of a plane wave. For a small amplitude undulation, the reflection is the sum of the wave reflected on the horizontal z = 0 plane according to Snell's law and a wave scattered by the undulation. The approximation of the harmonic scattered wave p_s(x,z) is obtained from p₀(x,z) (the wave propagating in the absence of perturbation) by a series expansion of the boundary condition (pressure = 0) on the undulating interface, by assuming P_s<<P₀.

incident+reflected harmonic wave on free surface z=0 pressure = 0 on free surface
z = h(x) = He^ik_hx wave = p₀ + scattered wave P_s<<P₀
k_z0H<<1 harmonic wave scattered in θ_s direction
p₀(x,z) = P₀e^ik_x0x(e^-ik_z0z-e^ik_z0z)
p₀(x,0) = 0
p(x,h(x)) = 0
p(x,0)+h(x)∂p/∂z(x,0) ≈ 0
p(x,z) = (p₀+p_s)(x,z)
p_s(x,0) = -h(x)∂(p₀+p_s)/∂z(x,0)
p_s(x,0) ≈ -h(x)∂p₀/∂z(x,0)
p_s(x,0) =P_se^ik_xsx
P_s = -2ik_z0HP₀
k_xs = (k_x0+k_h)
p_s(x,z) = P_se^ik_xsxe^ik_zsz
sinθ_s = sinθ₀+Vk_h/ω
k_zs = (ω²/V²-k_xs²)^½

Bibliography and links

ROYER D. et DIEULESAINT E., 1996, Elastic waves in solids 1 : Free and Guided Propagation, Masson, Springer

ACHENBACH J.D., 1973, Wave propagation in elastic solids, North-Holland
BREKHOVSKIKH L.M. ; LYSANOV Y.P., 1991, Fundamentals of Ocean Acoustics, Springer
CLAERBOUT J.F., 1976, Fundamentals of Geophysical Data Processing, McGraw-Hill
LAY T. , WALLACE T. C., 1995, Modern global seismology, Academic Press
MAVKO G. , MUKERJI T. , DVORKIN J., 1998, The rock physics handbook. Tools for seismic analysis in porous media, Cambridge
MULLER G., Theory of elastic waves, Samizdat Press
OFFICER C.B., 1958, Introduction to the theory of sound transmission, McGraw Hill
UDIAS A., 1999, Principles of seismology, Cambridge
WHITE J.E., 1965, Seismic waves : radiation, transmission and attenuation, McGraw-Hill

quantity	unit	notation	components	magnitude range
wavelength	m	λ	.	10^-3 m (grains in rocks) - 10⁷ m (size of the Earth)
period frequency angular frequency	s Hz rad/s	T f = 1/T ω = 2πf	.	10⁶ Hz (ultrasonic waves on rock samples) - 10^-3 Hz (free oscillations of the Earth)
propagation velocities P and S waves	m/s	V_P V_S		10² m/s (V_S in soils) - 10⁴ m/s (V_P in the mantle) V_P/V_S = 1.7 in consolidated rocks
density impedances	kg/m³ kg/m²/s	ρ I = ρV		10³ kg/m³ (water) - 10⁴ kg/m³ (Earth's core)
particle motion	m	u	u_x , u_y , u_z	U : 1 m (motion on a fault) - 10^-8m (ambient noise)
particle velocity	m/s	∂u/∂t	∂u_x/∂t , ∂u_y/∂t , ∂u_z/∂t	U/T
particle acceleration	m/s²	∂²u/∂t²	∂²u_x/∂t² , ∂²u_y/∂t² , ∂²u_z/∂t²	U/T² : 10 m/s²≈1g (at the epicenters of large earthquakes) - 10^-11m/s² (free oscillations of the Earth)
strains	.	ε	ε_xx , ε_yy , ε_zz ε_xy , ε_yz , ε_zx	U/λ : 10^-6
stresses	Pa	σ	σ_xx , σ_yy , σ_zz σ_xy , σ_yz , σ_zx	ρV²U/λ = IU/T : 10⁴ - 1 Pa (pressure recorded in water in seismic reflection experiments, f = 10-100Hz)
direction of propagation	.	e	e_x , e_y , e_z sinθ , 0 , cosθ	seismic sources emit in wide angular ranges
wave vector	m^-1	k = (2π/λ)e	k_x , k_y , k_z	characterizes an harmonic plane wave wavelength and direction of propagation
quality factor	.	Q	.	∞ - 10 (unconsolidated sediments)

.	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 = (1-2ν)ε_xx
stress-strain relation	coef. Lamé λ , μ bulk modulus K	σ_xx = λdivu+2με_xx	σ_xy = 2με_xy	divu = (σ_xx+σ_yy+σ_zz)/3K
relationships among moduli	E = 10¹⁰ - 10¹¹ Pa ν = 0.2 - 0.4 (rocks) ; 0.5 (water)	λ = Eν/((1+ν)(1-2ν))	μ = E/(2(1+ν))	K = λ+2μ/3 = E/(3(1-2ν))

	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³)	1	1000	920	2550	2700	2700	3320

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 = ∂²φ/∂x² , ε_yy , ε_zz ε_xy = ∂²φ/∂x∂y , ε_yz , ε_zx
stresses	σ_xx = E∂u_x/∂x	p = -Kdivu	σ_xx = λΔφ+2μ∂²φ/∂x² , σ_yy , σ_zz σ_xy = 2μ∂²φ/∂x∂y , σ_yz , σ_zx
equilibrium equations	ρ∂²u_x/∂t² = ∂σ_xx/∂x	ρ∂²u/∂t² = -gradp	ρ∂²u_x/∂t² = ∂σ_xx/∂x+∂σ_xy/∂y+∂σ_xz/∂z ρ∂/∂x(∂²φ/∂t²) = (λ+2μ)∂/∂x(Δφ) ρgrad(∂²φ/∂t²) = (λ+2μ)grad(Δφ)
wave equations	∂²u_x/∂t² = V_L²∂²u_x/∂x²	∂²p/∂t² = V_A²Δp	∂²φ/∂t² = V_P²Δφ
propagation velocity	V_L² = E/ρ	V_A² = K/ρ	V_P² = (λ+2μ)/ρ

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 = -∂²ψ_y/∂z∂x = -ε_zz ε_xz = (-∂²ψ_y/∂z²+∂²ψ_y/∂x²)/2
stresses	σ_xy = μ∂u_y/∂x	σ_xy = μ∂u_y/∂x , σ_yz	σ_xx = -2μ∂²ψ_y/∂z∂x = -σ_zz σ_xz = μ(-∂²ψ_y/∂z²+∂²ψ_y/∂x²)
equilibrium equations	ρ∂²u_y/∂t² = ∂σ_xy/∂x	ρ∂²u_y/∂t² = ∂σ_xy/∂x+∂σ_yz/∂z	ρ∂²u_x/∂t² = ∂σ_xx/∂x+∂σ_xz/∂z ρ∂/∂z(∂²ψ_y/∂t²) = μ∂/∂z(∂²ψ_y/∂x²+∂²ψ_y/∂z²)
wave equations	∂²u_y/∂t² = V_S²∂²u_y/∂x²	∂²u_y/∂t² = V_S²Δu_y	∂²ψ_y/∂t² = V_S²Δψ_y
propagation velocity	V_S² = μ/ρ	V_S² = μ/ρ	V_S² = μ/ρ

	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²/V_S²-2)/(2V_P²/V_S²-2)	0.32	0.24	0.31	0.22	0.29	0.29	0.25
ρ (kg/m³)	2100	2400	2400	2650	2880	3000	3300

1D wave	P waves	S waves
u→ = u(x/V-t) u← = u(-x/V-t)	φ(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/V-t) = 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)

∂²(rφ)/∂t² = V²∂²(rφ)/∂r²	particle motion	near/far field
φ = (A/r)f(±r/V-t)	u = (A/r)e_r(-f/r±f'/V)	.
φ = (A/r)e^i(kr-ωt)	u = (A/r)e_r(ik-1/r)e^i(kr-ωt)	kr>>1 u = (ikA/r)e_re^i(kr-ωt) kr<<1 u = (-A/r²)e_re^i(kr-ωt) A = -U₀r₀² for a compressional source with radius 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)²+(∂T/∂y)²+(∂T/∂z)² = 1/V² (1/V_ax)²+(1/V_ay)²+(1/V_az)² = 1/V²
spherical	T(x) = r/V	gradT = e_r/V
any surface	T(x)	gradT = e_M/V

Traveltime curves T(X) of direct, reflected and head waves for source-receiver offset SR=X and homogeneous media of velocity V₁ and V₂
dip α	direct wave	reflected wave total reflection for incidence θ₁>θ_c	head wave (only if V₂>V₁) critical incidence θ_c : sinθ_c = V₁/V₂
direco.m	T = X/V₁	T = (4d²+X²)^½/V₁	critical distance : X_c = 2dtgθ_c apparent horizontal et vertical velocities : V_ax = V₁/sinθ_c = V₂ ; V_az = V₁/cosθ_c traveltime : T = X/V_ax+2d/V_az = X/V₂+2dcosθ_c/V₁
	T = X/V₁	virtual source symmetrical of the source with respect to interface x_v = 2dsinα ; z_v = 2dcosα T = (4(dcosα)²+(X±2dsinα)²)^½/V₁ + for receivers located downdip of the source - for receivers located updip of the source	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 : V_aα = V₁/sinθ_c = V₂ ; V_a⊥α = V₁/cosθ_c distances along and orthogonally to interface : x_α = Xcosα ; d_⊥α = 2d±Xsinα traveltime : T = x_α/V_aα+d_⊥α/V_a⊥α T = Xcosα/V₂+(2d±Xsinα)cosθ_c/V₁ = Xsin(θ_c±α)/V₁+2dcosθ_c/V₁

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_itanθ_i = p∑d_iV_i/(1-p²V_i²)^½	T(p) = ∑d_i/cosθ_iV_i = ∑d_i/(1-p²V_i²)^½V_i τ(p) = T(p)-pX(p) = ∑d_icosθ_i/V_i = ∑d_i(1-p²V_i²)^½/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)/(1-p²V(z)²)^½	T(p) = ∫dz/(1-p²V(z)²)^½V(z) τ(p) = ∫dz(1-p²V(z)²)^½/V(z)
case V(z) = V₀+az the rays are circles of radius R(p)=1/(ap) centered at z = -V₀/a where V(z) = 0	sinθ(z) = (V₀/a+z)/R = pV(z)	X(p) = R[cosθ] = R[(1-p²V(z)²)^½]	T(p) = ∫Rdθ/V(z(θ)) = (1/a)∫dθ/sinθ = (1/a)[Log(tan(θ/2))] = (1/a)[Log(pV(z)/(1+(1-p²V(z)²)^½))]
case V(z) = V₀+az the wavefronts at times T = -(1/a)[Log(tan(θ₀/2))] are the circles orthogonal to the rays the equation of the wavefronts is x²+(z-Z_max(T))² = X(T)² hodovz.m	tan(θ₀/2) = e^-aT	X(T) = V₀/atanθ₀ = (V₀/2a)(e^aT-e^-aT) X(T) = (V₀/a)sinh(aT)	Z_max(T) = (V₀/a)(1/sinθ₀-1) = (V₀/2a)(e^aT+e^-aT-2) Z_max(T) = (V₀/a)(cosh(aT)-1)
V_RMS approximation for small incidence angles vrms.m	θ ≈ 0 p ≈ θ_i/V_i	X(p) ≈ p∑d_iV_i	T(p) ≈ ∑d_i(1+p²V_i²/2)/V_i = ∑d_i/V_i+(p²/2)∑d_iV_i = T(X=0)+X²/2∑d_iV_i
V_RMS approximation for small incidence angles vrms.m	θ ≈ 0 p ≈ θ_i/V_i	T(X)² ≈ T(0)²+X²/V_RMS² V_RMS² = (∑d_iV_i)/(∑d_i/V_i) = (∫V(z)dz)/(∫dz/V(z)) V_RMS(V₀+az) = (a(V₀z+az²/2)/Log(1+az/V₀))½

e = -cosθ₁n+sinθ₁m	e_r = cosθ₁n+sinθ₁m	e_r-e = 2cosθ₁n = -2(e.n)n	e_r = e-2(e.n)n
e/V₁ = -cosθ₁/V₁n+sinθ₁/V₁m	e_t/V₂ = -cosθ₂/V₂n+sinθ₂/V₂m	e_t/V₂-e/V₁ = (cosθ₁/V₁-cosθ₂/V₂)n	e_t = V₂/V₁(e-(e.n+((V₁/V₂)²-(1-e.n²))^½)n)

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/R-b/τ
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₀+a_xx+a_zz = V₀+a(xsinα+zcosα) by a rotation with angle α, one gets back to the case V(Z) = V₀+aZ seismic rays are circles centered on the straigth line of equation V(x,z) = 0	d(sinθ/V)/ds = -a_x/V² = -asinα/V² d(cosθ/V)/ds = -a_z/V² = -acosα/V²	1/R = dθ/ds = -a_xcosθ/V+a_zsinθ/V = asin(θ-α)/V

wave equation	quasi-plane wave	term in ω² → eikonal equation	term in ω → transport equation	ray tube with section dS(s)
Δφ+(ω/V(x))²φ = 0	φ = A(x)e^iω(T(x)-t)	gradT.gradT=1/V²	2gradA.gradT+AΔT = div(A²gradT) = 0	∫∫∫div(A²gradT)dυ = ∫∫(A²/V)e.dσ = 0 A²dS/V = cte

ω 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	P-SV 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η₁/V_S1 = ωsinη₂/V_S2	incident downgoing SH wave	reflected upgoing SH wave	transmitted downgoing SH wave
particle motion SH wave u_y(x,z,t) = U_ye^ik_zze^i(k_xx-ωt)	U_y^I = 1 k_z^I = ωcosη₁/V_S1	U_y^R = R k_z^R = -ωcosη₁/V_S1	U_y^T = T k_z^T = ωcosη₂/V_S2
stress on horizontal surface σ_zy(x,z,t) = μ∂u_y/∂z = ρV_S²ik_zu_y = iωS_zye^ik_zze^i(k_xx-ωt)	S_zy^I = ρ₁V_S1cosη₁	S_zy^R = -ρ₁V_S1cosη₁R	S_zy^T = ρ₂V_S2cosη₂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 (1-R)ρ₁V_S1cosη₁ = Tρ₂V_S2cosη₂	R = (ρ₁V_S1cosη₁-ρ₂V_S2cosη₂)/(ρ₁V_S1cosη₁+ρ₂V_S2cosη₂) T = 2ρ₁V_S1cosη₁/(ρ₁V_S1cosη₁+ρ₂V_S2cosη₂)
weak contrast ρ+δρ , V+δV δη = tanηδV/V (δp=0)	R ≈ -δ(ρVcosη)/2ρVcosη = -(1/2)(δρ/ρ+(1-tg²η)δV/V)		T ≈ 1
total reflection sinη_c=V_S1/V_S2 η₁>η_c , pV_S2>1	k_z^R = -k_z^I = -ω(1-p²V_S1²)^½/V_S1 R = (ρ₁V_S1(1-p²V_S1²)^½-iρ₂V_S2(p²V_S2²-1)^½)/(ρ₁V_S1(1-p²V_S1²)^½+iρ₂V_S2(p²V_S2²-1)^½) R = e^iχ(p) , χ(p) = -2Arctan(ρ₂V_S2(p²V_S2²-1)^½/ρ₁V_S1(1-p²V_S1²)^½) (u_y^I+u_y^R)(x,z,t) = (e^{ik_z^Iz}+e^{i(k_z^Rz+χ(p))})e^i(k_xx-ωt) = 2cos(k_z^Iz-χ(p)/2)e^{i(k_xx-ωt+χ(p)/2)}		transmitted evanescent wave : k_z^T = iω(p²V_S2²-1)^½/V_S2 T = 2ρ₁V_S1(1-p²V_S1²)^½/(ρ₁V_S1(1-p²V_S1²)^½+iρ₂V_S2(p²V_S2²-1)^½) T = \|\|T\|\|e^iχ(p)/2 u_y^T(x,z,t) = \|\|T\|\|e^{-ωz(p²V_S2²-1)^½/V_S2}e^{i(k_xx-ωt+χ(p)/2)}

Snell's law : k_x = ωp = ωsinθ/V_P = ωsinη/V_S	incident upgoing P wave	reflected downgoing P wave	reflected downgoing SV wave
displacement potentials P and SV φ(x,z,t) = Φ e^{ik_z^Pz}e^i(k_xx-ωt) ψ_y(x,z,t) = Ψ e^{ik_z^Sz}e^i(k_xx-ωt)	Φ^I = 1 k_z^P = -ωcosθ/V_P	Φ^R = R_φ k_z^P = ωcosθ/V_P	Ψ^R = R_ψ k_z^S = ωcosη/V_S
stresses due to P wave on a horizontal surface σ_zz^P = λΔφ+2μ∂²φ/∂z² = -ρω²(1-2p²V_S²)φ = ρω²S_zze^{ik_z^Pz}e^i(k_xx-ωt) σ_zx^P = 2μ∂²φ/∂z∂x = -2ρV_S²ωpk_zφ = ρω²S_zxe^{ik_z^Pz}e^i(k_xx-ωt)	S_zz = -(1-2p²V_S²) S_zx = 2V_S²pcosθ/V_P	S_zz = -(1-2p²V_S²)R_φ S_zx = -(2V_S²pcosθ/V_P)R_φ	.
stresses due to SV wave on a horizontal surface σ_zz^S = 2μ∂²ψ_y/∂z∂x = -2ρV_S²ωpk_zψ = ρω²S_zze^{ik_z^Sz}e^i(k_xx-ωt) σ_zx^S = μ(-∂²ψ_y/∂z²+∂²ψ_y/∂x²) = ρω²(1-2p²V_S²)ψ = ρω²S_zxe^{ik_z^Sz}e^i(k_xx-ωt)	.	.	S_zz = -(2V_S²pcosη/V_S)R_ψ S_zx = (1-2p²V_S²)R_ψ
incident P wave, stresses = 0 on free surface z = 0 (σ_zz^PI + σ_zz^PR + σ_zz^SR)(z=0) = 0 (σ_zx^PI + σ_zx^PR + σ_zx^SR)(z=0) = 0	(1-2p²V_S²)(1+R_φ)+(2V_S²pcosη/V_S)R_ψ = 0 (2V_S²pcosθ/V_P)(1-R_φ)+(1-2p²V_S²)R_ψ = 0
	D = (1-2p²V_S²)²+4V_S⁴p²cosθ/V_Pcosη/V_S R_φ = (-(1-2p²V_S²)²+4V_S⁴p²cosθ/V_Pcosη/V_S)/D R_ψ = (-4(1-2p²V_S²)V_S²pcosθ/V_P)/D
incident SV wave, stresses = 0 on free surface z = 0 (σ_zz^SI + σ_zz^PR + σ_zz^SR)(z=0) = 0 (σ_zx^SI + σ_zx^PR + σ_zx^SR)(z=0) = 0	(1-2p²V_S²)R'_φ+(2V_S²pcosη/V_S)(-1+R'_ψ) = 0 -(2V_S²pcosθ/V_P)R'_φ+(1-2p²V_S²)(1+R'_ψ) = 0
	R'_φ = (4(1-2p²V_S²)V_S²pcosη/V_S)/D R'_ψ = (-(1-2p²V_S²)²+4V_S⁴p²cosθ/V_Pcosη/V_S))/D
total reflection of SV wave, sinη > V_S/V_P (η > 36°) reflected evanescent P wave : k_z^P = iω(p²V_P²-1)^½/V_P	R'_ψ = e^iχ'(p) χ'(p) = -2Arctan((4V_S⁴p²(p²V_P²-1)^½/V_Pcosη/V_S)/(1-2p²V_S²)²)

Snell's law : k_x = ωp = ωsinθ₁/V_P₁ = ωsinη₁/V_S₁ = ωsinθ₂/V_P₂ = ωsinη₂/V_S₂
P incident downgoing Φ^I = 1 k_z^I = ωcosθ₁/V_P1	P reflected upgoing Φ^R = R_φ k_z^P = -ωcosθ₁/V_P1	SV reflected upgoing Ψ^R = R_ψ k_z^S = -ωcosη₁/V_S1	P transmitted downgoing Φ^T = T_φ k_z^P = ωcosθ₂/V_P2	SV transmitted downgoing Ψ^T = T_ψ k_z^S = ωcosη₂/V_S2
U_x = p	+ pR_φ	+ (cosη₁/V_S1)R_ψ	= pT_φ	-(cosη₂/V_S2)T_ψ
U_z = cosθ₁/V_P1	-(cosθ₁/V_P1)R_φ	+ pR_ψ	= (cosθ₂/V_P2)T_φ	+ pT_ψ
S_zz = ρ₁(1-2p²V_S1²)	+ ρ₁(1-2p²V_S1²)R_φ	- (2ρ₁V_S1²pcosη₁/V_S1)R_ψ	= ρ₂(1-2p²V_S2²)T_φ	+ (2ρ₂V_S2²pcosη₂/V_S2)T_ψ
S_zx = -2ρ₁V_S1²pcosθ₁/V_P1	+ (2ρ₁V_S1²pcosθ₁/V_P1)R_φ	+ ρ₁(1-2p²V_S1²)R_ψ	= -(2ρ₂V_S2²pcosθ₂/V_P2)T_φ	+ ρ₂(1-2p²V_S2²)T_ψ
matrix notation case liquid-solid	AR = B ; R = [R_φ , R_ψ , T_φ , T_ψ]^t ; B = [-p , -cosθ₁/V_P1 , -ρ₁(1-2p²V_S1²) , 2ρ₁V_S1²pcosθ₁/V_P1]^t A_LR_L = B_L ; R_L = [R_φ , T_φ , T_ψ]^t ; B_L = [-cosθ₁/V_P1 , -ρ₁ , 0]^t
A[4,4] A_L[3,3]	p	cosη₁/V_S1	-p	cosη₂/V_S2
	-cosθ₁/V_P1	p	-cosθ₂/V_P2	-p
	ρ₁(1-2p²V_S1²) (ρ₁)	-2ρ₁V_S1²pcosη₁/V_S1	-ρ₂(1-2p²V_S2²)	-2ρ₂V_S2²pcosη₂/V_S2
	2ρ₁V_S1²pcosθ₁/V_P1 (0)	ρ₁(1-2p²V_S1²)	2ρ₂V_S2²pcosθ₂/V_P2	-ρ₂(1-2p²V_S2²)
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) = P-Q R_φ = (-P-Q)/D R_ψ =2S/D	AR = B ; R = [R_φ , R_ψ , T_φ , T_ψ]^t ; B = [-p , -cosθ₁/V_P1 , -ρ₁(1-2p²V_S1²) , 2ρ₁V_S1²pcosθ₁/V_P1]^t = [-a₁₁ , a₂₁ , -a₃₁ , a₄₁]^t ; D = (a₁₁A₁₁+a₃₁A₃₁)-(a₂₁A₂₁+a₄₁A₄₁) ; DR_φ = -(a₁₁A₁₁+a₃₁A₃₁)-(a₂₁A₂₁+a₄₁A₄₁) DR_ψ/2 = a₁₁a₂₁(a₃₃a₄₄-a₄₃a₃₄) +a₁₁a₄₁(a₂₃a₃₄-a₃₃a₂₄) +a₃₁a₄₁(a₁₃a₂₄-a₂₃a₁₄) DT_φ/2 = -(a₁₁a₂₁(a₃₂a₄₄-a₄₂a₃₄) +a₁₁a₄₁(a₂₂a₃₄-a₃₂a₂₄) +a₃₁a₄₁(a₁₂a₂₄-a₂₂a₁₄)) DT_ψ/2 = (a₁₁a₂₁(a₃₂a₄₃-a₄₂a₃₃) +a₁₁a₄₁(a₂₂a₃₃-a₃₂a₂₃) +a₃₁a₄₁(a₁₂a₂₃-a₂₂a₁₃)
P incident upgoing k_z^I = -ωcosθ₂/V_P2	Ar = b ; r = [t^φ,t^ψ,r^φ,r^ψ]^t b = [p , -cosθ₂/V_P2 , ρ₂(1-2p²V_S2²) , 2ρ₂V_S2²pcosθ₂/V_P2]^t = [-a₁₃ , a₂₃ , -a₃₃ , a₄₃]^t D = (a₁₃A₁₃+a₃₃A₃₃)-(a₂₃A₂₃+a₄₃A₄₃) ; Dr_φ = -(a₁₃A₁₃+a₃₃A₃₃)-(a₂₃A₂₃+a₄₃A₄₃) Dr_ψ/2 = a₁₃a₂₃(a₃₁a₄₂-a₄₁a₃₂) +a₁₃a₄₃(a₂₁a₃₂-a₃₁a₂₂) +a₃₃a₄₃(a₁₁a₂₂-a₂₁a₁₂) Dt_φ/2 = -(a₁₃a₂₃(a₃₄a₄₂-a₄₄a₃₂) +a₁₃a₄₃(a₂₄a₃₂-a₃₄a₂₂) +a₃₃a₄₃(a₁₄a₂₂-a₂₄a₁₂)) Dt_ψ/2 = a₁₃a₂₃(a₃₄a₄₁-a₄₄a₃₁) +a₁₃a₄₃(a₂₄a₃₁-a₃₄a₂₁) +a₃₃a₄₃(a₁₄a₂₁-a₂₄a₁₁)
case liquid-solid reflpsl.m	A_LR_L = B_L ; R_L = [R_φL , T_φL , T_ψL]^t ; B_L = [-cosθ₁/V_P1 , -ρ₁ , 0]^t A_Lr_L = b_L ; r_L = [t_φL , r_φL , r_ψL]^t ; b_L = [-cosθ₂/V_P2 , ρ₂(1-2p²V_S2²) , 2ρ₂V_S2²pcosθ₂/V_P2]^t D_L = -ρ₂²cosθ₁/V_P1((1-2p²V_S2²)²+4V_S2⁴p²cosθ₂/V_P2cosη₂/V_S2) -ρ₁ρ₂cosθ₂/V_P2 D_LR_φL = -ρ₂²cosθ₁/V_P1((1-2p²V_S2²)²+4V_S2⁴p²cosθ₂/V_P2cosη₂/V_S2) +ρ₁ρ₂cosθ₂/V_P2 D_LT_φL = -2ρ₁ρ₂cosθ₁/V_P1(1-2p²V_S2²) ; D_Lt_φL = -2ρ₂²cosθ₂/V_P2(1-2p²V_S2²) D_LT_ψL = -4ρ₁ρ₂cosθ₁/V_P1V_S2²pcosθ₂/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η₁/V_S1 D = det(A) = -P'+Q' R'_ψ = (-P'-Q')/D	AR'=B' ; R' = [R'_φ , R'_ψ , T'_φ , T'_ψ]^t ; B' = [cosη₁/V_S1 , -p , -2ρ₁V_S1²pcosη₁/V_S1 , -ρ₁(1-2p²V_S1²)]^t = [a₁₂ , -a₂₂ , a₃₂ , -a₄₂]^t ; D = -(a₁₂A₁₂+a₃₂A₃₂)+(a₂₂A₂₂+a₄₂A₄₂) ; DR'_ψ = -(a₁₂A₁₂+a₃₂A₃₂)-(a₂₂A₂₂+a₄₂A₄₂)
SV incident upgoing k_z^I = -ωcosη₂/V_S2	Ar'=b' ; r' = [t'_φ , t'_ψ , r'_φ , r'_ψ]^t b' = [cosη₂/V_S2 , p , -2ρ₂V_S2²pcosη₂/V_S2 , ρ₂(1-2p²V_S2²)]^t = [a₁₄ , -a₂₄ , a₃₄ , -a₄₄]^t ; D = -(a₁₄A₁₄+a₃₄A₃₄)+(a₂₄A₂₄+a₄₄A₄₄) Dr'_ψ = -(a₁₄A₁₄+a₃₄A₃₄)-(a₂₄A₂₄+a₄₄A₄₄)
total reflection of downgoing SV wave p = sinη₁/V_S₁ 1/V_S1 > p > 1/V_P1 p > 1/V_S2 > 1/V_P2 R'_ψ = e^iχ'(p)	imaginary terms of A : a₁₄ = i(p²V_S2²-1)^½/V_S2 ; a₂₁ = -i(p²V_P1²-1)^½/V_P1 ; a₂₃ = -i(p²V_P2²-1)^½/V_P2 ; a₃₄ = -2iρ₂V_S2²p(p²V_S2²-1)^½/V_S2 ; a₄₁ = 2iρ₁V_S1²p(p²V_P1²-1)^½/V_P1 ; a₄₃ = 2iρ₂V_S2²p(p²V_P2²-1)^½/V_P2 ; ⇒ A₁₂, A₃₂ and P' are imaginary, A₁₂, A₃₂ and Q' are real R'_ψ = (-P'-Q')/(-P'+Q') = e^iχ'(p) χ'(p) = 2Arctan(P'/iQ') = 2Arctan((1/i)(a₁₂A₁₂+a₃₂A₃₂)/(a₂₂A₂₂+a₄₂A₄₂))

weak contrasts	ρ₂ = ρ+δρ ; V_P₂ = V_P+δV_P ; θ₂ = θ+δθ ; δθ = tanθδV_P/V_P ; V_S₂ = V_S+δV_S ; η₂ = η + δη ; δη = tanηδV_S/V_S a₁₁ = -a₁₃ = a₂₂ = -a₂₄ = p ; a₁₂ = cosη/V_S ; a₂₁ = -cosθ/V_P ; a₃₁ = a₄₂ = ρ(1-2p²V_S²) ; a₃₂ = -2ρV_S²pcosη/V_S ; a₄₁ = 2ρV_S²pcosθ/V_P a₁₄ = (a₁₂+δa₁₂) ; a₂₃ = (a₂₁+δa₂₁) ; a₃₃ = a₄₄ = -(a₃₁+δa₃₁) ; a₃₄ = (a₃₂+δa₃₂) ; a₄₃ = (a₄₁+δa₄₁)
D = det(A)	D = (p²-a₂₁a₁₂)(a₃₁²-a₄₁a₃₂) - (2pa₂₁)(-2a₃₁a₃₂) + (-p²-a₂₁a₁₂)(a₃₂a₄₁+a₃₁²) + (a₁₂a₂₁+p²)(-a₃₁²-a₄₁a₃₂) - (-2pa₁₂)(2a₃₁a₄₁) + (p²-a₂₁a₁₂)(a₃₁²-a₄₁a₃₂) D = 2(p²-a₂₁a₁₂)(a₃₁²-a₄₁a₃₂) -2(p²+a₂₁a₁₂)(a₃₂a₄₁+a₃₁²) +4pa₃₁(a₂₁a₃₂+a₁₂a₄₁) = -4p²a₄₁a₃₂ -4a₂₁a₁₂a₃₁² +4pa₃₁(a₂₁a₃₂+a₁₂a₄₁) D = 4cosθ/V_Pcosη/V_S [p²(2ρV_S²p)² +(ρ(1-2p²V_S²))² +4ρ²p²V_S²(1-2p²V_S²)] = 4ρ²cosθ/V_Pcosη/V_S
R_φ	DR_φ = (-p²-a₂₁a₁₂)((a₃₁+δa₃₁)²-(a₄₁+δa₄₁)(a₃₂+δa₃₂)) - (-p(a₂₁+δa₂₁)+pa₂₁)(-a₃₂(a₃₁+δa₃₁)-a₃₁(a₃₂+δa₃₂)) + (p²-a₂₁(a₁₂+δa₁₂))(a₃₂(a₄₁+δa₄₁)+a₃₁(a₃₁+δa₃₁)) + (a₁₂(a₂₁+δa₂₁)+p²)(a₃₁(a₃₁+δa₃₁)-a₄₁(a₃₂+δa₃₂)) - (-pa₁₂-p(a₁₂+δa₁₂))(-a₃₁(a₄₁+δa₄₁)+a₄₁(a₃₁+δa₃₁)) + (p²-(a₂₁+δa₂₁)(a₁₂+δa₁₂))(-a₃₁²-a₄₁a₃₂)) DR_φ = (-p²-a₂₁a₁₂)(2a₃₁δa₃₁-a₄₁δa₃₂-a₃₂δa₄₁) - 2pa₃₂a₃₁δa₂₁ + (p²-a₂₁a₁₂)(a₃₂δa₄₁+a₃₁δa₃₁)-a₂₁δa₁₂(a₃₂a₄₁+a₃₁²) + (p²+a₁₂a₂₁)(a₃₁δa₃₁-a₄₁δa₃₂)+a₁₂δa₂₁(a₃₁²-a₄₁a₃₂) + 2pa₁₂(-a₃₁δa₄₁+a₄₁δa₃₁) + (a₂₁δa₁₂+a₁₂δa₂₁)(a₃₁²+a₄₁a₃₂) DR_φ = 2(p(pa₃₂-a₁₂a₃₁)δa₄₁ + a₁₂(pa₄₁-a₂₁a₃₁)δa₃₁ + a₃₁(a₃₁a₁₂-pa₃₂)δa₂₁) DR_φ = 2ρ(-2p²cosη/V_Sδ(ρV_S²cosθ/V_P) +cosη/V_Scosθ/V_Pδ(ρ(1-2p²V_S²) + ρ(1-2p²V_S²)cosη/V_Sδ(-cosθ/V_P)) DR_φ = 2ρcosη/V_S(cosθ/V_P((1-4p²V_S²)δρ - 4ρp²δ(V_S²)) - ρδ(cosθ/V_P)) R_φ = (1/2)((1-4p²V_S²)δρ/ρ - 8p²V_S²δV_S/V_S - δ(cosθ/V_P)/(cosθ/V_P)) R_φ = (1/2)((1-4p²V_S²)δρ/ρ - 8p²V_S²δV_S/V_S + (1/cosθ)²δV_P/V_P) R_φ = (1/2)((δρ/ρ+δV_P/V_P) + (δV_P/V_P -4V_S²/V_P²(δρ/ρ+2δV_S/V_S)sin²θ + (tg²θ-sin²θ)δV_P/V_P) ≈ R₀+Gsin²θ
R_ψ	DR_ψ = (2pa₂₁)((a₃₁+δa₃₁)²-(a₄₁+δa₄₁)(a₃₂+δa₃₂)) - p((a₂₁+δa₂₁)+a₂₁)(a₃₁(a₃₁+δa₃₁)-a₄₁(a₃₂+δa₃₂)) + (-p²-a₂₁(a₁₂+δa₁₂))(-a₃₁(a₄₁+δa₄₁)+a₄₁(a₃₁+δa₃₁)) + p(-(a₂₁+δa₂₁)+a₂₁)(-a₃₁(a₃₁+δa₃₁)-a₄₁(a₃₂+δa₃₂)) - (p²-a₂₁(a₁₂+δa₁₂))(a₃₁(a₄₁+δa₄₁)+a₄₁(a₃₁+δa₃₁)) + (p²-(a₂₁+δa₂₁)(a₁₂+δa₁₂))(2a₃₁a₄₁) DR_ψ = 2pa₂₁(2a₃₁δa₃₁-a₄₁δa₃₂-a₃₂δa₄₁) - 2pa₂₁(a₃₁δa₃₁-a₄₁δa₃₂)-pδa₂₁(a₃₁²-a₄₁a₃₂) + (-p²-a₂₁a₁₂)(-a₃₁δa₄₁)+a₄₁δa₃₁)) + (-pδa₂₁)(-a₃₁²-a₄₁a₃₂) - (p²-a₂₁a₁₂)(a₃₁δa₄₁+a₄₁δa₃₁) +a₂₁δa₁₂2a₃₁a₄₁ + (-2a₃₁a₄₁(a₂₁δa₁₂+a₁₂δa₂₁)) DR_ψ = 2(p(a₂₁a₃₁-pa₄₁)δa₃₁ + a₂₁(-pa₃₂+a₁₂a₃₁)δa₄₁ + a₄₁(pa₃₂-a₃₁a₁₂)δa₂₁) DR_ψ = -2ρpcosθ/V_P(δ(ρ(1-2p²V_S²)) + cosη/V_Sδ(2ρV_S²cosθ/V_P) + 2ρV_S²cosη/V_Sδ(-cosθ/V_P)) DR_ψ = -2ρpcosθ/V_P ((1-2p²V_S²+2V_S²cosθ/V_Pcosη/V_S)δρ + 2ρ(-p²+cosη/V_Scosθ/V_P)δ(V_S²) R_ψ = -(1/2)pV_S/cosη ((1-2p²V_S²+2V_S²cosθ/V_Pcosη/V_S)δρ/ρ + 4V_S²(-p²+cosη/V_Scosθ/V_P)δV_S/V_S) R_ψ ≈ -(1/2)sinη ((1+2V_S/V_P)δρ/ρ + 4V_S/V_PδV_S/V_S) pour incidence faible

1/V_R = p = k_x/ω > 1/V_S > 1/V_P	evanescent P wave k_z^P = iω(p²V_P²-1)^½/V_P	evanescent S wave k_z^S = iω(p²V_S²-1)^½/V_S
σ_zz = σ_zx = 0 on free surface z = 0	(1-2p²V_S²)Φ+i2pV_S(p²V_S²-1)^½Ψ = 0 -i2V_S²p(p²V_P²-1)^½/V_PΦ+(1-2p²V_S²)Ψ = 0
	déterminant = 0 (1-2p²V_S²)²-4p²V_S⁴(p²V_S²-1)^½/V_S(p²V_P²-1)^½/V_P = 0 (2-V_R²/V_S²)²-4(1-V_R²/V_P²)^½(1-V_R²/V_S²)^½ = 0
(u_x , u_z) = (U_x(z) , U_z(z))e^iω(x/V_R-t)	U_x = iωpΦe^{-ωz(p²V_P²-1)^½/V_P}+ω(p²V_S²-1)^½/V_SΨe^{-ωz(p²V_S²-1)^½/V_S} U_z = -ω(p²V_P²-1)^½/V_PΦe^{-ωz(p²V_P²-1)^½/V_P}+iωpΨe^{-ωz(p²V_S²-1)^½/V_S} U_x = e^iπ/2ωpΦ(e^{-ωz(p²V_P²-1)^½/V_P}+(1/(2p²V_S²)-1)e^{-ωz(p²V_S²-1)^½/V_S}) U_z = e^iπω(p²V_P²-1)^½/V_PΦ(e^{-ωz(p²V_P²-1)^½/V_P}+(1/(2p²V_S²)-1)^-1e^{-ωz(p²V_S²-1)^½/V_S})

	standing wave in z k_z for normal modes	travelling wave in x , k_x = (ω²/V_S²-k_z²)^½ phase velocity V(ω) = V_ax = ω/k_x	ω ≥ high-frequency cut-off ω_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^ik_zz+e^-ik_zz = 2cos(k_zz) Hk_{z n} = nπ , nλ_{z n}/2 = H , n = 1,2,3...	k_{x n} = (ω²/V_S²-(nπ/H)²)^½ V_n(ω) = ω/k_{x n}	ω_{c n} = nπV_S/H
layer (-H ≤ z ≤ 0) over half-space (z≥ 0) R_surf = 1 ; R_interf = e^iχ(p)	e^ik_1zz+e^iχ(p)e^-ik_1zz = 2e^iχ(p)/2cos(k_1zz-χ(p)/2) Hk_{1z n}+χ(p)/2 = nπ , n = 0,1,2,3...	k_{x n} = (ω²/V_S1²-(nπ-χ(p)/2)²/H²)^½ 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²-1/V_S2²)^½)
displov.m	k_1z = (ω/V_S1)(1-p²V_S1²)^½ χ(p) = -2Arctan(ρ₂V_S2(p²V_S2²-1)^½/ρ₁V_S1(1-p²V_S1²)^½) = -2Arctan(I₂/I₁) dispersion relation of Love wave : ω(p_n) = (-χ(p_n)/2+nπ)V_S1/(H(1-p_n²V_S1²)^½) = (Arctan(I₂/I₁)+nπ)V_S1/(H(1-p_n²V_S1²)^½)

boundary conditions	normal modes in z	travelling wave in x	cut-off frequency
R_surf = e^iχ'₁(p) R_interf = e^iχ'₂(p) χ'₁₂(p) = χ'₁(p)+χ'₂(p)	(n-χ'₁₂(p)/2π)λ_{1z n}/2 = H Hk_{1z n}+χ'₁₂(p)/2 = nπ	k_{x n} = (ω²/V_S1²-(nπ-χ'₁₂(p)/2)²/H²)^½ V_{R n}(ω) = ω/k_{x n} = 1/p_n	ω_{c n} = (nπ-χ'₁₂(1/V_S2))/2)/(H(1/V_S1²-1/V_S2²)^½)
dispersion rel. of Rayleigh wave overtones n≥1 V_S1≤V_{R n}≤V_S2	χ'₁(p) = -2Arctan((4V_S1⁴p²(p²V_P1²-1)^½/V_P1(1-p²V_S1²)^½/V_S1)/(1-2p²V_S1²)²) = -2Arctan(S) χ'₂(p) = -2Arctan(Im(D)/Re(D)) = -2Arctan(I) tan(Hω(1-p²V_S1²)^½/V_S1) = -tan((χ'₁(p)+χ'₂(p))/2) = (I+S)/(1-IS) ω(p_n) = (-(χ'₁(p_n)+χ'₂(p_n))/2+nπ)V_S1/(H(1-p_n²V_S1²)^½) = (Arctan((I+S)/(1-IS))+nπ)V_S1/(H(1-p_n²V_S1²)^½)

boundary cond.	downgoing P1 , k_zP1 = ω(1-p²V_P1²)^½/V_P1 Z_P1 = e^-ik_zP1H	downgoing S1 , k_zS1 = ω(1-p²V_S1²)^½/V_S1 Z_S1 = e^-ik_zS1H	upgoing P1 , -k_zP1	upgoing S1 , -k_zS1	evanescent P2 , k_zP2 = iω(p²V_P2²-1)^½/V_P2	evanescent S2 , k_zS2 = iω(p²V_S2²-1)^½/V_S2
z = -H	σ_zz(k_zP1²)Z_P1	-2μ₁k_xk_zS1Z_S1	σ_zz(k_zP1²)/Z_P1	2μ₁k_xk_zS1/Z_S1	0	0
z = -H	-2μ₁k_xk_zP1Z_P1	σ_zx(k_zS1²)Z_S1	2μ₁k_xk_zP1/Z_P1	σ_zx(k_zS1²)/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μ₁k_xk_zS1	σ_zz(k_zP1²)	2μ₁k_xk_zS1	-σ_zz(k_zP2²)	2μ₂k_xk_zS2
	-2μ₁k_xk_zP1	σ_zx(k_zS1²)	2μ₁k_xk_zP1	σ_zx(k_zS1²)	2μ₂k_xk_zP2	-σ_zx(k_zS2²)

boundary cond.	waves guided by free surface in z = -H and prolongated in z = 0 by Z_P1 = e^{-ω(p²V_P1²-1)^½/V_P1H} (0) ; Z_S1 = e^{-ω(p²V_S1²-1)^½/V_S1H} (0)		waves guided by interf. in z = 0 and prolongated in z = -H by Z_P1 , Z_S1
z = -H	Rayleigh free surf.		(1-2p²V_S1²)Z_P1 (0)	-i2pV_S1(p²V_S1²-1)^½Z_S1 (0)	0	0
z = -H	Rayleigh free surf.		+i2V_S1²p(p²V_P1²-1)^½/V_P1Z_P1 (0)	+(1-2p²V_S1²)Z_S1 (0)	0	0
z = 0	pZ_P1 (0)	-i(p²V_S1²-1)^½/V_S1Z_S1 (0)	Stoneley at interface
	i(p²V_P1²-1)^½/V_P1Z_P1 (0)	pZ_S1 (0)
	ρ₁(1-2p²V_S1²)Z_P1 (0)	i2ρ₁V_S1²p(p²V_S1²-1)^½/V_S1Z_S1 (0)
	-i2ρ₁V_S1²p(p²V_P1²-1)^½/V_P1Z_P1 (0)	ρ₁(1-2p²V_S1²)Z_S1 (0)

z I_i = ρ_iV_Si²k_iz	downgoing wave e^k_izz k_iz = ω(1/V_Si²-p²)^½	upgoing wave e^-ik_izz	continuity of u_y and σ_yz at interface total reflection if free surface (upgoing = downgoing) : dispersion relation
-(h₃+h₂+h₁) layer 3	D₃e^-ik_3zh₃ = D₃(a₃-ib₃)	M₃e^ik_3zh₃ = M₃(a₃+ib₃)	A₃b₃+a₃B₃ = 0 tan(k_3zh₃)+(I₂/I₃)tan(k_2zh₂)+((I₁/I₃)-(I₁/I₂)tan(k_3zh₃)tan(k_2zh₂))tan(k_1zh₁+χ/2) = 0
-(h₂+h₁) layer 3	D₃ = A₃-iB₃	M₃ = A₃+iB₃	A₃ = A₂a₂-B₂b₂ = a₁a₂-(I₁/I₂)b₁b₂ B₃ = (I₂/I₃)(A₂b₂+B₂a₂) = (I₂/I₃)a₁b₂+(I₁/I₃)b₁a₂
-(h₂+h₁) layer 2	D₂e^-ik_2zh₂ = D₂(a₂-ib₂)	M₂e^ik_2zh₂ = M₂(a₂+ib₂)	A₂b₂+a₂B₂ = 0 tan(k_2zh₂)+(I₁/I₂)tan(k_1zh₁+χ/2) = 0
-h₁ layer 2	D₂ = A₂-iB₂	M₂ = A₂+iB₂	A₂ = a₁ B₂ = (I₁/I₂)b₁
-h₁ layer 1	D₁e^-ik_1zh₁ = a₁-ib₁	M₁e^ik_1zh₁ = a₁+ib₁	b₁ = 0 sin(k_1zh₁+χ/2) = 0
0 layer 1	D₁ = e^-iχ/2	M₁ = e^iχ/2	χ(p) = -2Arctan(ρV_S(p²V_S²-1)^½/ρ₁V_S1(1-p²V_S1²)^½) = -2Arctan(I/I₁)
> 0 1/2 space	evanescent wave e^-k_zz k_z = \|ω\|(p²-1/V²)^½

Bessel equation	solution	derivative
r²d²f/dr²+rdf/dr+ (k_r²r²-n²)f(r) = 0	f(r) = H_n(k_rr) = J_n(k_rr)+iY_n(k_rr) ≈ (2/πk_rr)^½e^{i(k_rr-π/4-nπ/2)} H_n(ik_rr) ≈ i^-(n+1)(2/πk_rr)^½e^-k_rr	dH_n(k_rr)/dr = k_rH_n-1(k_rr)-(n/r)H_n(k_rr) = (n/r)H_n(k_rr)-k_rH_n+1(k_rr)

Δφ(r,z) = ∂²φ/∂r²+(1/r)∂φ/∂r+∂²φ/∂z² = (1/V_P²)∂²φ/∂t²	φ(r,z) = f(r)e^i(k_zz-ωt)	k_p² = ω²/V_P²-k_z²	f(r) = AH₀(k_pr)
Δψ_θ(r,z) = (1/V_S²)∂²ψ_θ/∂t²	ψ_θ(r,z) = g(r)e^i(k_zz-ωt)	k_s² = ω²/V_S²-k_z²	g(r) = BH₁(k_sr)

Δφ(r,θ,z) = Δφ(r,z) + (1/r)²∂²φ/∂θ² = (1/V_P²)∂²φ/∂t²	φ(r,θ,z) =f(r)cos(nθ)e^i(k_zz-ωt)	k_p² = ω²/V_P²-k_z²	f(r) = AH_n(k_pr)
Δψ_r(r,θ,z) = (1/V_S²)∂²ψ_r/∂t²	ψ_r(r,z) = g(r)sin(nθ)e^i(k_zz-ωt)	k_s² = ω²/V_S²-k_z²	g(r) = BH_n+1(k_sr)
Δψ_θ(r,θ,z) = (1/V_S²)∂²ψ_θ/∂t²	ψ_θ(r,z) = -g(r)cos(nθ)e^i(k_zz-ωt)		g(r) = BH_n+1(k_sr)
Δψ_z(r,θ,z) = (1/V_S²)∂²ψ_z/∂t²	ψ_z(r,z) = h(r)sin(nθ)e^i(k_zz-ωt)		h(r) = CH_n(k_sr)

	u_r = cos(nθ)	u_θ = sin(nθ)	u_z = cos(nθ)	ε_rr = cos(nθ)	ε_rθ = sin(nθ)	ε_rz = cos(nθ)
Ae^i(k_zz-ωt)	n/rH_n(k_pr)-k_pH_n+1(k_pr)	-n/rH_n(k_pr)	ik_zH_n(k_pr)	(n(n-1)/r²-k_p²)H_n(k_pr)+k_p/rH_n+1(k_pr)	-n(n-0.5)/r²H_n(k_pr)+nk_p/rH_n+1(k_pr)	ik_z(n/rH_n(k_pr)-k_pH_n+1(k_pr))
Be^i(k_zz-ωt)	ik_zH_n+1(k_sr)	ik_zH_n+1(k_sr)	-k_sH_n(k_sr)	ik_z(k_sH_n(k_sr)-(n+1)/rH_n+1(k_sr))	ik_z(k_s/2H_n(k_sr)-(n+1)/rH_n+1(k_sr))	(-nk_s/rH_n(k_sr)+(k_s²-k_z²)H_n+1(k_sr))/2
Ce^i(k_zz-ωt)	n/rH_n(k_sr)	-n/rH_n(k_sr)+k_sH_n+1(k_sr)	0	n(n-1)/r²H_n(k_sr)-nk_s/rH_n+1(k_sr)	(-n(n-1)/r²+k_s²/2)H_n(k_sr)-k_s/rH_n+1(k_sr)	(ik_zn/2r)H_n(k_sr)