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)$.

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

← Structural Geology Next: Earth's Interior →