Wave Propagation

This section covers the mathematical treatment of seismic wave propagation through layered media, which forms the computational foundation of site response analysis.

Governing Equations

One-Dimensional Wave Equation

For vertically propagating shear waves in a layered medium, the equation of motion is:

\[\rho \frac{\partial^2 u}{\partial t^2} = G \frac{\partial^2 u}{\partial z^2}\]

where: - \(u(z,t)\) is the horizontal displacement - \(\rho\) is the mass density - \(G\) is the shear modulus - \(z\) is the vertical coordinate (positive upward)

Complex Shear Modulus

To incorporate material damping, the shear modulus is treated as complex:

\[G^* = G(1 + 2i\xi)\]

where \(\xi\) is the damping ratio.

Solution Methods

Frequency Domain Approach

Taking the Fourier transform of the wave equation:

\[-\omega^2 \rho \hat{u} = G^* \frac{\partial^2 \hat{u}}{\partial z^2}\]

This leads to the general solution:

\[\hat{u}(z,\omega) = A e^{ik^*z} + B e^{-ik^*z}\]

where \(k^* = \omega/V_s^*\) is the complex wavenumber and \(V_s^* = \sqrt{G^*/\rho}\).

Transfer Matrix Method

For layered systems, the Thomson-Haskell propagator matrix method relates the wave field at different depths. For a single layer of thickness \(h\), the transfer matrix is:

\[\begin{split}\mathbf{T} = \begin{pmatrix} \cos(k^* h) & \frac{i \sin(k^* h)}{Z^*} \\ i Z^* \sin(k^* h) & \cos(k^* h) \end{pmatrix}\end{split}\]

where \(Z^* = \rho V_s^*\) is the complex impedance.

Boundary Conditions

Free Surface

At the ground surface (\(z = 0\)), the shear stress must vanish:

\[\tau = G^* \frac{\partial u}{\partial z} = 0\]

Input Motion Specification

Three common approaches for specifying input motion:

  1. Outcrop Motion: Motion that would occur at a rock outcrop

  2. Within Motion: Motion recorded within the rock formation

  3. Incident Motion: Upward-traveling wave component only

Radiation Damping

For semi-infinite elastic bedrock, the radiation condition requires:

\[\tau = \rho V_s u̇\]

This represents energy radiation away from the site.

Numerical Implementation

Frequency Sampling

Adequate frequency resolution is critical for accurate results:

\[\Delta f \leq \frac{V_{s,min}}{8 H_{max}}\]

where \(V_{s,min}\) is the minimum shear wave velocity and \(H_{max}\) is the maximum depth of interest.

Stability Considerations

The method is stable for all frequencies, but numerical precision may be lost for: - Very thick, soft layers - High frequencies - Strong impedance contrasts

Efficiency Optimizations

  • Pre-compute layer matrices for repeated analyses

  • Use symmetry to reduce computation for real-valued inputs

  • Implement fast convolution for time domain conversion

Validation Examples

Analytical Solutions

  1. Uniform Half-Space

    For a uniform elastic half-space, the exact amplification function is:

    \[H(\omega) = \frac{2Z_0}{Z_0 + Z_r}\]

    where \(Z_0\) is the surface impedance and \(Z_r\) is the reference impedance.

  2. Single Layer over Rigid Base

    The transfer function has resonant peaks at:

    \[f_n = \frac{(2n-1)V_s}{4H}\]

    for \(n = 1, 2, 3, ...\)

Benchmark Comparisons

PyStrata wave propagation algorithms have been validated against: - Analytical solutions for simple geometries - SHAKE91 for equivalent linear analysis - DEEPSOIL for advanced nonlinear methods - Recorded earthquake data from instrumented sites

Implementation Details

The wave propagation calculations in PyStrata are implemented in the propagation module with careful attention to:

  • Numerical stability through appropriate algorithmic choices

  • Computational efficiency via optimized matrix operations

  • Physical accuracy through proper boundary condition treatment

  • Flexibility to accommodate various input motion types and analysis methods