Part IV: Applied Earth Science | Chapter 4

Natural Hazards

Earthquake statistics, volcanic hazards, landslides, and tsunami physics

4.1 Earthquake Hazard Assessment

Derivation 1: The Gutenberg-Richter Law

Beno Gutenberg and Charles Richter (1944) discovered that earthquake frequency follows a power-law distribution. The number of earthquakes with magnitude $\geq M$ in a given region and time period is:

$$\log_{10} N(M) = a - bM$$

where $a$ describes the overall seismicity rate and $b \approx 1.0$ globally (varying locally from ~0.7 to ~1.3). The $b$-value means that for each unit increase in magnitude, earthquakes become ~10 times less frequent. In terms of seismic moment:

$$N(M_0) \propto M_0^{-2b/3}$$

This power law reflects the self-similar (fractal) nature of faulting. The$b$-value is related to the stress state: lower $b$ indicates higher differential stress (and potentially larger maximum earthquakes). Volcanic regions often have $b > 1$ due to fluid-driven swarms of small earthquakes.

Derivation 2: Seismic Hazard and Return Periods

The Poisson model assumes earthquakes are independent in time. The probability of at least one event of magnitude $\geq M$ in time $t$ is:

$$P(\geq 1) = 1 - e^{-\nu(M) t}$$

where $\nu(M) = 10^{a-bM}$ is the annual rate. The return period is:

$$T_R = \frac{1}{\nu(M)} = 10^{bM - a}$$

For example, if a region has $a = 4$ and $b = 1$ (events per year with $M \geq 0$), then an M7 earthquake occurs on average once every $10^{7-4} = 1000$ years, with a ~5% chance in any 50-year window.

Historical Context

Probabilistic seismic hazard analysis (PSHA) was formalized by Allin Cornell (1968). It combines the G-R law with ground motion prediction equations (GMPEs) to estimate the probability of exceeding various levels of ground shaking at a site. PSHA underpins building codes worldwide and is critical for nuclear power plant siting and dam safety assessment.

4.2 Volcanic Hazards

Derivation 3: Eruption Column Height and Mass Flux

The height of a volcanic eruption column determines the extent of ash dispersal. For a buoyant plume rising through a stably stratified atmosphere, Morton, Taylor, and Turner (1956) showed that the maximum height scales as:

$$H = k \dot{Q}^{1/4}$$

where $\dot{Q}$ is the thermal flux at the vent and $k$ depends on atmospheric stratification. For volcanic plumes, Sparks (1986) derived:

$$H \approx 0.236 \dot{m}^{0.259}$$

where $H$ is in km and $\dot{m}$ is mass eruption rate in kg/s. The Volcanic Explosivity Index (VEI) classifies eruptions on a logarithmic scale of ejecta volume:

VEI 0-2: Effusive to small explosive ($< 10^7$ m$^3$)
VEI 3-4: Moderate explosive ($10^7$-$10^9$ m$^3$), e.g., Eyjafjallajokull 2010
VEI 5-6: Large explosive ($10^9$-$10^{11}$ m$^3$), e.g., Mt. St. Helens 1980, Pinatubo 1991
VEI 7-8: Supereruption ($> 10^{11}$ m$^3$), e.g., Tambora 1815, Toba ~74 ka

Pyroclastic density currents (PDCs) are the most lethal volcanic hazard, traveling at 100-700 km/hr at temperatures of 200-700$°$C. Their runout distance scales with collapse height and volume.

4.3 Landslide Mechanics

Derivation 4: The Infinite Slope Stability Analysis

The simplest slope stability analysis considers an infinite planar slope of angle $\beta$with a potential slip surface at depth $z$. The factor of safety (FS) is the ratio of resisting to driving forces:

$$FS = \frac{\tau_{\text{resist}}}{\tau_{\text{drive}}} = \frac{c + (\sigma_n - u)\tan\phi}{\tau}$$

For a soil block of unit area on a slope, the normal and shear stresses on the failure plane are:

$$\sigma_n = \gamma z \cos^2\beta, \quad \tau = \gamma z \sin\beta \cos\beta$$

With pore water pressure $u = \gamma_w z_w \cos^2\beta$ where $z_w$ is the height of water table above the slip surface:

$$FS = \frac{c}{\gamma z \sin\beta \cos\beta} + \frac{(\gamma - m\gamma_w)\tan\phi}{\gamma \tan\beta}$$

where $m = z_w/z$ is the ratio of water table height to slip surface depth. Failure occurs when $FS < 1$. For cohesionless soil ($c = 0$) on a dry slope ($m = 0$), failure occurs when $\beta > \phi$ (the angle of repose). Rainfall-induced landslides typically occur when rising water tables reduce the effective stress.

4.4 Tsunami Physics

Derivation 5: Shallow Water Wave Speed and Amplification

Tsunamis are long-wavelength gravity waves generated by seafloor displacement during earthquakes, volcanic collapses, or submarine landslides. Since their wavelength (~100-500 km) far exceeds ocean depth ($\leq$ 4 km), they behave as shallow water waves with phase speed:

$$c = \sqrt{gh}$$

where $h$ is the water depth. In the open ocean with $h \approx 4000$ m:$c \approx \sqrt{9.81 \times 4000} \approx 200$ m/s $\approx$ 720 km/hr, comparable to a jet aircraft. As the tsunami approaches shore and water shallows, the wave slows and its amplitude grows. Conservation of energy flux gives:

$$E c = \frac{1}{2}\rho g A^2 \sqrt{gh} = \text{const}$$

Therefore, the amplitude scales as:

$$\frac{A_2}{A_1} = \left(\frac{h_1}{h_2}\right)^{1/4}$$

This is Green's law. A 1 m wave in 4000 m depth amplifies to ~4.5 m in 10 m depth: $(4000/10)^{1/4} = 4.47$. Additional amplification occurs from focusing, resonance in bays, and nonlinear effects, leading to runup heights of 10-40 m in extreme cases.

The 2004 Indian Ocean tsunami ($M_w$ 9.1) generated waves up to 30 m high and killed ~230,000 people. The 2011 Tohoku tsunami reached 40 m runup height locally. These disasters led to expanded warning systems: the Pacific Tsunami Warning Center (est. 1949) and the Indian Ocean Tsunami Warning System (est. 2006).

4.5 Multi-Hazard Risk Assessment

Risk is the product of hazard, exposure, and vulnerability:

$$\text{Risk} = \text{Hazard} \times \text{Exposure} \times \text{Vulnerability}$$

Hazard is the probability of a given intensity of the natural phenomenon occurring. Exposure is the number of people and value of assets in the affected area. Vulnerability is the fraction of exposure that would be lost. Modern risk assessment uses probabilistic methods integrating over all possible event magnitudes and locations, combined with fragility curves for structures and population exposure models.

Cascading hazards are particularly dangerous: earthquakes trigger tsunamis and landslides, volcanic eruptions produce lahars and PDCs, and climate change amplifies flood and landslide frequency. The Fukushima disaster (2011) exemplified cascading risk: earthquake damage to infrastructure compounded by tsunami inundation leading to nuclear meltdown.

Flood Hazard and Return Periods

River floods follow similar statistical frameworks to earthquakes. The annual exceedance probability for a flood of magnitude $Q$ is estimated from the frequency analysis of historical records, often using the Gumbel (extreme value type I) distribution:

$$P(Q > q) = 1 - \exp\left(-\exp\left(-\frac{q - \mu}{\sigma}\right)\right)$$

The "100-year flood" has a 1% annual exceedance probability, meaning a 26% chance of occurring within a 30-year mortgage period. Climate change is altering flood statistics: what was historically a 100-year event may become a 25-50 year event in many regions due to intensified precipitation from increased atmospheric moisture content (Clausius-Clapeyron scaling: ~7% more moisture per degree of warming).

Induced Seismicity

Human activities can trigger earthquakes through fluid injection (wastewater disposal, hydraulic fracturing, enhanced geothermal systems), reservoir impoundment, and mining. The 2011 Oklahoma earthquake sequence (up to M5.7) was linked to deep wastewater injection. The mechanism is pore pressure increase reducing effective stress on pre-existing faults, consistent with the Mohr-Coulomb analysis. Managing induced seismicity through traffic-light protocols (monitoring and adjusting injection rates) is an active area of applied geophysics.

4.6 Ground Motion Prediction and Attenuation

Ground motion prediction equations (GMPEs) relate earthquake magnitude and distance to the expected intensity of ground shaking. A typical GMPE for peak ground acceleration (PGA) takes the form:

$$\ln(\text{PGA}) = c_1 + c_2 M + c_3 \ln(R + c_4 e^{c_5 M}) + c_6 S$$

where $R$ is distance, $M$ is magnitude, and $S$ is a site class term. The logarithmic distance dependence reflects geometric spreading ($\sim 1/R$ for body waves) and anelastic attenuation. Site effects amplify shaking on soft soils (basin amplification), as demonstrated tragically in Mexico City (1985, M8.0) where the lake sediments amplified shaking by factors of 5-10.

The Modified Mercalli Intensity (MMI) scale describes observed effects from I (not felt) to XII (total destruction). The empirical relationship between PGA and MMI is:

$$\text{MMI} \approx 3.66 \log_{10}(\text{PGA}) + 1.66$$

where PGA is in cm/s$^2$. Building codes specify design ground motions based on PSHA: the typical design basis in the USA is the 2% probability of exceedance in 50 years (return period ~2475 years). Japan, with higher seismicity, designs for even rarer events following the Tohoku earthquake experience.

Earthquake Early Warning

Early warning systems exploit the speed difference between electromagnetic signals (speed of light) and seismic waves (~6-8 km/s for P-waves). Detection of P-waves near the source can provide seconds to tens of seconds of warning before destructive S-waves and surface waves arrive at distant cities. Japan's system issued a warning 8 seconds before the Tohoku earthquake S-waves reached Tokyo, allowing bullet trains to brake and factories to shut down. The ShakeAlert system in the western USA began public alerting in 2019.

Computational Lab: Natural Hazards

Gutenberg-Richter Law, Tsunami Propagation, and Slope Stability

Python

script.py349 lines

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

# ============================================================
# NATURAL HAZARDS SIMULATION
# ============================================================
print("NATURAL HAZARDS SIMULATION")
print("=" * 70)

# ============================================================
# 1. Gutenberg-Richter Frequency-Magnitude Relation
# ============================================================
print()
print("1. Gutenberg-Richter Law: Earthquake Statistics")
print("-" * 65)

a_val = 5.0  # log10(N) for M>=0
b_val = 1.0

print(f"  Parameters: a = {a_val}, b = {b_val}")
print(f"  log10(N) = {a_val} - {b_val}*M")
print()
print(f"  {'Magnitude':12} {'N/year':14} {'Return (yr)':14} {'P(50yr) %':12} {'Energy (J)'}")
print("  " + "-" * 60)

for M in range(0, 10):
    N = 10**(a_val - b_val * M)
    T_ret = 1.0 / N if N > 0 else float('inf')
    P_50 = (1 - np.exp(-N * 50)) * 100
    energy = 10**(1.5 * M + 4.8)
    print(f"  {M:10} {N:12.2f} {T_ret:12.1f} {P_50:10.1f} {energy:12.2e}")

# Cumulative energy
print()
total_E = 0
for M in range(0, 10):
    N = 10**(a_val - b_val * M)
    energy = 10**(1.5 * M + 4.8)
    cum_E = N * energy
    total_E += cum_E
    print(f"  M{M}: {N:.1f} events x {energy:.2e} J = {cum_E:.2e} J/yr ({cum_E/total_E*100:.1f}% cumulative)")

# ============================================================
# 2. Seismic Hazard: Different b-values
# ============================================================
print()
print()
print("2. Effect of b-value on Seismic Hazard")
print("-" * 55)

print(f"  {'b-value':10} {'N(M>=5)/yr':14} {'N(M>=7)/yr':14} {'T_ret(M7) yr':14} {'Max M (50yr)'}")
print("  " + "-" * 55)

for b in [0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3]:
    N5 = 10**(a_val - b * 5)
    N7 = 10**(a_val - b * 7)
    T7 = 1.0 / N7
    # Max M with return period <= 50 yr
    M_max = (a_val - np.log10(1.0/50)) / b
    print(f"  {b:8.1f} {N5:12.2f} {N7:12.4f} {T7:12.0f} {M_max:10.1f}")

# ============================================================
# 3. Tsunami Propagation
# ============================================================
print()
print()
print("3. Tsunami Wave Speed and Amplification")
print("-" * 65)

g = 9.81

print(f"  {'Depth (m)':14} {'Speed (m/s)':14} {'Speed (km/hr)':16} {'Amp factor':14} {'A for 1m deep'}")
print("  " + "-" * 65)

h_ref = 4000
A_ref = 1.0  # 1 m in deep ocean

for h in [5000, 4000, 2000, 1000, 500, 200, 100, 50, 20, 10, 5, 2]:
    c = np.sqrt(g * h)
    c_kmhr = c * 3.6
    amp_factor = (h_ref / h)**0.25
    A_shore = A_ref * amp_factor
    print(f"  {h:12} {c:12.1f} {c_kmhr:14.0f} {amp_factor:12.2f} {A_shore:10.2f}")

# ============================================================
# 4. Slope Stability Analysis
# ============================================================
print()
print()
print("4. Infinite Slope Stability: Factor of Safety")
print("-" * 65)

gamma = 18.0  # kN/m^3 (soil unit weight)
gamma_w = 9.81  # kN/m^3 (water)
phi = 30.0  # friction angle (degrees)
c = 5.0  # cohesion (kPa)
z = 3.0  # slip depth (m)

print(f"  gamma = {gamma} kN/m^3, phi = {phi} deg, c = {c} kPa, z = {z} m")
print()
print(f"  {'Slope (deg)':14} {'m=0 (dry)':12} {'m=0.5':12} {'m=1.0 (sat)':14} {'Critical m'}")
print("  " + "-" * 55)

for beta in [10, 15, 20, 25, 30, 35, 40, 45]:
    beta_rad = np.radians(beta)
    results = []
    for m in [0, 0.5, 1.0]:
        FS = c / (gamma * z * np.sin(beta_rad) * np.cos(beta_rad)) +              (gamma - m * gamma_w) * np.tan(np.radians(phi)) / (gamma * np.tan(beta_rad))
        results.append(FS)

# Find critical m (FS = 1)
    # FS = c/(gamma*z*sin*cos) + (gamma - m*gamma_w)*tan(phi)/(gamma*tan(beta))
    # 1 = c_term + (gamma - m*gamma_w)*tan_phi_term
    c_term = c / (gamma * z * np.sin(beta_rad) * np.cos(beta_rad))
    tan_term = np.tan(np.radians(phi)) / (gamma * np.tan(beta_rad))
    # 1 - c_term = (gamma - m*gamma_w) * tan_term
    # (1 - c_term)/tan_term = gamma - m*gamma_w
    if tan_term > 0:
        m_crit = (gamma - (1 - c_term) / tan_term) / gamma_w
        m_crit = np.clip(m_crit, 0, 10)
    else:
        m_crit = 0
    m_str = f"{m_crit:.2f}" if m_crit < 5 else ">5"

print(f"  {beta:12} {results[0]:10.2f} {results[1]:10.2f} {results[2]:12.2f} {m_str}")

# ============================================================
# 5. Volcanic Eruption Column Height
# ============================================================
print()
print()
print("5. Volcanic Eruption Column Height vs Mass Flux")
print("-" * 60)

print(f"  {'Mass rate (kg/s)':18} {'H (km)':12} {'VEI':8} {'Example'}")
print("  " + "-" * 55)

eruptions = [
    (1e3, "Stromboli-like"),
    (1e4, "Vulcanian"),
    (1e5, "Sub-Plinian"),
    (1e6, "St. Helens 1980"),
    (1e7, "Plinian (Pinatubo)"),
    (1e8, "Ultra-Plinian"),
    (1e9, "Supervolcano"),
]

for mdot, example in eruptions:
    H = 0.236 * mdot**0.259  # Sparks formula
    if H < 1:
        vei = 0
    elif H < 5:
        vei = 2
    elif H < 15:
        vei = 3
    elif H < 25:
        vei = 4
    elif H < 35:
        vei = 5
    else:
        vei = 6
    print(f"  {mdot:16.0e} {H:10.1f} {vei:6} {example}")

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

# Plot 1: Gutenberg-Richter law
ax = axes[0, 0]
M_range = np.linspace(0, 9, 100)

for b, color in [(0.7, 'red'), (1.0, 'blue'), (1.3, 'green')]:
    N_range = 10**(a_val - b * M_range)
    ax.semilogy(M_range, N_range, color=color, linewidth=2, label=f'b = {b}')

ax.set_xlabel('Magnitude (M)')
ax.set_ylabel('Cumulative number per year (N >= M)')
ax.set_title('Gutenberg-Richter Frequency-Magnitude')
ax.legend()
ax.grid(True, alpha=0.3, which='both')
ax.set_xlim(0, 9)
ax.set_ylim(1e-5, 1e6)

# Mark notable earthquakes
notable_eq = [(5.0, 'M5 (common)'), (7.0, 'Kobe 1995'),
              (8.8, 'Chile 2010'), (9.1, 'Tohoku 2011')]
for m_n, label in notable_eq:
    n_n = 10**(a_val - 1.0 * m_n)
    ax.plot(m_n, n_n, 'ko', markersize=6)
    ax.annotate(label, (m_n, n_n), textcoords="offset points",
                xytext=(5, 5), fontsize=7)

# Plot 2: Tsunami propagation and amplification
ax = axes[0, 1]
h_range = np.logspace(0.3, 4, 200)  # 2 to 10000 m
c_range = np.sqrt(g * h_range) * 3.6  # km/hr
A_range = (4000 / h_range)**0.25

ax2 = ax.twinx()
ax.semilogx(h_range, c_range, 'b-', linewidth=2, label='Wave speed (km/hr)')
ax2.semilogx(h_range, A_range, 'r-', linewidth=2, label='Amplification factor')

ax.set_xlabel('Water depth (m)')
ax.set_ylabel('Wave speed (km/hr)', color='blue')
ax2.set_ylabel('Amplitude amplification', color='red')
ax.set_title('Tsunami Speed and Amplification')
ax.tick_params(axis='y', labelcolor='blue')
ax2.tick_params(axis='y', labelcolor='red')
ax.legend(loc='upper left', fontsize=9)
ax2.legend(loc='upper right', fontsize=9)
ax.grid(True, alpha=0.3, which='both')
ax.invert_xaxis()

# Plot 3: Slope stability
ax = axes[1, 0]
beta_arr = np.linspace(5, 50, 200)
beta_rad_arr = np.radians(beta_arr)

for m, color, label in [(0, 'green', 'Dry (m=0)'), (0.5, 'orange', 'm=0.5'),
                         (1.0, 'red', 'Saturated (m=1)')]:
    FS_arr = c / (gamma * z * np.sin(beta_rad_arr) * np.cos(beta_rad_arr)) +              (gamma - m * gamma_w) * np.tan(np.radians(phi)) / (gamma * np.tan(beta_rad_arr))
    ax.plot(beta_arr, FS_arr, color=color, linewidth=2, label=label)

ax.axhline(y=1.0, color='black', linewidth=2, linestyle='--', label='FS = 1 (failure)')
ax.axhline(y=1.5, color='gray', linestyle=':', alpha=0.5, label='FS = 1.5 (design)')
ax.fill_between(beta_arr, 0, 1, alpha=0.1, color='red')
ax.set_xlabel('Slope angle (degrees)')
ax.set_ylabel('Factor of Safety')
ax.set_title('Slope Stability vs Slope Angle')
ax.legend(fontsize=8)
ax.grid(True, alpha=0.3)
ax.set_xlim(5, 50)
ax.set_ylim(0, 5)

# Plot 4: Tsunami 1D propagation simulation
ax = axes[1, 1]
# Simple 1D shallow water wave propagation
nx = 1000
x = np.linspace(0, 500, nx)  # km
dx = x[1] - x[0]

# Bathymetry: constant depth then shoaling
h_bathy = np.ones(nx) * 4000
shelf_start = 350  # km
for i in range(nx):
    if x[i] > shelf_start:
        h_bathy[i] = 4000 * (1 - 0.99 * (x[i] - shelf_start) / (500 - shelf_start))
        h_bathy[i] = max(h_bathy[i], 5)

# Initial condition: Gaussian wave
eta0 = 1.0 * np.exp(-((x - 100)**2) / (50**2))  # 1m wave at x=100km

# Apply Green's law for amplitude
eta_green = eta0 * (4000 / h_bathy)**0.25

# Plot several "snapshots" by shifting the wave
for t_label, x_shift in [("t=0", 0), ("t=30 min", 120), ("t=60 min", 240),
                          ("t=75 min", 330), ("t=85 min", 380)]:
    eta_shifted = 1.0 * np.exp(-((x - 100 - x_shift)**2) / (50**2))
    eta_amp = eta_shifted * (4000 / h_bathy)**0.25
    ax.plot(x, eta_amp, linewidth=1.5, label=t_label)

# Bathymetry (secondary axis)
ax_bathy = ax.twinx()
ax_bathy.fill_between(x, -h_bathy/1000, 0, alpha=0.15, color='blue')
ax_bathy.set_ylabel('Depth (km)', color='blue')
ax_bathy.set_ylim(-5, 0)
ax_bathy.tick_params(axis='y', labelcolor='blue')

ax.set_xlabel('Distance (km)')
ax.set_ylabel('Wave amplitude (m)')
ax.set_title('Tsunami Propagation and Shoaling')
ax.legend(fontsize=8, loc='upper left')
ax.grid(True, alpha=0.3)
ax.set_xlim(0, 500)

plt.tight_layout()
plt.savefig('output.png', dpi=150, bbox_inches='tight')
plt.close()
print()
print()
print("6. Ground Motion Attenuation")
print("-" * 60)

print(f"  PGA (cm/s^2) vs Distance for Different Magnitudes:")
print(f"  {'Distance (km)':16}", end="")
for M_gm in [5, 6, 7, 8]:
    print(f"{'M=' + str(M_gm):14}", end="")
print()
print("  " + "-" * 60)

# Simplified GMPE: ln(PGA) = c1 + c2*M + c3*ln(R)
c1_gm, c2_gm, c3_gm = -0.5, 1.5, -1.5

for R_dist in [5, 10, 20, 50, 100, 200, 500]:
    print(f"  {R_dist:14}", end="")
    for M_gm in [5, 6, 7, 8]:
        ln_pga = c1_gm + c2_gm * M_gm + c3_gm * np.log(R_dist)
        pga = np.exp(ln_pga)
        print(f"  {pga:12.1f}", end="")
    print(" cm/s^2")

# MMI conversion
print()
print("  PGA to MMI Conversion:")
print(f"  {'PGA (cm/s^2)':16} {'MMI':8} {'Description'}")
print("  " + "-" * 45)

mmi_descriptions = {
    1: "Not felt", 2: "Weak", 3: "Weak", 4: "Light",
    5: "Moderate", 6: "Strong", 7: "Very strong", 8: "Severe",
    9: "Violent", 10: "Extreme", 11: "Extreme", 12: "Total destruction"
}

for pga_val in [0.5, 1, 2, 5, 10, 20, 50, 100, 200, 500, 1000]:
    mmi_val = 3.66 * np.log10(pga_val) + 1.66
    mmi_int = max(1, min(12, int(round(mmi_val))))
    desc = mmi_descriptions.get(mmi_int, "")
    print(f"  {pga_val:14.0f} {mmi_val:6.1f} {desc}")

# Early warning time calculation
print()
print()
print("7. Earthquake Early Warning: Available Time")
print("-" * 55)

vp_ew = 6.0  # km/s P-wave
vs_ew = 3.5  # km/s S-wave
v_light = 300000  # km/s (effectively instant)

print(f"  {'Distance (km)':16} {'P arrival (s)':16} {'S arrival (s)':16} {'Warning time (s)'}")
print("  " + "-" * 60)

for dist_ew in [10, 20, 30, 50, 75, 100, 150, 200, 300, 500]:
    tp = dist_ew / vp_ew
    ts = dist_ew / vs_ew
    # Assume P detected at epicenter, warning sent at light speed
    detect_time = 3.0  # seconds to detect P and issue warning
    warning = ts - detect_time  # simplified
    warning = max(0, warning)
    print(f"  {dist_ew:14} {tp:14.1f} {ts:14.1f} {warning:10.1f}")

print()
print()
print("[Chart saved: G-R law, tsunami speed, slope stability, tsunami propagation]")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Seismic Hazard Analysis and VEI Statistics

Python

script.py64 lines

import numpy as np

print("PROBABILISTIC SEISMIC HAZARD ANALYSIS")
print("=" * 70)

# ============================================================
# 1. Poisson Probability of Exceedance
# ============================================================
print()
print("1. Earthquake Probability by Return Period and Exposure Time")
print("-" * 65)

print(f"  {'Return period':18}", end="")
for t_exp in [1, 10, 30, 50, 100, 200, 500]:
    print(f"{'t=' + str(t_exp) + 'yr':10}", end="")
print()
print("  " + "-" * 80)

for T_ret in [50, 100, 200, 475, 1000, 2500, 5000, 10000]:
    nu = 1.0 / T_ret
    print(f"  {T_ret:16} yr", end="")
    for t_exp in [1, 10, 30, 50, 100, 200, 500]:
        P = (1 - np.exp(-nu * t_exp)) * 100
        print(f"  {P:8.1f}", end="")
    print(" %")

# ============================================================
# 2. VEI Statistics
# ============================================================
print()
print()
print("2. Volcanic Explosivity Index: Frequency and Impact")
print("-" * 65)

vei_data = [
    (0, "Effusive", 1e4, 1.0, "Kilauea"),
    (1, "Gentle", 1e5, 5, "Stromboli"),
    (2, "Explosive", 1e6, 12, "Galeras 1993"),
    (3, "Severe", 1e7, 3, "Nevado del Ruiz 1985"),
    (4, "Cataclysmic", 1e8, 1, "Eyjafjallajokull 2010"),
    (5, "Paroxysmal", 1e9, 0.17, "Mt St Helens 1980"),
    (6, "Colossal", 1e10, 0.033, "Pinatubo 1991"),
    (7, "Super-colossal", 1e11, 0.002, "Tambora 1815"),
    (8, "Mega-colossal", 1e12, 0.00001, "Yellowstone 640ka"),
]

print(f"  {'VEI':6} {'Type':16} {'Volume (m3)':14} {'Per yr':10} {'Return (yr)':14} {'Example'}")
print("  " + "-" * 75)

for vei, vtype, vol, freq, example in vei_data:
    ret = 1.0/freq if freq > 0 else float('inf')
    ret_str = f"{ret:.0f}" if ret < 1e6 else f"{ret:.0e}"
    print(f"  {vei:4} {vtype:14} {vol:12.0e} {freq:8.3f} {ret_str:12} {example}")

# Climate impact
print()
print("  Volcanic climate impacts:")
print("    VEI 6+: Global cooling 0.1-0.5 C for 1-3 years")
print("    VEI 7 (Tambora): 'Year without summer' 1816, ~0.5 C cooling")
print("    VEI 8 (Toba): possible 3-5 C cooling for decades")

print()
print("Complete: hazard analysis finished.")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

← Hydrogeology Course Overview →