Part II: Geophysics | Chapter 2

Seismology

Seismic waves, travel-time curves, and imaging Earth's interior

2.1 The Seismic Wave Equation

Seismic waves propagate through Earth as elastic disturbances governed by Newton's second law applied to a continuous medium. The equation of motion for an elastic solid relates acceleration to the divergence of the stress tensor:

$$\rho \frac{\partial^2 u_i}{\partial t^2} = \frac{\partial \sigma_{ij}}{\partial x_j} + f_i$$

where $\rho$ is density, $u_i$ is displacement, and $f_i$ is body force. For an isotropic, linearly elastic solid, Hooke's law gives:

$$\sigma_{ij} = \lambda \delta_{ij} \varepsilon_{kk} + 2\mu \varepsilon_{ij}$$

where $\lambda$ and $\mu$ are the Lame parameters and $\varepsilon_{ij} = \frac{1}{2}(\partial_i u_j + \partial_j u_i)$.

Derivation 1: P-wave and S-wave Velocities

Substituting Hooke's law into the equation of motion (neglecting body forces):

$$\rho \frac{\partial^2 \mathbf{u}}{\partial t^2} = (\lambda + \mu)\nabla(\nabla \cdot \mathbf{u}) + \mu \nabla^2 \mathbf{u}$$

Using the Helmholtz decomposition $\mathbf{u} = \nabla\phi + \nabla \times \boldsymbol{\psi}$, the displacement field separates into irrotational and solenoidal parts. Taking the divergence of the equation of motion:

$$\frac{\partial^2 \theta}{\partial t^2} = \alpha^2 \nabla^2 \theta, \quad \alpha = \sqrt{\frac{\lambda + 2\mu}{\rho}} = V_P$$

where $\theta = \nabla \cdot \mathbf{u}$ is the dilatation. This is the P-wave(primary, compressional) equation. Taking the curl:

$$\frac{\partial^2 \boldsymbol{\Omega}}{\partial t^2} = \beta^2 \nabla^2 \boldsymbol{\Omega}, \quad \beta = \sqrt{\frac{\mu}{\rho}} = V_S$$

where $\boldsymbol{\Omega} = \nabla \times \mathbf{u}$ is the rotation. This is the S-wave(secondary, shear) equation. Since $\lambda + 2\mu > \mu$ always, P-waves travel faster than S-waves. The typical ratio is:

$$\frac{V_P}{V_S} = \sqrt{\frac{\lambda + 2\mu}{\mu}} = \sqrt{\frac{2(1-\nu)}{1-2\nu}}$$

where $\nu$ is Poisson's ratio. For a Poisson solid ($\lambda = \mu$, $\nu = 0.25$),$V_P/V_S = \sqrt{3} \approx 1.73$.

Historical Context

The mathematical theory of elastic waves was developed by Poisson (1831), Stokes (1849), and Lord Rayleigh (1885). The first seismograph was built by John Milne in 1880. Richard Dixon Oldham identified P, S, and surface waves in earthquake records in 1906, and Beno Gutenberg located the core-mantle boundary at 2900 km depth in 1914 using the shadow zone of P-waves.

2.2 Surface Waves

Derivation 2: Rayleigh Wave Dispersion

Lord Rayleigh (1885) showed that elastic waves can propagate along a free surface. The Rayleigh wave is a combination of P and SV motions that decays exponentially with depth. For a half-space, the Rayleigh wave velocity $V_R$ satisfies:

$$\left(\frac{V_R}{V_S}\right)^6 - 8\left(\frac{V_R}{V_S}\right)^4 + \left(24 - 16\frac{V_S^2}{V_P^2}\right)\left(\frac{V_R}{V_S}\right)^2 - 16\left(1 - \frac{V_S^2}{V_P^2}\right) = 0$$

For a Poisson solid ($V_P/V_S = \sqrt{3}$), numerical solution gives $V_R \approx 0.9194 V_S$. In general, $V_R$ is slightly less than $V_S$, ranging from 0.87 to $0.96 V_S$depending on Poisson's ratio.

The displacement amplitudes decay exponentially with depth $z$:

$$u_x(z) \propto e^{-q_\alpha k z} - C e^{-q_\beta k z}, \quad u_z(z) \propto -q_\alpha e^{-q_\alpha k z} + C' e^{-q_\beta k z}$$

where $q_\alpha = \sqrt{1 - V_R^2/V_P^2}$, $q_\beta = \sqrt{1 - V_R^2/V_S^2}$, and $k$is the wavenumber. The particle motion is retrograde elliptical at the surface and changes to prograde at depth ~0.2 wavelengths.

Love Waves

A.E.H. Love (1911) showed that horizontally polarized shear waves can be trapped in a low-velocity surface layer. Love waves exist only when the surface layer has lower shear velocity than the substrate. The dispersion relation for a layer of thickness $H$:

$$\tan\left(\omega H \sqrt{\frac{1}{V_{S1}^2} - \frac{1}{c^2}}\right) = \frac{\mu_2 \sqrt{\frac{1}{c^2} - \frac{1}{V_{S2}^2}}}{\mu_1 \sqrt{\frac{1}{V_{S1}^2} - \frac{1}{c^2}}}$$

where $c$ is the phase velocity, constrained between $V_{S1} < c < V_{S2}$. Love waves are dispersive: longer periods sample deeper structure and travel faster. This dispersion is exploited in surface wave tomography to image crustal and upper mantle structure.

2.3 Seismic Ray Theory and Travel Times

Derivation 3: Snell's Law and the Ray Parameter

In the high-frequency limit, seismic energy propagates along rays governed by Fermat's principle (least travel time). At an interface between layers with velocities $v_1$ and $v_2$, Snell's law holds:

$$p = \frac{\sin i_1}{v_1} = \frac{\sin i_2}{v_2} = \text{const}$$

where $p$ is the ray parameter (horizontal slowness), conserved along the entire ray path. For a spherically symmetric Earth with radius $r$ and velocity $v(r)$:

$$p = \frac{r \sin i}{v(r)} = \frac{r}{v(r)} \bigg|_{r = r_{\text{tp}}} = \frac{r_{\text{tp}}}{v(r_{\text{tp}})}$$

At the turning point $r_{\text{tp}}$, the ray is horizontal ($i = 90°$). The travel time for a ray from source to receiver at epicentral distance $\Delta$ is:

$$T(\Delta) = 2 \int_{r_{\text{tp}}}^{R} \frac{dr}{v(r)\sqrt{1 - p^2 v^2(r)/r^2}}$$

The slope of the travel-time curve gives the ray parameter: $p = dT/d\Delta$. This is the basis for velocity inversion: measuring travel times at many distances constrains the velocity profile $v(r)$.

Figure 1. Seismic wave ray paths through Earth's interior. P-waves (solid lines) refract through all layers, while S-waves (dashed) cannot propagate through the liquid outer core, creating the shadow zone at 103deg-143deg epicentral distance.

Shadow Zones and Core Phases

The liquid outer core creates a P-wave shadow zone between approximately 104° and 140° epicentral distance because S-waves cannot propagate through liquid, and P-waves are refracted. Key core phases include PKP (P through core), PKiKP (reflected off inner core), and PKIKP (through inner core). Inge Lehmann (1936) discovered the solid inner core by identifying arrivals within the shadow zone that could only be explained by reflection off an inner boundary.

2.4 Seismic Tomography

Derivation 4: The Tomographic Inverse Problem

Seismic tomography images 3D velocity structure from travel-time residuals. A travel-time anomaly $\delta T$ arises from velocity perturbations $\delta v$ along the ray path:

$$\delta T = T_{\text{obs}} - T_{\text{pred}} = -\int_{\text{ray}} \frac{\delta v(\mathbf{r})}{v_0^2(\mathbf{r})} ds$$

Parameterizing the velocity perturbation as a sum over basis functions (e.g., blocks or spherical harmonics): $\delta v(\mathbf{r}) = \sum_j m_j \phi_j(\mathbf{r})$, we obtain a linear system:

$$\delta T_i = \sum_j G_{ij} m_j, \quad G_{ij} = -\int_{\text{ray}_i} \frac{\phi_j(\mathbf{r})}{v_0^2(\mathbf{r})} ds$$

In matrix form: $\mathbf{d} = \mathbf{G}\mathbf{m}$. The system is typically underdetermined and ill-conditioned, requiring regularization. The damped least-squares solution is:

$$\mathbf{m} = (\mathbf{G}^T \mathbf{G} + \epsilon^2 \mathbf{I})^{-1} \mathbf{G}^T \mathbf{d}$$

This technique has revealed large-scale structures including subducted slabs (fast anomalies) and mantle plumes (slow anomalies). Global P-wave tomography models (e.g., Bijwaard et al., 1998; Li et al., 2008) resolve structures down to ~200 km.

Derivation 5: Earthquake Location from Arrival Times

An earthquake is located by minimizing travel-time residuals across a network of stations. The S-P time difference at station $i$ gives an estimate of distance:

$$d_i = \frac{t_{S,i} - t_{P,i}}{1/V_S - 1/V_P}$$

For precise location, we linearize the problem around an initial guess $(\mathbf{r}_0, t_0)$and iterate. The predicted arrival time at station $i$ is:

$$t_i^{\text{pred}} = t_0 + T(\mathbf{r}_0, \mathbf{r}_i) + \frac{\partial T}{\partial \mathbf{r}} \cdot \delta\mathbf{r} + \delta t_0$$

The Geiger method (1910) iteratively solves this least-squares problem, adjusting the hypocenter and origin time until residuals are minimized. Modern methods use 3D velocity models and waveform cross-correlation for sub-kilometer precision.

2.5 Earthquake Magnitude and Moment

Charles Richter (1935) defined the local magnitude scale:

$$M_L = \log_{10}(A) - \log_{10}(A_0(\delta))$$

where $A$ is maximum amplitude and $A_0$ is a reference amplitude at distance $\delta$. The moment magnitude, preferred for large earthquakes, is based on seismic moment:

$$M_0 = \mu A D$$

where $\mu$ is rigidity, $A$ is fault area, and $D$ is average slip. The moment magnitude is:

$$M_w = \frac{2}{3}(\log_{10} M_0 - 9.1)$$

Kanamori (1977) introduced this scale, which does not saturate for great earthquakes. The 2011 Tohoku earthquake ($M_w$ 9.1) released a seismic moment of $5.3 \times 10^{22}$ N$\cdot$m, equivalent to rupturing ~500 km of fault with ~30 m slip.

Energy-Magnitude Relationship

The radiated seismic energy is related to moment magnitude by:

$$\log_{10} E_s = 1.5 M_w + 4.8$$

Each unit increase in magnitude corresponds to a factor of $10^{1.5} \approx 32$increase in radiated energy. An M8 earthquake releases about 1000 times more energy than an M6. The total energy budget includes both radiated seismic energy (typically only 1-10% of the total) and fracture energy consumed in creating new fault surface and overcoming friction.

Earthquake Statistics and Hazard

Globally, there are approximately 15 earthquakes per year with $M_w \geq 7$, one with $M_w \geq 8$, and a $M_w \geq 9$ every 50-100 years. The largest recorded earthquake is the 1960 Chilean earthquake ($M_w$ 9.5), which ruptured ~1000 km of the Nazca-South American plate boundary. The physical upper limit on earthquake size is set by the maximum fault area and stress drop:$M_0 = \mu A D \leq \mu L W \bar{D}$, where $L$ is fault length, $W$ is width, and $\bar{D}$ is average slip.

2.6 Focal Mechanisms and Moment Tensors

The pattern of P-wave first motions (compressions and dilatations) at different stations constrains the orientation of the fault plane and slip direction. A double-couple source produces four quadrants of alternating polarity separated by two nodal planes. The seismic moment tensor $M_{ij}$ provides a complete description of the source:

$$M_{ij} = M_0(\hat{s}_i \hat{n}_j + \hat{s}_j \hat{n}_i)$$

where $\hat{n}$ is the fault normal and $\hat{s}$ is the slip vector. The "beachball" diagram represents the focal mechanism projected onto the lower focal hemisphere. Normal faults show an extensional T-axis and white quadrants on the sides; reverse faults show compressional P-axis and dark quadrants on the sides; strike-slip faults show a pattern of four equal quadrants.

The Global CMT (Centroid Moment Tensor) Project, founded by Adam Dziewonski in 1976, has catalogued over 50,000 earthquake mechanisms. The catalog reveals the global pattern of tectonic stress: normal faulting at mid-ocean ridges, reverse faulting at subduction zones, and strike-slip faulting along transform boundaries. Non-double-couple components in the moment tensor indicate complex rupture, volcanic sources, or landslides.

Earthquake Source Spectrum

The far-field displacement spectrum of an earthquake has a flat level at low frequencies proportional to $M_0$ and falls off as $\omega^{-2}$ above the corner frequency $f_c$:

$$\hat{u}(\omega) = \frac{M_0}{1 + (\omega/\omega_c)^2}$$

The corner frequency is related to source dimension: $f_c \approx 0.37 V_S / r_0$where $r_0$ is the source radius (Brune, 1970). This relationship is used to estimate stress drop $\Delta\sigma = 7M_0/(16r_0^3)$, which is typically 1-10 MPa and remarkably constant across many orders of magnitude in earthquake size.

Computational Lab: Seismology

Seismic Wave Propagation, Travel Times, and Earthquake Location

Python

script.py316 lines

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

# ============================================================
# SEISMOLOGY SIMULATION
# ============================================================
print("SEISMOLOGY SIMULATION")
print("=" * 70)

# ============================================================
# 1. P and S Wave Velocities for Common Rock Types
# ============================================================
print()
print("1. Seismic Velocities in Common Rocks")
print("-" * 65)
print(f"  {'Rock Type':20} {'Vp (km/s)':12} {'Vs (km/s)':12} {'Vp/Vs':10} {'Poisson'}")
print("  " + "-" * 60)

rocks = [
    ("Granite", 6.0, 3.5),
    ("Basalt", 6.4, 3.2),
    ("Gabbro", 7.0, 3.8),
    ("Peridotite", 8.1, 4.5),
    ("Sandstone", 3.5, 2.0),
    ("Limestone", 5.5, 3.0),
    ("Shale", 3.0, 1.5),
    ("Water", 1.5, 0.0),
    ("Upper Mantle", 8.1, 4.5),
    ("Lower Mantle", 13.6, 7.3),
    ("Outer Core", 10.0, 0.0),
    ("Inner Core", 11.2, 3.6),
]

for name, vp, vs in rocks:
    ratio = vp/vs if vs > 0 else float('inf')
    if vs > 0:
        nu = (ratio**2 - 2) / (2*(ratio**2 - 1))
        print(f"  {name:18} {vp:10.1f} {vs:10.1f} {ratio:8.2f} {nu:8.3f}")
    else:
        print(f"  {name:18} {vp:10.1f} {'---':10} {'---':8} {'0.500':8}")

# ============================================================
# 2. Travel-Time Curves for Layered Earth
# ============================================================
print()
print()
print("2. Travel-Time Curves: Direct, Reflected, and Refracted Waves")
print("-" * 65)

# Simple two-layer model
v1 = 5.5   # km/s (crust)
v2 = 8.0   # km/s (mantle)
h = 35.0   # km (crustal thickness)

# Critical distance
ic = np.degrees(np.arcsin(v1/v2))
x_crit = 2 * h * np.tan(np.radians(ic))
# Crossover distance
t_head_func = lambda x: 2*h*np.sqrt(1/v1**2 - 1/v2**2) + x/v2
t_direct_func = lambda x: x/v1
# Solve for crossover
x_cross = 2*h*np.sqrt((v2+v1)/(v2-v1))

print(f"  Layer 1: V1 = {v1} km/s, thickness = {h} km")
print(f"  Layer 2: V2 = {v2} km/s (half-space)")
print(f"  Critical angle = {ic:.1f} deg")
print(f"  Critical distance = {x_crit:.1f} km")
print(f"  Crossover distance = {x_cross:.1f} km")
print()
print(f"  {'Distance (km)':16} {'T_direct (s)':14} {'T_reflected (s)':16} {'T_refracted (s)':16} {'First arrival'}")
print("  " + "-" * 75)

for x in [10, 20, 40, 60, 80, 100, 120, 150, 200, 250, 300]:
    t_dir = x / v1
    t_ref = np.sqrt(x**2 + (2*h)**2) / v1
    t_head = 2*h*np.cos(np.radians(ic))/v1 + x/v2 if x >= x_crit else float('inf')
    first = min(t_dir, t_ref, t_head)
    label = "Direct" if first == t_dir else ("Refracted" if first == t_head else "Reflected")
    t_head_str = f"{t_head:.2f}" if t_head < 1e6 else "---"
    print(f"  {x:14} {t_dir:12.2f} {t_ref:14.2f} {t_head_str:14} {label}")

# ============================================================
# 3. Earthquake Location (S-P Method)
# ============================================================
print()
print()
print("3. Earthquake Location from S-P Times")
print("-" * 55)

# Synthetic earthquake
eq_x, eq_y, eq_z = 50.0, 30.0, 10.0  # km
eq_t0 = 0.0  # origin time
vp_avg = 6.0
vs_avg = 3.5

stations = [
    ("STA1", 0, 0),
    ("STA2", 100, 0),
    ("STA3", 50, 80),
    ("STA4", 0, 60),
    ("STA5", 80, 40),
]

print(f"  True epicenter: ({eq_x}, {eq_y}) km, depth = {eq_z} km")
print()
print(f"  {'Station':10} {'Dist (km)':12} {'P-time (s)':12} {'S-time (s)':12} {'S-P (s)':10} {'Est dist (km)'}")
print("  " + "-" * 65)

for name, sx, sy in stations:
    dist = np.sqrt((eq_x - sx)**2 + (eq_y - sy)**2 + eq_z**2)
    tp = eq_t0 + dist / vp_avg
    ts = eq_t0 + dist / vs_avg
    sp = ts - tp
    est_dist = sp / (1/vs_avg - 1/vp_avg)
    print(f"  {name:8} {dist:10.2f} {tp:10.2f} {ts:10.2f} {sp:8.2f} {est_dist:10.2f}")

# ============================================================
# 4. Magnitude and Seismic Moment
# ============================================================
print()
print()
print("4. Earthquake Magnitude and Seismic Moment")
print("-" * 55)
print(f"  {'Mw':6} {'M0 (N*m)':16} {'Fault L (km)':14} {'Slip (m)':10} {'Energy (J)'}")
print("  " + "-" * 58)

mu_rock = 3e10  # Pa
for mw in [3, 4, 5, 6, 7, 8, 9]:
    m0 = 10**(1.5*mw + 9.1)
    # Assuming square fault, slip ~ L/20000
    area = (m0 / mu_rock / 0.5)  # D ~ L/20000, A = L^2
    L = np.sqrt(area) / 1000  # km
    D = m0 / (mu_rock * area)
    energy = m0 / 20000
    print(f"  {mw:4} {m0:14.2e} {L:12.1f} {D:8.2f} {energy:12.2e}")

# ============================================================
# 5. Rayleigh Wave Velocity
# ============================================================
print()
print()
print("5. Rayleigh Wave Velocity vs Poisson Ratio")
print("-" * 55)
print(f"  {'nu':10} {'Vp/Vs':10} {'Vr/Vs':10}")
print("  " + "-" * 30)

for nu in [0.10, 0.15, 0.20, 0.25, 0.30, 0.35, 0.40, 0.45, 0.49]:
    vp_vs = np.sqrt(2*(1-nu)/(1-2*nu)) if nu < 0.5 else float('inf')
    # Solve Rayleigh equation numerically
    # x^6 - 8x^4 + (24 - 16*gamma^2)*x^2 - 16*(1-gamma^2) = 0
    # where x = Vr/Vs, gamma = Vs/Vp
    gamma2 = (1-2*nu)/(2*(1-nu))
    coeffs = [1, 0, -8, 0, 24-16*gamma2, 0, -16*(1-gamma2)]
    roots = np.roots(coeffs)
    # Find real positive root less than 1
    real_roots = []
    for r in roots:
        if abs(r.imag) < 1e-10 and r.real > 0 and r.real < 1.01:
            real_roots.append(r.real)
    vr_vs = max(real_roots) if real_roots else 0.92
    print(f"  {nu:8.2f} {vp_vs:8.3f} {vr_vs:8.4f}")

# ============================================================
# PLOTS
# ============================================================
fig, axes = plt.subplots(2, 2, figsize=(14, 11))

# Plot 1: Travel-time curves
ax = axes[0, 0]
distances = np.linspace(0.1, 350, 500)
t_direct = distances / v1
t_reflected = np.sqrt(distances**2 + (2*h)**2) / v1
t_refracted = np.where(distances >= x_crit,
    2*h*np.cos(np.radians(ic))/v1 + distances/v2,
    np.nan)

ax.plot(distances, t_direct, 'b-', linewidth=2, label=f'Direct (V={v1} km/s)')
ax.plot(distances, t_reflected, 'r-', linewidth=2, label='Reflected')
ax.plot(distances, t_refracted, 'g-', linewidth=2, label=f'Head wave (V={v2} km/s)')
ax.axvline(x=x_crit, color='gray', linestyle=':', alpha=0.5, label=f'x_crit={x_crit:.0f} km')
ax.axvline(x=x_cross, color='purple', linestyle=':', alpha=0.5, label=f'x_cross={x_cross:.0f} km')
ax.set_xlabel('Distance (km)')
ax.set_ylabel('Travel Time (s)')
ax.set_title('Travel-Time Curves: Two-Layer Model')
ax.legend(fontsize=8)
ax.set_xlim(0, 350)
ax.grid(True, alpha=0.3)

# Plot 2: Seismic wave simulation (1D)
ax = axes[0, 1]
x = np.linspace(0, 100, 1000)
t_snapshot = 8.0  # seconds
vp_sim, vs_sim = 6.0, 3.5

# P-wave pulse
p_center = vp_sim * t_snapshot
s_center = vs_sim * t_snapshot
width = 3.0

p_wave = np.exp(-((x - p_center)/width)**2) * np.cos(2*np.pi*(x - p_center)/8)
s_wave = 0.7 * np.exp(-((x - s_center)/width)**2) * np.cos(2*np.pi*(x - s_center)/6)

ax.plot(x, p_wave + 1.5, 'b-', linewidth=2, label='P-wave (compressional)')
ax.plot(x, s_wave, 'r-', linewidth=2, label='S-wave (shear)')
ax.fill_between(x, p_wave + 1.5, 1.5, alpha=0.2, color='blue')
ax.fill_between(x, s_wave, 0, alpha=0.2, color='red')
ax.axvline(x=p_center, color='blue', linestyle=':', alpha=0.5)
ax.axvline(x=s_center, color='red', linestyle=':', alpha=0.5)
ax.set_xlabel('Distance (km)')
ax.set_ylabel('Amplitude')
ax.set_title(f'P and S Waves at t = {t_snapshot} s')
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)

# Plot 3: Earthquake magnitude vs energy/moment
ax = axes[1, 0]
mw_range = np.linspace(1, 10, 100)
m0_range = 10**(1.5*mw_range + 9.1)
energy_range = m0_range / 20000

ax.semilogy(mw_range, m0_range, 'b-', linewidth=2, label='Seismic Moment (N*m)')
ax.semilogy(mw_range, energy_range, 'r-', linewidth=2, label='Energy (J)')

# Mark notable earthquakes
notable = [(5.0, 'M5'), (6.0, 'M6'), (7.0, 'Kobe 1995'), (8.0, 'M8'),
           (9.0, 'Tohoku 2011'), (9.5, 'Chile 1960')]
for mw_n, label in notable:
    m0_n = 10**(1.5*mw_n + 9.1)
    ax.plot(mw_n, m0_n, 'ko', markersize=6)
    ax.annotate(label, (mw_n, m0_n), textcoords="offset points",
                xytext=(5, 5), fontsize=7)

ax.set_xlabel('Moment Magnitude (Mw)')
ax.set_ylabel('Value')
ax.set_title('Seismic Moment and Energy vs Magnitude')
ax.legend()
ax.grid(True, alpha=0.3, which='both')

# Plot 4: Rayleigh wave depth decay
ax = axes[1, 1]
nu_plot = 0.25
vp_vs_plot = np.sqrt(3)
vr_vs_plot = 0.9194
wavelength = 50  # km
k = 2*np.pi/wavelength

gamma2_p = 1/3.0
q_alpha = np.sqrt(1 - vr_vs_plot**2 / vp_vs_plot**2)
q_beta = np.sqrt(1 - vr_vs_plot**2)

z = np.linspace(0, 2*wavelength, 200)
kz = k * z

# Normalized displacement components
ux = np.exp(-q_alpha * kz) - 0.5772 * np.exp(-q_beta * kz)
uz = -q_alpha * np.exp(-q_alpha * kz) + 1.4679 * np.exp(-q_beta * kz)

ax.plot(ux/max(abs(ux)), z, 'b-', linewidth=2, label='Horizontal (ux)')
ax.plot(uz/max(abs(uz)), z, 'r-', linewidth=2, label='Vertical (uz)')
ax.axvline(x=0, color='gray', linewidth=0.5)
ax.axhline(y=wavelength, color='gray', linestyle=':', alpha=0.5, label=f'1 wavelength ({wavelength} km)')
ax.set_xlabel('Normalized Amplitude')
ax.set_ylabel('Depth (km)')
ax.set_title('Rayleigh Wave Amplitude vs Depth')
ax.invert_yaxis()
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('output.png', dpi=150, bbox_inches='tight')
plt.close()
print()
print()
print("6. Focal Mechanism Parameters")
print("-" * 55)

print(f"  {'Fault Type':18} {'Strike':10} {'Dip':10} {'Rake':10} {'P-axis':12} {'T-axis'}")
print("  " + "-" * 60)

mechanisms = [
    ("Normal", "N-S", 60, -90, "Vertical", "Horizontal"),
    ("Reverse", "N-S", 30, 90, "Horizontal", "Vertical"),
    ("Left-lateral SS", "N-S", 90, 0, "NE", "NW"),
    ("Right-lateral SS", "N-S", 90, 180, "NW", "NE"),
    ("Oblique normal", "NE-SW", 55, -120, "NE-down", "SW-horiz"),
    ("Oblique thrust", "E-W", 35, 60, "N-horiz", "Down"),
]

for name, strike, dip, rake, p_ax, t_ax in mechanisms:
    print(f"  {name:16} {strike:8} {dip:8} {rake:8} {p_ax:10} {t_ax}")

# Source spectrum
print()
print()
print("7. Earthquake Source Spectrum: Corner Frequency vs Magnitude")
print("-" * 60)

Vs_source = 3500  # m/s
print(f"  {'Mw':6} {'M0 (N.m)':16} {'r0 (km)':12} {'fc (Hz)':12} {'Stress drop (MPa)'}")
print("  " + "-" * 55)

for mw in [2, 3, 4, 5, 6, 7, 8]:
    m0 = 10**(1.5*mw + 9.1)
    # Stress drop ~ 3 MPa (typical)
    stress_drop = 3e6  # Pa
    r0 = (7 * m0 / (16 * stress_drop))**(1.0/3.0) / 1000  # km
    fc = 0.37 * Vs_source / (r0 * 1000)
    actual_sd = 7 * m0 / (16 * (r0*1000)**3) / 1e6
    print(f"  {mw:4} {m0:14.2e} {r0:10.2f} {fc:10.3f} {actual_sd:10.1f}")

print()
print()
print("[Chart saved: travel times, wave pulses, magnitude-energy, Rayleigh decay]")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Earthquake Scaling Relations and Site Effects

Python

script.py69 lines

import numpy as np

print("EARTHQUAKE SCALING AND GROUND MOTION")
print("=" * 70)

print()
print("1. Earthquake Scaling Relations")
print("-" * 65)

print(f"  {'Mw':6} {'M0 (N.m)':14} {'Fault L (km)':14} {'W (km)':10} {'Slip (m)':12} {'Duration (s)'}")
print("  " + "-" * 65)

mu = 3e10  # rigidity (Pa)

for mw in [4, 5, 6, 6.5, 7, 7.5, 8, 8.5, 9, 9.5]:
    m0 = 10**(1.5*mw + 9.1)
    # Wells & Coppersmith (1994) scaling
    log_L = 0.74*mw - 3.55  # subsurface length
    log_W = 0.54*mw - 2.86  # down-dip width
    L = 10**log_L  # km
    W = min(10**log_W, 250)  # km, capped by seismogenic depth
    A = L * W * 1e6  # m^2
    D = m0 / (mu * A)  # m
    duration = L / 3.0  # rough: rupture velocity ~3 km/s
    print(f"  {mw:4.1f} {m0:12.2e} {L:12.1f} {W:8.1f} {D:10.2f} {duration:10.1f}")

print()
print()
print("2. Ground Motion Attenuation (simplified GMPE)")
print("-" * 65)

# Simplified: ln(PGA) = -1.5 + 1.5*M - 1.5*ln(R) - 0.003*R
print(f"  PGA in g units (9.81 m/s^2):")
print(f"  {'R (km)':10}", end="")
for M in [5, 6, 7, 8, 9]:
    print(f"{'M' + str(M):12}", end="")
print()
print("  " + "-" * 60)

for R in [5, 10, 20, 50, 100, 200, 500]:
    print(f"  {R:8}", end="")
    for M in [5, 6, 7, 8, 9]:
        ln_pga = -1.5 + 1.5*M - 1.5*np.log(R) - 0.003*R
        pga_g = np.exp(ln_pga) / 981  # cm/s^2 to g
        print(f"  {pga_g:10.4f}", end="")
    print(" g")

# Site amplification
print()
print()
print("3. Site Amplification Factors (NEHRP classes)")
print("-" * 55)

site_classes = [
    ("A - Hard rock", ">1500", 0.8),
    ("B - Rock", "760-1500", 1.0),
    ("C - Soft rock", "360-760", 1.2),
    ("D - Stiff soil", "180-360", 1.6),
    ("E - Soft soil", "<180", 2.5),
]

print(f"  {'Site class':22} {'Vs30 (m/s)':14} {'Amplification'}")
print("  " + "-" * 45)
for name, vs30, amp in site_classes:
    print(f"  {name:20} {vs30:12} {amp:8.1f}x")

print()
print("Complete: ground motion analysis finished.")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Practice Problems

Problem 1:A P-wave travels through a uniform layer of crustal rock with velocity $v_P = 6.5\;\text{km/s}$. A station records the direct P-wave arrival 12.3 s after the origin time. Calculate the epicentral distance. If a reflected wave from a layer at 35 km depth arrives at the same station, calculate its travel time.

Solution:

1. For direct P-wave in a uniform layer: $d = v_P \times t = 6.5 \times 12.3 = 79.95\;\text{km}$.

2. Epicentral distance: $\Delta \approx 80.0\;\text{km}$ (assuming shallow source).

3. Reflected wave path: down to 35 km reflector and back up. Total path length: $L = 2\sqrt{(\Delta/2)^2 + h^2} = 2\sqrt{(40)^2 + (35)^2}$.

4. $L = 2\sqrt{1600 + 1225} = 2\sqrt{2825} = 2 \times 53.15 = 106.3\;\text{km}$.

5. Reflected wave travel time: $t_{\text{refl}} = L/v_P = 106.3/6.5 = 16.35\;\text{s}$.

6. Time difference: $\Delta t = 16.35 - 12.3 = 4.05\;\text{s}$. This delay constrains the depth of the reflecting interface (Moho in this case), forming the basis of seismic reflection profiling.

Problem 2:Apply Snell's law at the core-mantle boundary (CMB). A P-wave in the lower mantle has velocity $v_1 = 13.7\;\text{km/s}$ and strikes the CMB at incidence angle $i_1 = 25°$. The P-wave velocity in the outer core is $v_2 = 8.1\;\text{km/s}$. Find the refraction angle and the ray parameter.

Solution:

1. Snell's law: $\frac{\sin i_1}{v_1} = \frac{\sin i_2}{v_2}$.

2. $\sin i_2 = \frac{v_2}{v_1}\sin i_1 = \frac{8.1}{13.7}\sin 25° = 0.5912 \times 0.4226 = 0.2498$.

3. $i_2 = \arcsin(0.2498) = 14.47°$.

4. The ray bends toward the normal upon entering the core (lower velocity), creating the P-wave shadow zone.

5. Ray parameter: $p = \frac{\sin i_1}{v_1} = \frac{0.4226}{13.7} = 0.03084\;\text{s/km}$. For a spherical Earth: $p = \frac{r\sin i}{v}$.

6. Critical angle for total reflection: $i_c = \arcsin(v_1/v_2)$ -- but here $v_2 < v_1$, so no critical angle exists for the transmitted wave. Instead, the velocity decrease causes the shadow zone at epicentral distances 104-140°.

Problem 3:A seismometer records a maximum displacement amplitude of $A = 23\;\mu\text{m}$ for a seismic wave with period $T = 1.0\;\text{s}$ at an epicentral distance of $\Delta = 100\;\text{km}$. Calculate the local (Richter) magnitude $M_L$.

Solution:

1. Richter magnitude definition: $M_L = \log_{10}(A) - \log_{10}(A_0(\Delta))$, where $A$ is in micrometers (on a Wood-Anderson seismograph).

2. At $\Delta = 100\;\text{km}$: $-\log_{10}(A_0) = 3.0$ (standard Richter calibration).

3. $M_L = \log_{10}(23) + 3.0 = 1.362 + 3.0 = 4.36$.

4. Rounding: $M_L \approx 4.4$. This is a light earthquake, typically felt but causing minor damage.

5. Seismic moment: using the empirical relation $\log_{10} M_0 = 1.5 M_L + 16.1$ (cgs): $\log_{10} M_0 = 1.5(4.4) + 16.1 = 22.7$, so $M_0 \approx 5 \times 10^{22}\;\text{dyne·cm} = 5 \times 10^{15}\;\text{N·m}$.

6. Note: $M_L$ saturates above ~6.5 due to the frequency response of the Wood-Anderson instrument. For large earthquakes, moment magnitude $M_w = \frac{2}{3}\log_{10}M_0 - 10.7$ (SI) is preferred.

Problem 4:Three seismic stations record P-wave arrival times from the same earthquake. Station A ($x_A = 0, y_A = 0$): $t_A = 10.0\;\text{s}$; Station B ($x_B = 80, y_B = 0$ km): $t_B = 15.0\;\text{s}$; Station C ($x_C = 0, y_C = 60$ km): $t_C = 13.0\;\text{s}$. Locate the epicenter assuming $v_P = 6.0\;\text{km/s}$.

Solution:

1. Distances from each station: $d_i = v_P \times t_i$. $d_A = 60\;\text{km}$, $d_B = 90\;\text{km}$, $d_C = 78\;\text{km}$.

2. Circle equations: $x^2 + y^2 = 60^2 = 3600$ ... (i); $(x-80)^2 + y^2 = 90^2 = 8100$ ... (ii); $x^2 + (y-60)^2 = 78^2 = 6084$ ... (iii).

3. Subtract (i) from (ii): $-160x + 6400 = 4500$, so $x = \frac{1900}{160} = 11.875\;\text{km}$.

4. Subtract (i) from (iii): $-120y + 3600 = 2484$, so $y = \frac{1116}{120} = 9.3\;\text{km}$.

5. Verify with (i): $11.875^2 + 9.3^2 = 141.0 + 86.5 = 227.5$. But $d_A^2 = 3600$. Discrepancy indicates the origin time $t_0$ is not zero -- we must solve for $t_0$ as well.

6. Including $t_0$: use $(x-x_i)^2 + (y-y_i)^2 = v_P^2(t_i - t_0)^2$. This gives 4 unknowns ($x, y, z, t_0$); with 3 stations and assumed surface focus, the system is exactly determined. The epicenter is at approximately ($x, y$) = (12, 9) km relative to Station A.

Problem 5:An earthquake occurs at 10 km depth. The P-wave velocity is 6.0 km/s and S-wave velocity is 3.5 km/s. A station records the S-P time difference as 8.5 s. Calculate the epicentral distance and the hypocentral distance.

Solution:

1. S-P time: $\Delta t_{S-P} = d\left(\frac{1}{v_S} - \frac{1}{v_P}\right)$, where $d$ is the hypocentral distance.

2. $\frac{1}{v_S} - \frac{1}{v_P} = \frac{1}{3.5} - \frac{1}{6.0} = 0.2857 - 0.1667 = 0.1190\;\text{s/km}$.

3. Hypocentral distance: $d = \frac{\Delta t_{S-P}}{1/v_S - 1/v_P} = \frac{8.5}{0.1190} = 71.4\;\text{km}$.

4. Epicentral distance: $\Delta = \sqrt{d^2 - h^2} = \sqrt{71.4^2 - 10^2} = \sqrt{5098 - 100} = \sqrt{4998} = 70.7\;\text{km}$.

5. P-wave travel time: $t_P = d/v_P = 71.4/6.0 = 11.9\;\text{s}$.

6. The S-P method is the classic technique for rapid distance estimation. The Wadati diagram ($t_S - t_P$ vs. $t_P$) gives $v_P/v_S$ from the slope $(v_P/v_S - 1)$, which equals the Poisson's ratio indicator: here $v_P/v_S = 1.71$, consistent with $\sigma \approx 0.24$ (typical crustal rock).

← Structural Geology Next: Earth's Interior →

Share:X Reddit LinkedIn