Past exam of the mathematics course of the University of Cambridge / 2004 / iii / Paper 79 (Seismic Waves)

Choose depth $z$ positive downwards and harmonic time dependence $e^{-i\omega t}$. At normal incidence, P waves and the two polarizations of S waves are independent scalar problems: a P wave uses the longitudinal modulus $\rho {\alpha }^2$, and an S wave uses the shear modulus $\rho {\beta }^2$. In a uniform layer with wave speed $c$ and seismic impedance $q=\rho c$, write velocity and compressive traction asAt a perfectly bonded interface the velocity and traction are continuous. For incident velocity wave amplitude one from medium 1 into medium 2, these two conditions are $1+r_v=t_v$ and $q_1(1-r_v)=q_2{t}_v$. ThusReflection is determined by seismic impedance, not by speed alone. A step to larger seismic impedance reverses the reflected velocity/displacement field waveform; the pressure reflection coefficient has the opposite sign. Its transmitted pressure coefficient is $t_p=2q_2/(q_1+q_2)$. Equal seismic impedance eliminates the reflection even when speed and wavelength change. At a free surface traction vanishes and velocity reflects with coefficient $+1$; at a rigid surface velocity reflects with coefficient $-1$.

The mean elastic-wave energy flux is $\frac{1}{2}\mathrm{Re}(P{V}^{\ast })=q(|A{|}^2-|B{|}^2)/2$ for positive real $q$. Therefore the reflected and transmitted energy fractions areA displacement field transmission coefficient greater than one is compatible with conservation of energy; its flux must be weighted by the receiving seismic impedance.

For many layers, use the seismic layer transfer matrix. A layer of thickness $h_j$ has wave phase thickness ${\delta }_j=\omega h_j/c_j$ andThis follows either by eliminating $A,B$ between the two faces or by integrating $V^{\prime }=i\omega P/(\rho c^2)$ and $P^{\prime }=i\omega \rho V$. For the complete stack put $M=M_N\cdots M_1=\left(\begin{array}{cc}A& B\\ C& D\end{array}\right)$. With upper and lower half-space impedances $q_a,q_b$, an outgoing lower wave satisfies $P_b=q_b{V}_b$. Solving this condition gives the input seismic impedanceHere $T$ is lower-face velocity divided by upper-face incident velocity, including the propagation wave phase. For lossless layers and real frequency, $|R{|}^2+(q_b/q_a)|T{|}^2=1$. This supplies a sensitive check on any spreadsheet or numerical implementation; it is not correct simply to multiply the individual interface transmissions and discard the elastic waves.

The physical reason is repeated wave interference. In one layer between media 0 and 2, let $r_{ij},t_{ij}$ denote the velocity coefficients above and $\delta =\omega h/c_1$. The first transmitted wave has wave amplitude $t_{01}{t}_{12}{e}^{i\delta }$. Every extra round trip supplies $r_{10}{r}_{12}{e}^{2i\delta }$. Summing the geometric series givesIn the time domain the successive arrivals are delayed by $2h/c_1$, producing reverberation. In the frequency domain the same delayed terms produce resonant maxima and minima. For a layer between identical half-spaces of impedance $q_0$,At $\delta =n\pi$ the layer transmits all the incident energy: the internally elastic waves cancel the direct reflection. Halfway between these resonances, a large impedance contrast gives strong reflection. For unequal exterior impedances, a quarter-wave layer with $q_1^2=q_a{q}_b$ transforms its lower load to $Z_{\mathrm{in}}=q_1^2/q_b=q_a$, again giving zero reflection at its design frequency. This matching is narrow-band.

Periodic alternation produces a periodic seismic multilayer stop band. For a cell of two layers with total thickness $d$, the cell matrix has reciprocal eigenvalues $e^{\pm iKd}$, soWhen the absolute value of this expression exceeds one, the infinite periodic medium has no propagating cell mode; a long finite stack has exponentially small transmission. The energy is reflected, not absorbed. Pass bands instead contain oscillatory transmission peaks, becoming more numerous as the stack length increases. Disorder removes the exact cell periodicity but not multiple scattering or frequency-selective interference.

For a concrete reproducible investigation, take $q_H/q_L=4$, equal one-way travel times $\tau$ in both layer types, exterior impedance $q_L$, and frequency coordinate $\nu =2\omega \tau /\pi$. At $\nu =1$ each layer is a quarter wavelength; one cell has diagonal matrix entries $-4,-1/4$. For $N$ cells,It is about $0.22145$ for one cell, $6.10\times 10^{-5}$ for four cells and $9.31\times 10^{-10}$ for eight cells. The infinite-cell stop band is $2\arcsin (0.8)/\pi <\nu <2-2\arcsin (0.8)/\pi$. The accompanying calculation also shows the first transmitted arrivals for a single high-impedance slab: their wave amplitudes are $0.64(0.36{)}^n$ at times $(2n+1)\tau$. Their sum is one, consistent with unit zero-frequency transmission. These calculations illustrate both the spectrum and the pulse train without assuming access to the course's spreadsheet.

Finally, as $\omega \to 0$ for a fixed finite stack, every $M_j\to I$, so the response approaches that of a direct interface between the exterior half-spaces. Equal exterior impedances give $R\to 0,T\to 1$. For wavelengths much longer than a repeating cell, expansion of $M_j$ to first order instead gives an effective medium with averaged mass density ${\rho }_{\mathrm{eff}}=\sum h_j{\rho }_j/d$ and harmonic-mean modulus $M_{\mathrm{eff}}^{-1}=\sum h_j/M_j^{\mathrm{mod}}/d$, where $M_j^{\mathrm{mod}}={\rho }_j{c}_j^2$. Its speed is $\sqrt{M_{\mathrm{eff}}/{\rho }_{\mathrm{eff}}}$, not generally the arithmetic mean of the layer speeds. Attenuation makes the energy balance an inequality, while for a fixed oblique horizontal slowness the scalar SH calculation uses $q_j={\mu }_j{p}_{\beta ,j}$ and ${\delta }_j=\omega p_{\beta ,j}{h}_j$. Oblique P-SV waves require coupled displacement field and traction matrices because of mode conversion at a planar elastic interface. Thick evanescent layers require scattering or impedance recursion rather than exponentially ill-conditioned transfer-matrix products.

Work in a meridional $R,z$ plane and denote the interface height by $Z(R)$, to distinguish it from an observation coordinate. PutThe interface is $Z=(\alpha s+{\alpha }_{f}D)/d^2$, so $Z^{\prime }={\alpha }_{f}R/D$. Its unit normal into the fluid is proportional to $(-Z^{\prime },1)$. Consequently the angle between this normal and the vertical satisfiesThe outgoing source P wave has ray direction proportional to $(R,Z)$, tilting outwards by $\arctan (R/Z)$ from the vertical, while the normal tilts inwards. Thus ${\theta }_P=\arctan (R/Z)+{\theta }_z$. The tangent addition formula givesDirect substitution simplifies its denominator to $s(\alpha D+{\alpha }_{f}s)/(d^{2}D)$ and its numerator to $\alpha R(\alpha D+{\alpha }_{f}s)/(d^{2}D)$. HenceThe denominator is positive, so these are the physical incidence angles in $[0,\pi /2)$; the transmitted P wave reaches the fluid rather than propagating back into the solid.

The angles from part (i) have sinesBy Snell law for elastic and acoustic waves, $\sin {\theta }_f/{\alpha }_f=\sin {\theta }_P/\alpha$. Both angles lie in $[0,\pi /2)$, so ${\theta }_f={\theta }_z$. The outgoing fluid ray tilts away from the inward-tilted normal by exactly its tilt from vertical; it therefore travels vertically upwards. This is hyperbolic interface collimation of P-waves.

The source-to-interface distance satisfiesOne can verify this by squaring, using $D^2-s^2=d^2{R}^2$. ConsequentlyAdding the vertical fluid travel time from the interface to height $z$ givesThus surfaces of equal arrival time in the fluid are horizontal planes. Their leading wave amplitude may vary with $R$: a planar wavefront does not imply a spatially uniform plane-wave solution at every order.

Use local coordinates with the interface $z=0$, solid below and fluid above, and choose the P-SV displacement potentials convention $u_x={\varphi }_x-{\psi }_z$, $u_z={\varphi }_z+{\psi }_x$. Let the reflected P and SV potential wave amplitudes be $A,B$, while the transmitted P wave amplitude is $T_P$. All four waves have the same tangential slownessThe incident and transmitted retarded times have negative $q_{P}z,q_{f}z$, while elastic waves have positive $q_{P}z,q_{S}z$. At the interface all waveform derivatives therefore have the same argument. displacement field continuity givesWith Lamé parameters $\mu =\rho {\beta }^2$ and $\lambda =\rho ({\alpha }^2-2{\beta }^2)$, the shear and normal stresses follow directly from differentiation:The fluid has no shear traction, so the coefficient of $f^{\prime \prime }$ in ${\sigma }_{xz}$ must vanish:The coefficient of $f^{\prime \prime }$ in the solid stress is $\rho [\cos 2{\theta }_S(1+A)-\sin 2{\theta }_{S}B]$; in the fluid it is ${\rho }_f{T}_P$. ThusThese are the three required matching conditions; tangential displacement field need not be continuous at a solid–inviscid-fluid boundary.

Put $C=\cos 2{\theta }_S$ and $X=1-A$. From (2), $B=({\beta }^2/{\alpha }^2)(\sin 2{\theta }_P/C)X$. In (1), the coefficient of $X$ simplifies because $C+2({\beta }^2/{\alpha }^2){\sin }^2{\theta }_P=1$. Hence $X=(\alpha C\cos {\theta }_f)/({\alpha }_f\cos {\theta }_P)T_P$. Substitute this into (3) to obtainTherefore the solid-fluid P-wave transmission coefficient isDivision by $C$ was only an elimination device; the final expression extends continuously to $C=0$. The normal-incidence value is $2\rho {\alpha }_f/(\rho \alpha +{\rho }_f{\alpha }_f)$. This is a potential coefficient; longitudinal displacement field transmission additionally multiplies it by $\alpha /{\alpha }_f$. The original PDF's bracket lies in the denominator and contains ${\beta }^2/{\alpha }^2$ and $\cos {\theta }_f$. The converted TeX misplaces the bracket, substitutes ${\theta }^2$ for ${\beta }^2$, and loses the final cosine.

The leading source displacement field isThe second term is one derivative smoother at the leading arrival. At each surface point, use its tangent plane, the angles from (i)–(ii), and the local solid-fluid P-wave transmission coefficient. The incident potential coefficient is $1/r(R)$, so the transmitted leading potential coefficient isTransport this coefficient along the vertical fluid rays. Their transverse area does not change with height, so ray-tube conservation for a P-wave jump gives no further geometrical spreading. The leading displacement field coefficient is consequentlyThis is the transmitted jump at a collimating solid-fluid interface. Radial derivatives of $A_{\varphi }$ multiply $f$, and therefore do not contribute to the first derivative jump carried by $f^{\prime }$. A varying transverse wave amplitude does generate smoother later terms, so this construction is a leading-front result, not a complete global plane-wave field.

Since the source is zero before $t=0$ and is $C^2$ there, $f(0)=f^{\prime }(0)=f^{\prime \prime }(0)=0$. Thus displacement field itself does not jump. If $J=[f^{(3)}(0)]$ exists and is the first nonzero jump, the leading discontinuity is the acceleration jump $[{\partial }_t^2\mathbf{u}]=-A_u(R)J{\mathbf{e}}_z$. More generally, a first jump in $f^{(m)}$ gives $[{\partial }_t^{m-1}\mathbf{u}]=-A_u(R)[f^{(m)}]{\mathbf{e}}_z$. The supplied regularity alone does not guarantee a nonzero discontinuity: a source smooth to all orders may have none. The geometrical transfer factor $A_u$ is nevertheless determined.

At $R=0$, $r=z_0$ and all three angles vanish. HenceFor large $R$, put $C_{\infty }=1-2{\beta }^2/{\alpha }^2$ and retain $s=z_0(\alpha -{\alpha }_f)$, $d=\sqrt{{\alpha }^2-{\alpha }_f^2}$. The relevant limits areWhen $C_{\infty }\ne 0$, the square-bracket factor in the denominator of $T_P$ is $C_{\infty }^2+O(R^{-1})$, soThese are signed waveform coefficients; their absolute values give wave amplitude magnitudes. The leading displacement field decays as $R^{-2}$ rather than $R^{-1}$ because grazing incidence adds a small transmission factor to spherical spreading. The leading coefficient is independent of fluid and solid mass densities, although the exact wave amplitude and its higher-order terms depend on them.

The assumptions do not exclude ${\beta }^2={\alpha }^2/2$, where $C_{\infty }=0$ and the preceding division is invalid. In this special case $\cos 2{\theta }_S={\cos }^2{\theta }_P$ exactly, and the square-bracket factor is $\cos {\theta }_P+O({\cos }^3{\theta }_P)$. Keeping both denominator terms givesIt has no $R^{-2}$ term. This exceptional case matters when requesting an asymptotic formula without an additional restriction on the elastic moduli.

Write $\omega =u+iv$, with $v>0$, and a chosen square root $k_{\beta }=a+ib$. Squaring givesOne cannot have $b=0$: the second relation would force $u=0$, and then the first would require the square of a real number to be negative. ConsequentlyThe bracket has the sign of $b$, proving the causal vertical wavenumber branch equivalenceThis selects the square root analytically in the upper frequency half-plane. At $\omega =iv$ it is $k_{\beta }=i\sqrt{v^2/{\beta }^2+k^2}$; continuation supplies the outgoing branch for other frequencies. The expression $\sqrt{{\beta }^{-2}-p^2}$ must use this inherited branch rather than an independently chosen square root. For positive real frequency its limit is positive real for a propagating wave and positive imaginary for an evanescent wave.

Let the shear-horizontal wave displacement field be $U(z)e^{ikx-i\omega t}$. Its velocity is $V=-i\omega U$, and its traction is $T=\mu U^{\prime }$. The scalar shear-wave equation givesWith $p=k/\omega$ and $p_{\beta }^2={\beta }^{-2}-p^2$, the first-order velocity–traction system isThese equations remain valid for depth-dependent $\mu ,\rho$; differentiating $\mu U^{\prime }$ is essential and one must not incorrectly use $\mu U^{\prime \prime }$ when $\mu$ varies.

Where $V\ne 0$, define the SH input impedance $Z=-T/V$. Differentiating $T=-ZV$ and substituting the two equations gives the SH impedance Riccati equation:The second relation is a one-way scalar equation once the complete SH input impedance is known. It does not discard reflections: $Z$ describes the actual superposition of both wave directions. At a zero of $V$, $Z$ may have a pole while the velocity–traction system stays regular.

Use the quadratic spectral flux $S_z=-\frac{1}{2}\mathrm{Re}(V^{\ast }T)$. For positive real material mass density and shear modulus, differentiation of the velocity–traction system givesBecause $k$ is real, $\mu p_{\beta }^2=\rho -\mu k^2/{\omega }^2$. If $v=\mathrm{Im}\omega >0$, then $\mathrm{Re}(i{\omega }^{\ast })=v$ and $\mathrm{Re}(-i\omega \mu p_{\beta }^2)=v(\rho +\mu k^2/|\omega {|}^2)$. This proves complex-frequency SH flux monotonicity:Strictness holds for a nontrivial wavefield: if $V,T$ vanished together at one point, uniqueness of the first-order system would make the entire solution zero. The zero solution is the necessary exception to the literal strict inequality. Across a bonded material boundary $V,T$, and hence $S_z$, are continuous; integrating the negative derivative gives strict decrease across any positive-depth interval as well.

Where the SH input impedance is finite, $S_z=\frac{1}{2}\mathrm{Re}Z|V{|}^2$. If $\mathrm{Re}Z(z_1)>0$, then $S_z(z_1)>0$ and all shallower fluxes are larger and positive. They cannot have $V=0$, since that would make $S_z=0$. Thus $\mathrm{Re}Z(z)>0$ for every $z<z_1$. Likewise, negative real SH input impedance at $z_0$ gives negative flux throughout $z>z_0$, excludes velocity zeros there and implies $\mathrm{Re}Z(z)<0$ for every $z>z_0$.

In a uniform region let $q=\mu p_{\beta }$ on the causal vertical wavenumber branch. The constant Riccati solutions are $Z=\pm q$, and integrating the one-way equation givesSince $\mathrm{Re}q>0$, the plus solution carries positive downward flux and the minus solution negative upward flux. Writing $k_{\beta }=a+ib$ with $b>0$, their moduli are proportional to $e^{-bz}$ and $e^{bz}$ respectively: each decays in its direction of energy flow.

For complex frequency the displayed quadratic form is the spectral flux diagnostic specified in the question; a genuinely periodic time average is recovered on the real-frequency axis. The decrease here comes from the upper-half-plane convergence parameter, not from material dissipation. At real frequency in a lossless layer the corresponding flux is constant.

Let $z_r$ be either endpoint at which $Z_r$ is prescribed, and put $d_z=z-z_r$, $q=\mu p_{\beta }$ and $\kappa =\omega p_{\beta }$. Integrating the constant-coefficient velocity–traction system givesTaking the ratio $-T/V$ yields the uniform-layer SH impedance update and its accompanying velocity:The initial value $V_r$ sets an arbitrary normalization; specifying impedance alone cannot specify the overall wave wave amplitude. Positive $d_z$ propagates data from $z_0$, while negative $d_z$ propagates a specified load from $z_1$, so the same expressions solve both requested initial-value problems.

Expanding sine and cosine givesAbsorbing the reference-point phases into the two constants gives the requested forward/backward exponential superposition. Equivalently, wherever the ratios are defined,The exceptional values $Z_r=q,-q$ give the pure directional solutions. If a denominator in the impedance formula vanishes, propagate the finite pair $(V,T)$ through that point; this is an impedance-chart pole, not a failure of the physical solution. The upper-half-plane assumptions ensure $q\ne 0$.

In the lower uniform half-space the causal vertical wavenumber branch has $\mathrm{Im}{k}_{\beta ,b}>0$. The downward factor $e^{i{k}_{\beta ,b}(z-z_b)}$ decays as $z\to +\infty$, whereas the upward factor grows there and represents an incoming wave from the lower infinity. The outgoing radiation condition therefore retains only the downward factor and givesFor each uniform layer $j$ of thickness $h_j$, work upwards from its known bottom load $Z_{j,b}$. With $q_j={\mu }_j{p}_{\beta ,j}$ and ${\delta }_j=\omega p_{\beta ,j}{h}_j$, part (iv) givesvelocity and traction are continuous across a bonded interface, so this top impedance becomes the bottom load of the next shallower layer. Repeat up to $z_a$ to obtain $Z_{\mathrm{in}}=Z(z_a)$. Poles can be handled by carrying $(V,T)$ or the reciprocal ratio. Since the lower load has positive real part, part (iii) actually ensures that the outgoing whole-stack impedance stays finite with positive real part at all shallower depths in the upper frequency half-plane.

The upper half-space has directional impedance $q_a={\mu }_a{p}_{\beta ,a}$. At $z_a$, an incident velocity $V_0$ and reflected velocity $R{V}_0$ giveTaking their ratio and solving establishes reflection from a layered SH impedance:Thus no individual deeper reflection coefficients are needed after $Z_{\mathrm{in}}$ is known: it contains all underlying interfaces and their coherent multiples. A matched load gives $R=0$; the real-frequency limits of zero traction and zero velocity give $R=+1$ and $R=-1$. For real frequency and propagating upper and lower half-space waves, the usual energy-flux balance provides a further check on the computed reflection and transmission. A stable implementation for thick layers uses the ratio update or its decaying-exponential version instead of multiplying growing transfer matrices.