Past exam of the mathematics course of the University of Cambridge / 2003 / iii / Paper 78 (Seismic Waves)

At normal incidence the P wave and the two S wave polarizations decouple. In each uniform layer the relevant displacement satisfies $\rho u_{tt}=C{u}_{zz}$, with $C=\lambda +2\mu$ for a longitudinal wave and $C=\mu$ for a shear wave. Its speed is $v=\sqrt{C/\rho }$ and its seismic impedance is $Z=\rho v$. A welded interface preserves particle velocity and traction. For displacement or velocity amplitudes, incidence from medium 0 onto medium 1 therefore givesReflection reverses the displacement sign when the receiving seismic impedance is larger. The energy-flux fractions areThe impedance factor matters: a displacement transmission coefficient larger than one does not imply energy amplification. Equal impedance eliminates reflection even when the densities and speeds differ.

For a stack, every interface creates reflected waves which can return to all earlier interfaces. A seismic layer transfer matrix incorporates those multiple reflections exactly. Use harmonic time dependence $e^{-i\omega t}$, particle velocity $V$ and compressive traction $P=-{\sigma }_{zz}$ for a P-wave, or the negative relevant shear traction for an S-wave. Through layer $j$ of thickness $h_j$,This follows by solving $V_z=i\omega P/C_j$ and $P_z=i\omega {\rho }_{j}V$. Continuity means that the state passes unchanged across each welded interface. Put $M=M_N\cdots M_1=\left(\begin{array}{cc}a& b\\ c& d\end{array}\right)$, with $ad-bc=1$, and let the exterior impedances be $Z_L,Z_R$. A unit incident velocity gives top state $(1+r,Z_L(1-r))$ and bottom state $(t,Z_{R}t)$. Eliminating these states yieldsFor real moduli and no attenuation, conservation of elastic-wave energy flux gives $|r{|}^2+(Z_R/Z_L)|t{|}^2=1$. Reciprocal incidence gives equal energy transmission in the reverse direction; flux-normalized amplitudes express that reciprocity without unequal exterior impedances obscuring it.

The wave phase accumulated in each layer determines constructive or destructive interference. For one layer between identical exterior media, an integral number of half-wavelengths makes $M_j=\pm I$ and gives perfect transmission, regardless of the layer impedance. At quarter-wave phase thickness, a single matching layer with $Z_j=\sqrt{Z_L{Z}_R}$ also gives zero reflection: the layer transforms the terminal impedance to $Z_j^2/Z_R=Z_L$. Away from these matching conditions there are transmission peaks, reflection maxima and rapid frequency-dependent phase changes. In the time domain these are a direct transmitted pulse followed by delayed internal reverberations and reflected echoes. Long residence in a resonant stack produces narrow peaks and long ringing, without creating energy.

For a periodic two-layer stack, the half-trace of one cell isThis is the periodic seismic multilayer stop band criterion. When the absolute right side exceeds one, the cell eigenvalues are a reciprocal growing/decaying pair and a long finite stack has exponentially small transmission. The missing transmitted energy is reflected, not absorbed. Large impedance contrast and repeated quarter-wave layers strengthen such stop bands. In pass bands the cell modes propagate, and a finite stack superposes transmission resonances on that propagation. An irregular stack need not have equally spaced peaks: the individual travel times and impedance contrasts set its interference pattern.

At sufficiently low frequency a finite stack has $M\to I$, so its leading reflection is that of the two exterior media alone. For many thin repeated layers with wavelength much longer than a cell, their effective density is the volume-weighted density and their normal modulus is the harmonic mean, $C_{\mathrm{eff}}^{-1}=⟨C^{-1}⟩$. This follows by expanding each propagation matrix to first order and averaging its two coefficients. At high frequency, sharp interfaces still reflect: there is no general disappearance of reflection, only increasingly rapid interference fringes. Smooth impedance grading instead suppresses short-wavelength reflection when its variation scale is long compared with the wavelength.

These calculations provide a reproducible investigation: vary the layer travel times, impedances and number of cells; compute the matrix product; plot reflected/transmitted fluxes and check their sum. Adding attenuation makes the flux sum smaller than one and reduces reverberation peaks, while oblique incidence introduces Snell law for elastic waves and P-SV mode conversion absent from this one-dimensional setting. Impedance contrasts generate reflections, propagation phases organize their interference, and energy conservation separates reflection from genuine absorption.

Let $u(x,z,t)$ be displacement along the source line and let $F(t)$ be its force per unit length. The shear-horizontal wave equation isRotational symmetry about the line makes $u=u(r,t)$. Integrate over a disk of radius $r$ and use the divergence theorem:At a fixed regular time, a logarithmic singularity and its time derivatives are locally integrable; the inertial integral is $O(r^2|\log r|)$ and vanishes as $r\downarrow 0$. The leading balance is consequently $r{u}_r=-F(t)/(2\pi \mu )$, giving the SH line-source near-field singularityThe reference length $r_{\ast }$ only changes the bounded term; boundary conditions and source history determine that term. This is a local fixed-time statement, interpreted distributionally in time if the source force has jumps. It establishes logarithmic spatial dependence directly, without extrapolating far-field ray spreading into the source.

Away from the line source, insert a leading wavefront form $u\sim a(\mathbf{x})W(t-T(\mathbf{x}))$ into the shear-horizontal wave equation and match the two most singular orders. The first gives the eikonal equation $|\nabla T|=1/\beta$. The next givesEquivalently, SH ray-tube amplitude transport conserves $\rho \beta a^{2}J$, where $J$ is the transverse width of a two-dimensional ray tube. A line source has straight radial elastic rays and $J=rd\theta$ for fixed angular width, soThe same geometrical spreading applies to the first nonzero derivative discontinuity or leading wavefront singularity; the source history determines which waveform $W$ and normalization $C$ occur. In three dimensions a point source instead has $J\propto r^2$ and amplitude proportional to $r^{-1}$.

This does not contradict the SH line-source near-field singularity. Ray transport describes an arrival-front limit with $t-r/\beta$ small, or the high-frequency field outside the immediate source neighborhood. The logarithm describes $r\downarrow 0$ at fixed elapsed time and includes the tail behind the two-dimensional wavefront. For example, a causal step force $F_{0}H(t)$ has the exact radial displacementFor $r\ll \beta t$, this is $F_0\log (2\beta t/r)/(2\pi \mu )$, the near-source logarithm. For $t=r/\beta +\tau$ with $0<\tau \ll r/\beta$, it isThe coefficient of the leading arrival singularity therefore has exactly the ray-theory dependence. A force impulse gives the corresponding $r^{-1/2}H(\tau )/\sqrt{\tau }$ leading field. The limits describe different regions of one causal solution; neither is a uniform approximation to the other.

Choose interface $z=0$, source $(0,-h)$ in medium 1 and transmitted medium 2 at $z>0$. Write ${\mu }_j={\rho }_j{\beta }_j^2$, take $p$ to be signed tangential slowness, and defineThe propagating range is $|p|<1/\max ({\beta }_1,{\beta }_2)$; the absolute value is needed if both azimuthal directions are included. Snell law for elastic waves preserves $p$. The incident ray reaches the interface after time $t_1$ at horizontal position $x_I$, withIn the transmitted medium its velocity is ${\beta }_2(\sin {\theta }_2,\cos {\theta }_2)=({\beta }_2^{2}p,{\beta }_2^2{q}_2)$. Therefore the refraction of a cylindrical SH wavefront givesTime is measured from source emission. For a later emission add its time to the clock. These straight transmitted rays and their common-time endpoints determine the refracted wavefront; critical incidence and an evanescent transmitted field are excluded by the strict slowness bound.

At a welded interface, elastic-interface boundary conditions require continuous SH displacement and shear traction. Incident and reflected vertical slownesses have signs $q_1,-q_1$, givingThus the local displacement coefficient isIt multiplies the leading waveform discontinuity as well as a harmonic displacement. At the interface the incident amplitude from SH ray-tube amplitude transport is $C/\sqrt{r_1}$, where $r_1=h/({\beta }_1{q}_1)$ is the incident path length and $C$ is the source normalization from the previous part.

The transmitted geometrical spreading is not simply its total path length. Differentiate the ray position with respect to $p$ at fixed travel time and project onto the unit direction transverse to the transmitted ray, ${\mathbf{m}}_2=({\beta }_2{q}_2,-{\beta }_{2}p)$. Since $d{x}_I/dp=h/({\beta }_1^2{q}_1^3)$, the transverse width per unit $p$ isThe derivative of $t_1$ contributes only in the ray direction and drops out of this projection. At entry,Within uniform medium 2 the leading amplitude varies as $J_2^{-1/2}$, so the required refraction of a cylindrical SH wavefront amplitude isThe local transmission coefficient and ray-tube spreading together incorporate both impedance contrast and refraction. Indeed the normal energy transmission fraction is ${\mu }_2{q}_2|T_{\mathrm{SH}}{|}^2/({\mu }_1{q}_1)$, while the change of transverse width supplies the refraction factor in the conserved tube flux. If the media coincide, $q_1=q_2$ and $T_{\mathrm{SH}}=1$, reducing this expression to $C/\sqrt{{\beta }_{1}t}$ as required. The formula holds for the stated propagating slownesses, without a critical or head-wave contribution.

Let the layer occupy $0<z<h$, with a free surface at $z=0$ and rigid base at $z=h$. Write an SH-wave as $u_y=Y(z)e^{i(kx-\omega t)}$. The equation and boundary conditions giveThus $Y=A\cos (qz)$ and the guided SH modes between rigid and free planes haveThe propagating branch has $\omega >{\omega }_n$ and $k=\sqrt{{\omega }^2-{\omega }_n^2}/{\beta }^{\prime }$. Hence the rigid-base approximation to Love waves predictsAt each threshold $c_n\to \infty$ and $U_n\to 0$. As frequency rises, the phase velocity decreases towards ${\beta }^{\prime }$ and the group velocity increases towards ${\beta }^{\prime }$. Frequencies below ${\omega }_n$ have imaginary horizontal wavenumber and do not give propagating guided modes. The infinite threshold phase speed does not violate causality: the limiting group speed is zero. These quarter-wave thresholds differ from the actual Love-wave cutoff frequencies of an elastic substrate.

The full Love wave problem replaces the rigid base by a decaying substrate displacement. With $Y=A\cos (qz)$ in the layer and $Y=B{e}^{-\kappa (z-h)}$ below it,where ${\mu }^{\prime }={\rho }^{\prime }{\beta }^{\prime 2}$ and $\mu =\rho {\beta }^2$. Displacement continuity gives $B=A\cos (qh)$ and traction continuity gives ${\mu }^{\prime }q\sin (qh)=\mu \kappa \cos (qh)$. The rigid base is therefore a good approximation only when $\mu \kappa \gg {\mu }^{\prime }q$: then the substrate displacement is small and $qh$ lies close to ${\chi }_n=(n+1/2)\pi$.

Put $ϵ={\beta }^{\prime }/\beta$, $\eta ={\mu }^{\prime }/\mu$, $y=qh$, $W=\omega h/{\beta }^{\prime }$ and $K=kh$. The complete dispersion branch can be parameterized without root finding:Then $c={\beta }^{\prime }W/K$ and $U={\beta }^{\prime }{W}^{\prime }/K^{\prime }$. At the lower end, $c\to \beta$ and $U\to \beta$. The Love-wave cutoff frequencies are ${\omega }_{c,n}=n\pi {\beta }^{\prime }/[h\sqrt{1-{ϵ}^2}]$, with zero cutoff for the fundamental. Thus the rigid-layer formula cannot describe the fundamental's low-frequency limit or the low-frequency part of an overtone: it would either give no propagating mode or predict $c>\beta$, outside the trapping interval ${\beta }^{\prime }<c<\beta$.

For comparable densities and a fixed mode index, $\eta =({\rho }^{\prime }/\rho ){ϵ}^2\ll 1$. Write $y={\chi }_n-\delta$ near the strong-dispersion transition. In the region controlling the group-velocity minimum of a high-contrast Love wave, $\delta \ll 1$ and $\eta /\delta \gg ϵ$. The leading derivatives giveTheir minimum is atThe frequency estimate follows by inserting ${\delta }_{\min }$ into $W\simeq ({\chi }_n-\delta )[1+{\eta }^2/(2{\delta }^2)+{ϵ}^2/2]$: the two leading terms of size $\delta$ cancel. The real group velocity has a finite minimum and turns upwards near the rigid-layer threshold, whereas the rigid formula has zero group speed at its artificial threshold.

For a useful estimate of the validity range, put $\Delta =W-{\chi }_n>0$. The rigid prediction has $K_R^2=W^2-{\chi }_n^2\simeq 2{\chi }_n\Delta$. The finite-load phase correction is $\delta \simeq \eta {\chi }_n/(\kappa h)$; its changes to both $K$ and $dK/dW$ are small ifThe first condition keeps $K_R$ away from the substrate propagation boundary; the second controls the rapidly changing interface phase. Equivalently, the rigid formulas are reliable sufficiently above the actual group-velocity minimum, beyond its narrow transition interval. At fixed $\omega /{\omega }_n>1$ separated from one, they become asymptotically accurate as the contrast grows. They are particularly simple at $\omega \gg {\omega }_n$, where both velocities approach ${\beta }^{\prime }$. Within the transition and below it the full dispersion relation is required.

The small speed ratio alone is insufficient if the layer density grows so strongly that ${\mu }^{\prime }/\mu$ is not small. For such materials the exact load criterion $\mu \kappa \gg {\mu }^{\prime }q$ and its frequency derivatives should be used; the preceding contrast estimates assume a bounded order-one density ratio and fixed $n$.

With vertical shear modulus ${\rho }^{\prime }{\beta }_z^{\prime 2}$ and horizontal shear modulus ${\rho }^{\prime }{\beta }_x^{\prime 2}$, force balance therefore gives the anisotropic shear-horizontal wave equationFor $u_y=Y(z)e^{i(kx-\omega t)}$, the layer vertical wavenumber and substrate decay exponent areThe free top still gives $Y=A\cos (qz)$; the substrate still has $Y=B{e}^{-\kappa (z-h)}$. Welded displacement continuity and the now anisotropic traction continuity giveHence the anisotropic Love-wave dispersion relation isEquivalently, writing $c=\omega /k$,This replaces the isotropic speed in the layer equation by separate horizontal and vertical speeds and replaces the layer traction modulus by ${\rho }^{\prime }{\beta }_z^{\prime 2}$. No additional SH/P-SV coupling arises for the aligned symmetry axes and constitutive law here.

A rigid base still imposes $q_n={\chi }_n/h$, but nowThusThe vertical speed sets the resonance frequencies, while the horizontal speed sets the high-frequency limits and $c_n{U}_n={\beta }_x^{\prime 2}$.

For the elastic substrate, put ${ϵ}_x={\beta }_x^{\prime }/\beta$, $\zeta ={\rho }^{\prime }{\beta }_x^{\prime }{\beta }_z^{\prime }/(\rho {\beta }^2)$, $W=\omega h/{\beta }_z^{\prime }$ and $K=kh{\beta }_x^{\prime }/{\beta }_z^{\prime }$. The exact parameterization isIt reproduces the isotropic analysis with ${\beta }^{\prime }$ replaced by ${\beta }_z^{\prime }$ in the frequency scale, by ${\beta }_x^{\prime }$ in the velocities, and by the two effective contrast parameters above. In particular the true cutoffs are ${\omega }_{c,n}=n\pi {\beta }_z^{\prime }/[h\sqrt{1-{ϵ}_x^2}]$, the group minima occur near ${\omega }_{n,z}$, and a sufficient rigid-limit validity range isFor comparable horizontal and vertical layer speeds, a bounded density ratio and fixed mode index, the minimum estimate becomes ${\delta }_{\min }\sim ({\zeta }^2{\chi }_n/2{)}^{1/3}$ and $U_{\min }/{\beta }_x^{\prime }\sim 3(\zeta /{\chi }_n{)}^{1/3}/2^{2/3}$. If their ratio is extreme, that more specific minimum-depth estimate additionally needs $\zeta /{\delta }_{\min }\gg {ϵ}_x$; the exact parameterization and load criterion remain valid without it. As in the isotropic case, neither a slow horizontal nor a slow vertical speed alone justifies the rigid boundary at cutoff or within the strong-dispersion transition.

Let $d(\mathbf{y},t)=[u^{[n]}]$ denote the signed opening, positive when the faces separate, and let $S$ be the crack surface. Its displacement jump is $\mathbf{b}=d\mathbf{n}$. The singular strain associated with this jump has moment density $c_{ijkl}{b}_k{n}_l$. For the isotropic elastic stiffness tensor, contraction givesThe equivalent distributional body force is the negative divergence of this prescribed moment density:The resulting forced elastic wave in an isotropic solid equation is $\rho u_{i,tt}={\partial }_j{\sigma }_{ij}+f_i$. This is a tensile-crack seismic source; it need not be represented as an ordinary smooth volume force. For a small nearly planar patch of area $S_0$, uniform normal $\mathbf{n}$ and opening $d(t)$, define its volume change $V(t)=S_{0}d(t)$. Its seismic moment tensor and point-source force areThe sign convention makes a positive opening a positive extensional moment. Reversing the convention for the displacement jump changes both signs together; the force is independent of reversing the orientation of the crack normal.

Complete the crack normal to an orthonormal basis $({\mathbf{e}}_1,{\mathbf{e}}_2,\mathbf{n})$ with the first two vectors in the crack plane. The seismic moment tensor is diagonal in this basis:Each term gives an elastic force dipole, $f_i=-D{e}_i{e}_j{\partial }_j\delta$: two equal and opposite forces in direction $\mathbf{e}$, separated along $\mathbf{e}$, with moment $D$. Hence there are two identical in-plane dipoles of strength $\lambda V$ and one normal dipole of strength $(\lambda +2\mu )V$. For positive opening and ordinary positive Lamé parameters, all are extensional and the normal one is larger by $2\mu V$.

More generally, stable isotropic linear elasticity requires $\mu >0$ and $3\lambda +2\mu >0$. These imply $\lambda +2\mu >|\lambda |$, so the normal dipole has greater magnitude even if an auxetic material has negative $\lambda$ and compressive in-plane dipoles. The point source has zero net force, because it differentiates a delta function, and zero net torque, because its moment tensor is symmetric. Its normal moment distinguishes opening from a pure shear dislocation.

Treat a small patch as the point tensile-crack seismic source at ${\mathbf{y}}_0$, with $M(t)=KV(t)$ and $K=\lambda I+2\mu \mathbf{n}{\mathbf{n}}^{\mathsf{T}}$. Let $\mathbf{e}=(\mathbf{x}-{\mathbf{y}}_0)/R$, $\gamma =\mathbf{e}\cdot \mathbf{n}$, $\mathbf{a}=K\mathbf{e}$, $m=\mathbf{e}\cdot K\mathbf{e}$ and $s=\mathrm{tr}K$. Source convolution and an integration by parts give the point-moment elastodynamic displacementThe derivatives act on angular factors and on retarded times as well as on powers of $R$; keeping only the explicit $1/R$ Green-function terms would miss the static displacement.

Set $V_P=V(t-R/\alpha )$, $V_S=V(t-R/\beta )$ andConvolving the elastodynamic Green tensor with the opening history givesUse ${\partial }_{j}R=e_j$, ${\partial }_j{e}_i=({\delta }_{ij}-e_i{e}_j)/R$, ${\partial }_j{V}_P=-e_j{\dot{V}}_P/\alpha$, andFor example, $K_{kj}{\partial }_j(e_i{e}_k/R)=[a_i+e_i(s-3m)]/R^2$, while differentiating the near-field angular tensor gives $[6a_i+3e_{i}s-15e_{i}m]/R^4$. Combining these with the endpoint terms in ${\partial }_{R}I$ yieldsFor the normal crack, $\mathbf{a}=\lambda \mathbf{e}+2\mu \gamma \mathbf{n}$, $m=\lambda +2\mu {\gamma }^2$, and $s=3\lambda +2\mu$. Explicitly, the three nonradiative vector coefficients in the last two lines areThus those lines are ${\mathbf{C}}_P{V}_P/({\alpha }^2{R}^2)+{\mathbf{C}}_S{V}_S/({\beta }^2{R}^2)+{\mathbf{C}}_{I}I/R^4$. The small-source approximation assumes the patch size is small relative to both $R$ and the wavelengths of interest; a finite crack is obtained by integrating the same point-patch expression over its surface with its local opening and normal.

The first two terms are the radiative P wave and S wave displacements:They decay as $R^{-1}$ and are purely transient: once the opening rate vanishes, each disappears after its respective arrival time. The P displacement is longitudinal and the S displacement transverse. An abrupt opening history makes these formulas distributions; a smooth finite-duration intrusion gives ordinary pulses.

Suppose the source starts from zero opening and ends with volume change $V_0$. At times later than the end of intrusion plus $R/\beta$, both retarded openings equal $V_0$ andThe remaining terms combine to the static displacement fieldIt decays as $R^{-2}$, so it must not be identified with the transient far-field pulses.

For a Poisson solid, $\lambda =\mu$, ${\alpha }^2=3{\beta }^2$ and $\mu =\rho {\beta }^2$. Let $\theta$ be polar angle from the crack normal, with unit vector ${\mathbf{e}}_{\theta }$ in the increasing-$\theta$ direction. Then $\gamma =\cos \theta$ and $\mathbf{n}-\gamma \mathbf{e}=-\sin \theta {\mathbf{e}}_{\theta }$. The tensile-crack radiation pattern becomesAt fixed $R$ and positive opening rate, the P pattern is outward, has no angular node, and is three times larger along the normal than in the crack plane. The signed S amplitude has four alternating lobes: it vanishes along the normal and in the crack plane, and its magnitude is maximal at $45^{\circ }$ to them. Neither field has an azimuthal SH component about the crack normal.

The final static field simplifies toIts radial part is outward near the poles, changes sign at ${\cos }^2\theta =1/8$, and is inward in the crack plane. The tangential part has the same signed angular factor as the S radiation, but is permanent and has a different distance dependence. Its radial sign changes are not nodes of the whole vector because the tangential component generally remains nonzero there. Rotating a meridional sketch around the crack normal gives the complete three-dimensional axisymmetric pattern.

The arrows show displacement at points on a circle of constant observation radius. Each panel uses its own scale; the separate formulas retain the physical P/S prefactors and arrival times. The normal is vertical and the crack lies horizontally.