7.5 Tsunamis

Tsunamis are long-wavelength gravity waves generated by impulsive seafloor displacement. With wavelengths of hundreds of kilometers and periods of 10–60 minutes, they behave as shallow water waves even in the deepest ocean, propagating at jet-aircraft speeds and amplifying catastrophically upon reaching coastlines.

Shallow Water Wave Physics

Tsunamis qualify as shallow water waves because their wavelength (λ ~ 100–500 km) vastly exceeds the ocean depth (h ~ 4 km), so kh « 1. The phase speed depends only on depth:

$$c = \sqrt{gh}$$

At h = 4000 m: c ≈ 198 m/s ≈ 713 km/hr (jet aircraft speed)

Derivation: Shallow Water Phase Speed

Starting from the linearized shallow water equations for a fluid of uniform depth $H$ with surface elevation $\eta$ and depth-averaged velocity $u$:

$$\frac{\partial u}{\partial t} = -g\frac{\partial \eta}{\partial x} \quad \text{(momentum equation)}$$

$$\frac{\partial \eta}{\partial t} = -H\frac{\partial u}{\partial x} \quad \text{(continuity equation)}$$

Differentiate the continuity equation with respect to $t$ and substitute the momentum equation:

$$\frac{\partial^2 \eta}{\partial t^2} = -H\frac{\partial}{\partial x}\frac{\partial u}{\partial t} = -H\frac{\partial}{\partial x}\left(-g\frac{\partial \eta}{\partial x}\right) = gH\frac{\partial^2 \eta}{\partial x^2}$$

This is the classical wave equation $\partial^2\eta/\partial t^2 = c^2\,\partial^2\eta/\partial x^2$ with phase speed:

$$c = \sqrt{gH}$$

Because these waves are non-dispersive, all frequency components travel at the same speed, and the wave form is preserved over transoceanic distances. The wave energy is distributed over the entire water column:

$$E = \frac{1}{2}\rho g \eta^2 \cdot \lambda, \quad \text{Energy flux} = E \cdot c = \frac{1}{2}\rho g \eta^2 \lambda \sqrt{gh}$$

100–500 km

Typical wavelength

10–60 min

Period range

<1 m

Open ocean height

Generation Mechanisms

Earthquake Generation (Okada Model)

Subduction zone earthquakes (M ≥ 7.5) produce vertical seafloor displacement over fault areas spanning hundreds of kilometers. The Okada (1985) dislocation model computes the seafloor deformation from fault parameters (length, width, dip, slip, depth):

$\Delta h = f(\text{strike}, \text{dip}, \text{rake}, \text{slip}, L, W, d)$

The 2011 Tohoku earthquake (M9.0) displaced the seafloor by up to 10 m over a 500 km × 200 km rupture area.

Submarine Landslide Tsunamis

Submarine landslides displace water volume proportional to the slide thickness and area. They can generate locally devastating tsunamis (1998 Papua New Guinea, >15 m run-up) but attenuate faster than seismic tsunamis due to shorter wavelengths. The Storegga Slide (~8,200 years ago) produced a major tsunami in the North Atlantic.

Volcanic Tsunamis

Volcanic flank collapse, caldera collapse, or pyroclastic flows entering the sea can generate tsunamis. The 1883 Krakatoa eruption produced waves reaching 30 m. The 2022 Hunga Tonga-Hunga Ha'apai eruption generated a unique atmospheric-pressure-driven tsunami detected globally.

Shoaling Amplification: Green's Law

As a tsunami enters shallow water, conservation of energy flux requires the wave height to increase as the depth decreases. Green's law provides the amplification factor:

$$\frac{H_2}{H_1} = \left(\frac{h_1}{h_2}\right)^{1/4} \cdot \left(\frac{b_1}{b_2}\right)^{1/2}$$

where b is the width of the wave ray tube (accounts for focusing/spreading)

Derivation: Green's Law

Energy flux must be conserved along a wave ray tube of width $b$:

$$E \cdot c_g \cdot b = \text{const along a ray tube}$$

The wave energy per unit area is $E = \tfrac{1}{2}\rho g H^2$, and in shallow water the group velocity equals the phase speed $c_g \approx c = \sqrt{gh}$. Equating energy flux at two locations:

$$\tfrac{1}{2}\rho g H_1^2 \sqrt{g h_1}\, b_1 = \tfrac{1}{2}\rho g H_2^2 \sqrt{g h_2}\, b_2$$

Cancel common factors ($\tfrac{1}{2}\rho g \sqrt{g}$) and solve for the height ratio:

$$H_1^2 \, h_1^{1/2} \, b_1 = H_2^2 \, h_2^{1/2} \, b_2$$

$$\frac{H_2^2}{H_1^2} = \frac{h_1^{1/2}}{h_2^{1/2}} \cdot \frac{b_1}{b_2}$$

Taking the square root:

$$\frac{H_2}{H_1} = \left(\frac{h_1}{h_2}\right)^{1/4} \cdot \left(\frac{b_1}{b_2}\right)^{1/2}$$

For a 1 m tsunami in 4000 m depth shoaling to 10 m depth (ignoring width changes):

$$H_2 = 1 \text{ m} \times \left(\frac{4000}{10}\right)^{1/4} = 1 \times 4.47 \approx 4.5 \text{ m}$$

Additional amplification from coastal geometry (funnel-shaped bays, headlands) and resonance can increase run-up further. The 2011 Tohoku tsunami reached >40 m run-up in some V-shaped valleys.

Run-up Estimation

The maximum run-up $R$ on a uniform slope is given by $R/H_0 = 2.83\sqrt{\cot\beta} \cdot (H_0/h_0)^{1/4}$ (Synolakis, 1987), where $\beta$ is the beach slope. Steeper beaches produce higher run-up for the same offshore wave height.

Dispersion Effects

Landslide-generated tsunamis have shorter wavelengths and are dispersive (waveform spreads out). For these events, Boussinesq equations capturing frequency dispersion are needed: $c^2 = gh(1 - k^2h^2/3)$. Dispersion reduces peak heights over transoceanic distances.

Tsunami Bore Formation

In very shallow water, the leading edge of a tsunami can steepen into a turbulent bore (hydraulic shock). The speed of the bore front is:

$$c_{\text{bore}} = \sqrt{g h_2 \cdot \frac{h_2 + h_1}{2 h_1}}$$

where $h_1$ is the undisturbed depth ahead of the bore and $h_2$ is the depth behind

Derivation: Bore Speed

Apply conservation of mass and momentum across the bore front. In a reference frame moving with the bore at speed $c$:

$$\text{Mass: } \quad h_1\, c = h_2\,(c - u_2)$$

$$\text{Momentum: } \quad h_1 c^2 + \tfrac{1}{2}g h_1^2 = h_2(c - u_2)^2 + \tfrac{1}{2}g h_2^2$$

From the mass equation: $c - u_2 = c\, h_1 / h_2$. Substitute into the momentum equation:

$$h_1 c^2 + \tfrac{1}{2}g h_1^2 = \frac{h_1^2 c^2}{h_2} + \tfrac{1}{2}g h_2^2$$

Rearrange and solve for $c^2$:

$$c^2 = g h_2 \cdot \frac{h_2 + h_1}{2 h_1}$$

Tsunami Warning Systems

Seismic Detection

Earthquake magnitude, location, and focal mechanism determined within minutes. Initial tsunami warning issued based on seismic parameters alone. Magnitude ≥ 7.5 at shallow depth triggers alert.

DART Buoys

Deep-ocean Assessment and Reporting of Tsunamis. Bottom pressure recorders detect tsunamis with 1 cm accuracy in 6000 m depth. ~60 DART stations in the Pacific, Atlantic, and Indian oceans.

Modern numerical tsunami models (MOST, GeoClaw, COMCOT) can produce inundation forecasts within 10–15 minutes of an earthquake, using pre-computed propagation databases and real-time DART data assimilation.

Coastal Resilience

Japan's multi-layered approach includes offshore breakwaters, coastal seawalls, elevated roads as secondary barriers, vertical evacuation buildings, and regular community drills. The 2011 Tohoku event demonstrated that engineered defenses alone are insufficient when events exceed design parameters — community preparedness and land-use planning are equally critical.

Major Historical Tsunamis

2004 Indian Ocean (M9.1)

230,000+ fatalities. 1,300 km rupture along Sunda Trench. Run-up >30 m in Sumatra. No warning system existed in the Indian Ocean. Led to establishment of Indian Ocean Tsunami Warning System.

2011 Tohoku, Japan (M9.0)

~20,000 fatalities. Max run-up ~40 m. Triggered Fukushima nuclear disaster. Tsunami exceeded seawall design heights (5–10 m). Arrived in 30 minutes locally, reached California in 10 hours.

Python: Tsunami Propagation, Green's Law & 1D Shallow Water Solver

Python

!/usr/bin/env python3

script.py103 lines

#!/usr/bin/env python3
"""tsunamis.py - Propagation time, Green's law, 1D shallow water solver"""
import numpy as np
import matplotlib.pyplot as plt

g = 9.81

# --- 1. Tsunami propagation time across Pacific ---
distances = np.array([0, 1000, 3000, 5000, 8000, 10000, 15000])  # km
labels = ['Source', 'Near coast', 'Hawaii', 'Mid-Pacific',
          'US West Coast', 'South America', 'Indian Ocean']
avg_depth = 4000  # m
c = np.sqrt(g * avg_depth)  # m/s
times_hr = (distances * 1000) / c / 3600

fig, axes = plt.subplots(2, 2, figsize=(14, 10))

ax1 = axes[0, 0]
ax1.plot(distances, times_hr, 'r-o', lw=2, ms=6)
for d, t, lab in zip(distances[1:], times_hr[1:], labels[1:]):
    ax1.annotate(lab, (d, t), textcoords='offset points',
                 xytext=(5,5), fontsize=7, color='white')
ax1.set_xlabel('Distance (km)'); ax1.set_ylabel('Travel time (hours)')
ax1.set_title('Tsunami Propagation Time (h=4000 m)')
ax1.grid(True, alpha=0.3)

# --- 2. Green's law shoaling ---
ax2 = axes[0, 1]
h_deep = 4000
H_deep = np.array([0.1, 0.5, 1.0, 2.0])  # initial height (m)
h_shallow = np.linspace(4000, 1, 500)

for H0 in H_deep:
    H_shoal = H0 * (h_deep / h_shallow)**0.25
    ax2.semilogy(h_shallow, H_shoal, lw=2, label=f'H0 = {H0} m')

ax2.axhline(10, color='r', ls='--', alpha=0.5, label='10 m threshold')
ax2.set_xlabel('Water depth (m)'); ax2.set_ylabel('Wave height (m)')
ax2.set_title("Green's Law Shoaling Amplification")
ax2.legend(fontsize=8); ax2.grid(True, alpha=0.3)
ax2.invert_xaxis()

# --- 3. Speed vs depth ---
ax3 = axes[1, 0]
depths = np.linspace(1, 6000, 500)
speeds = np.sqrt(g * depths) * 3.6  # km/hr
ax3.plot(depths, speeds, 'c-', lw=2)
ax3.set_xlabel('Depth (m)'); ax3.set_ylabel('Speed (km/hr)')
ax3.set_title('Tsunami Speed vs Depth')
ax3.grid(True, alpha=0.3)
ax3.axhline(900, color='orange', ls='--', label='Commercial jet')
ax3.legend()

# --- 4. 1D shallow water equation solver ---
ax4 = axes[1, 1]
nx = 2000; Lx = 500e3  # 500 km domain
dx = Lx / nx
x = np.linspace(0, Lx, nx)

# Bathymetry: deep ocean + continental shelf
h = np.where(x < 400e3, 4000.0, 4000 - (x - 400e3)/100e3 * 3990)
h = np.maximum(h, 5.0)

# Initial tsunami (Gaussian)
eta = 2.0 * np.exp(-((x - 50e3) / 20e3)**2)
u_vel = eta * np.sqrt(g / h)  # rightward propagation

dt = 0.4 * dx / np.sqrt(g * 4000)
nsteps = int(1800 / dt)  # 30 min simulation

for step in range(nsteps):
    # Flux form of shallow water equations
    flux_mass = (h + eta) * u_vel
    flux_mom = u_vel**2 * (h + eta) + 0.5 * g * (h + eta)**2

# Finite difference (upwind-like Lax-Friedrichs)
    eta_new = eta.copy()
    u_new = u_vel.copy()
    for i in range(1, nx-1):
        eta_new[i] = 0.5*(eta[i+1]+eta[i-1]) - dt/(2*dx) * \
                     (flux_mass[i+1] - flux_mass[i-1])
        H_total = h[i] + eta_new[i]
        if H_total > 0.01:
            mom = (h[i]+eta[i])*u_vel[i]
            mom = 0.5*((h[i+1]+eta[i+1])*u_vel[i+1] + \
                       (h[i-1]+eta[i-1])*u_vel[i-1]) - \
                  dt/(2*dx)*(flux_mom[i+1]-flux_mom[i-1]) + \
                  dt * g * H_total * (h[i+1]-h[i-1])/(2*dx)
            u_new[i] = mom / H_total
    eta_new[0] = eta_new[1]; eta_new[-1] = eta_new[-2]
    u_new[0] = u_new[1]; u_new[-1] = 0.0
    eta = eta_new; u_vel = u_new

ax4.plot(x/1e3, eta, 'b-', lw=1.5, label='Surface')
ax4.fill_between(x/1e3, -h/100, alpha=0.3, color='brown',
                 label='Seafloor (scaled)')
ax4.set_xlabel('Distance (km)'); ax4.set_ylabel('Elevation (m)')
ax4.set_title('1D Tsunami Propagation (t=30 min)')
ax4.legend(); ax4.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('tsunamis.png', dpi=150)
plt.show()

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Fortran: 2D Shallow Water Equation Solver for Tsunami Propagation

Fortran

2D shallow water equations for tsunami propagation

program.f9092 lines

program tsunami_2d
  ! 2D shallow water equations for tsunami propagation
  ! Finite difference method on a regular grid
  implicit none
  integer, parameter :: dp = selected_real_kind(15)
  real(dp), parameter :: g = 9.81_dp, pi = 3.141592653589793_dp
  integer, parameter :: nx = 200, ny = 200
  real(dp) :: Lx, Ly, dx, dy, dt
  real(dp) :: eta(nx,ny), u(nx,ny), v(nx,ny)
  real(dp) :: eta_new(nx,ny), u_new(nx,ny), v_new(nx,ny)
  real(dp) :: h(nx,ny), H_total
  real(dp) :: x, y, r, max_eta
  integer :: i, j, step, nsteps

Lx = 1000.0e3_dp; Ly = 1000.0e3_dp   ! 1000 km domain
  dx = Lx / real(nx-1, dp)
  dy = Ly / real(ny-1, dp)

! Setup bathymetry: flat ocean + island shelf
  do j = 1, ny
    do i = 1, nx
      h(i,j) = 4000.0_dp   ! uniform 4 km depth
      ! Add continental shelf on right boundary
      x = real(i-1, dp) * dx
      if (x > 800.0e3_dp) then
        h(i,j) = 4000.0_dp - (x - 800.0e3_dp)/200.0e3_dp * 3950.0_dp
        h(i,j) = max(h(i,j), 5.0_dp)
      end if
    end do
  end do

! Initial tsunami: Gaussian displacement at source
  eta = 0.0_dp; u = 0.0_dp; v = 0.0_dp
  do j = 1, ny
    do i = 1, nx
      x = real(i-1, dp) * dx - 200.0e3_dp
      y = real(j-1, dp) * dy - 500.0e3_dp
      r = sqrt(x**2 + y**2)
      eta(i,j) = 3.0_dp * exp(-(r / 50.0e3_dp)**2)
    end do
  end do

! CFL condition
  dt = 0.3_dp * min(dx, dy) / sqrt(g * 4000.0_dp)
  nsteps = nint(3600.0_dp / dt)  ! 1 hour simulation

print *, '2D Tsunami Propagation Solver'
  print '(A,I4,A,I4)', ' Grid: ', nx, ' x ', ny
  print '(A,F6.1,A)', ' dt = ', dt, ' s'
  print '(A,I6)', ' Steps for 1 hour: ', nsteps

do step = 1, nsteps
    ! Update velocities (momentum equations)
    do j = 2, ny-1
      do i = 2, nx-1
        H_total = h(i,j) + eta(i,j)
        if (H_total > 1.0_dp) then
          u_new(i,j) = u(i,j) - dt * g * &
                       (eta(i+1,j) - eta(i-1,j)) / (2.0_dp * dx)
          v_new(i,j) = v(i,j) - dt * g * &
                       (eta(i,j+1) - eta(i,j-1)) / (2.0_dp * dy)
        end if
      end do
    end do

! Update surface elevation (continuity)
    do j = 2, ny-1
      do i = 2, nx-1
        H_total = h(i,j) + eta(i,j)
        eta_new(i,j) = eta(i,j) - dt * ( &
          H_total * (u_new(i+1,j) - u_new(i-1,j)) / (2.0_dp*dx) + &
          H_total * (v_new(i,j+1) - v_new(i,j-1)) / (2.0_dp*dy) )
      end do
    end do

! Boundary conditions (open/absorbing)
    eta_new(1,:) = eta_new(2,:); eta_new(nx,:) = 0.0_dp
    eta_new(:,1) = eta_new(:,2); eta_new(:,ny) = eta_new(:,ny-1)
    u_new(1,:) = u_new(2,:); u_new(nx,:) = 0.0_dp
    v_new(:,1) = 0.0_dp; v_new(:,ny) = 0.0_dp

eta = eta_new; u = u_new; v = v_new

if (mod(step, nsteps/5) == 0) then
      max_eta = maxval(abs(eta))
      print '(A,I6,A,F8.3,A)', ' Step ', step, &
            '  max|eta| = ', max_eta, ' m'
    end if
  end do

print *, 'Simulation complete.'
end program tsunami_2d

Click Run to execute the Fortran code

Code will be compiled with gfortran and executed on the server

Numerical Tsunami Modeling

Modern tsunami forecasting relies on numerical solution of the nonlinear shallow water equations (NLSW):

$$\frac{\partial \eta}{\partial t} + \frac{\partial}{\partial x}[(h+\eta)u] + \frac{\partial}{\partial y}[(h+\eta)v] = 0$$

$$\frac{\partial u}{\partial t} + u\frac{\partial u}{\partial x} + v\frac{\partial u}{\partial y} + g\frac{\partial \eta}{\partial x} = -\frac{\tau_x^b}{\rho(h+\eta)}$$

Bottom friction $\tau_b$ becomes important in shallow water and during inundation

MOST (NOAA)

Method of Splitting Tsunami. Uses pre-computed propagation database (unit sources) combined in real-time with earthquake parameters. Provides forecast within 10–15 minutes of earthquake detection.

GeoClaw

Adaptive mesh refinement solver. Automatically increases resolution near coastlines and wave fronts. Handles wetting/drying for inundation. Open-source (Clawpack library).

COMCOT

Cornell Multi-grid Coupled Tsunami Model. Nested grids from ocean basin (~1 km) to coast (~10 m). Linear equations offshore, nonlinear nearshore. Validated against 2004 and 2011 events.

Boussinesq Models

Include frequency dispersion (important for landslide tsunamis with shorter wavelengths). FUNWAVE, COULWAVE. Higher computational cost but more accurate for dispersive waves.

Probabilistic Tsunami Hazard Assessment (PTHA)

PTHA provides the scientific basis for coastal planning and emergency preparedness by estimating the probability of tsunami inundation at a given location over a specified time period. The framework parallels probabilistic seismic hazard analysis:

$$P(\eta > \eta_0) = 1 - \exp\left(-\nu T \cdot P(\eta > \eta_0 | \text{event})\right)$$

where ν is the annual rate of tsunamigenic events, T is the exposure period, and the conditional probability integrates over all source scenarios

PTHA involves three key steps:

Source Characterization

Develop earthquake source models: magnitude-frequency distribution (Gutenberg-Richter), fault geometry, slip distribution. Include epistemic uncertainties with logic trees. Typically thousands of scenarios per source zone.

Propagation Modeling

Run numerical simulations for each source scenario. Use pre-computed Green's functions or unit-source databases to rapidly estimate wave heights at coastal points. Account for bathymetric focusing and refraction.

Hazard Aggregation

Combine all sources to produce exceedance probability curves and hazard maps. The 500-year return period wave height is commonly used for design standards. Results inform building codes, evacuation zones, and insurance.

Key Concepts Summary

Physics

Tsunamis are shallow water waves: $c = \sqrt{gh}$. Non-dispersive propagation preserves waveform. Speed ~700 km/hr in deep ocean, ~30 km/hr at coast.

Green's Law

Shoaling amplification: $H_2/H_1 = (h_1/h_2)^{1/4}$. A 1 m deep-ocean wave can amplify to 5+ m at the coast. Focusing in bays further increases height.

Warning Systems

Seismic detection (<5 min), DART buoys (real-time confirmation), numerical models (forecast in 10–15 min). Pacific TWS operational since 1949; Indian Ocean since 2006.

Hazard Assessment

Probabilistic tsunami hazard analysis (PTHA) maps inundation risk for coastal planning. Run-up depends on local bathymetry, coastal morphology, and wave direction.

← 7.4 Tidal Patterns 7.6 Coastal Processes →