7.1 Wave Dynamics
Ocean surface waves are the primary manifestation of energy transfer from the atmosphere to the ocean. Linear wave theory (Airy wave theory) provides the foundational mathematical framework for describing wave propagation, energy transport, and particle kinematics across all water depths.
Linear (Airy) Wave Theory
The linearized boundary value problem for surface gravity waves yields a sinusoidal free surface displacement and a velocity potential that satisfies Laplace's equation within the fluid domain. For small-amplitude waves on an inviscid, irrotational fluid, the surface elevation is:
$$\eta(x,t) = \frac{H}{2}\cos(kx - \omega t)$$
where $H$ is the wave height, $k = 2\pi/\lambda$ is the wavenumber, and $\omega = 2\pi/T$ is the angular frequency
The velocity potential satisfying Laplace's equation with the linearized free surface and bottom boundary conditions is:
$$\phi(x,z,t) = \frac{gH}{2\omega}\frac{\cosh k(z+h)}{\cosh kh}\sin(kx - \omega t)$$
Assumptions
Inviscid fluid, irrotational flow, incompressible, constant depth, small amplitude (H/ฮป ยซ 1), no surface tension
Particle Orbits
Circular in deep water, elliptical in intermediate/shallow water. Orbital radius decays as $e^{kz}$ with depth.
Dispersion Relation
The dispersion relation links wave frequency to wavenumber and water depth. It arises from applying the dynamic and kinematic boundary conditions at the free surface:
$$\omega^2 = gk\tanh(kh)$$
This transcendental equation must generally be solved numerically for k given ฯ and h
Derivation: Dispersion Relation
Step 1. Begin with the governing equation for an inviscid, irrotational, incompressible fluid. The velocity potential $\phi$ satisfies Laplace's equation throughout the fluid domain:
$$\nabla^2 \phi = \frac{\partial^2 \phi}{\partial x^2} + \frac{\partial^2 \phi}{\partial z^2} = 0 \quad \text{for } -h \leq z \leq 0$$
Step 2. Apply the boundary conditions. At the impermeable bottom, the vertical velocity vanishes:
$$\frac{\partial \phi}{\partial z} = 0 \quad \text{at } z = -h$$
At the free surface ($z = 0$ after linearization), we have the kinematic and dynamic boundary conditions:
$$\text{Kinematic: } \frac{\partial \eta}{\partial t} = \frac{\partial \phi}{\partial z} \quad \text{at } z = 0$$
$$\text{Dynamic: } \frac{\partial \phi}{\partial t} + g\eta = 0 \quad \text{at } z = 0$$
Step 3. Assume a separable, traveling-wave solution of the form:
$$\phi(x,z,t) = Z(z)\sin(kx - \omega t)$$
Substituting into Laplace's equation gives $Z'' - k^2 Z = 0$, whose general solution is:
$$Z(z) = A\cosh k(z+h) + B\sinh k(z+h)$$
Step 4. Apply the bottom BC: $Z'(-h) = 0$ requires $B = 0$. Therefore:
$$Z(z) = A\cosh k(z+h)$$
Step 5. Combine the two surface BCs. Eliminate $\eta$ by differentiating the dynamic condition with respect to $t$ and substituting the kinematic condition:
$$\frac{\partial^2 \phi}{\partial t^2} + g\frac{\partial \phi}{\partial z} = 0 \quad \text{at } z = 0$$
Step 6. Substitute $\phi = A\cosh k(z+h)\sin(kx - \omega t)$ into the combined surface condition at $z = 0$:
$$-\omega^2 A\cosh(kh) + gkA\sinh(kh) = 0$$
For a non-trivial solution ($A \neq 0$), we require:
$$\boxed{\omega^2 = gk\tanh(kh)}$$
The phase speed (celerity) and group velocity follow directly:
$$c = \frac{\omega}{k} = \sqrt{\frac{g}{k}\tanh(kh)}, \quad c_g = \frac{\partial \omega}{\partial k} = \frac{c}{2}\left(1 + \frac{2kh}{\sinh 2kh}\right)$$
Derivation: Group Velocity
Step 1. The group velocity is defined as the rate of change of frequency with respect to wavenumber:
$$c_g = \frac{\partial \omega}{\partial k}$$
Step 2. Start from the dispersion relation $\omega^2 = gk\tanh(kh)$ and differentiate both sides implicitly with respect to $k$:
$$2\omega\frac{\partial \omega}{\partial k} = g\tanh(kh) + gk \cdot h\,\mathrm{sech}^2(kh)$$
Step 3. Solve for $c_g = \partial\omega/\partial k$. Note that $c = \omega/k$, so $\omega = ck$ and $\tanh(kh) = \omega^2/(gk)$:
$$c_g = \frac{g}{2\omega}\left[\tanh(kh) + kh\,\mathrm{sech}^2(kh)\right]$$
Step 4. Factor out $\tanh(kh)$ from the bracket. Using $c = (g/\omega)\tanh(kh)$ (from the dispersion relation divided by $k$):
$$c_g = \frac{c}{2}\left[1 + \frac{kh\,\mathrm{sech}^2(kh)}{\tanh(kh)}\right]$$
Step 5. Simplify the second term using $\mathrm{sech}^2(kh)/\tanh(kh) = 1/(\sinh(kh)\cosh(kh)) = 2/\sinh(2kh)$:
$$\boxed{c_g = \frac{c}{2}\left(1 + \frac{2kh}{\sinh 2kh}\right)}$$
In deep water ($kh \to \infty$), $2kh/\sinh(2kh) \to 0$, giving $c_g = c/2$. In shallow water ($kh \to 0$), $2kh/\sinh(2kh) \to 1$, giving $c_g = c$.
Deep Water (kh ยป 1)
When depth exceeds half the wavelength, tanh(kh) โ 1:
$c = \frac{g}{\omega} = \frac{gT}{2\pi}, \quad c_g = \frac{c}{2}$
Waves are dispersive: longer waves travel faster. Group velocity is half the phase speed.
Shallow Water (kh ยซ 1)
When depth is much less than wavelength, tanh(kh) โ kh:
$c = \sqrt{gh}, \quad c_g = c$
Non-dispersive: all wavelengths travel at the same speed, depending only on depth.
Wave Energy and Energy Flux
The total energy per unit surface area of a linear wave is equally partitioned between kinetic and potential energy:
$$E = \frac{1}{2}\rho g a^2 = \frac{1}{8}\rho g H^2$$
where $a = H/2$ is the wave amplitude
Derivation: Wave Energy and Equipartition
Step 1. The potential energy per unit horizontal area is found by computing the work done against gravity to deform the free surface from its rest position. Averaging over one wavelength:
$$\langle PE \rangle = \frac{1}{\lambda}\int_0^{\lambda} \frac{1}{2}\rho g \eta^2 \,dx$$
Substituting $\eta = a\cos(kx - \omega t)$ and using $\langle\cos^2\rangle = 1/2$:
$$\langle PE \rangle = \frac{1}{2}\rho g a^2 \cdot \frac{1}{2} = \frac{1}{4}\rho g a^2$$
Step 2. The kinetic energy per unit horizontal area is found by integrating the squared velocity over the water column. The horizontal and vertical velocity components from linear theory are:
$$u = a\omega\frac{\cosh k(z+h)}{\sinh(kh)}\cos(kx-\omega t), \quad w = a\omega\frac{\sinh k(z+h)}{\sinh(kh)}\sin(kx-\omega t)$$
Step 3. Compute the wavelength-averaged kinetic energy:
$$\langle KE \rangle = \frac{1}{\lambda}\int_0^{\lambda}\int_{-h}^{0} \frac{1}{2}\rho(u^2 + w^2)\,dz\,dx$$
After evaluating the depth and wavelength integrals (using $\langle\cos^2\rangle = \langle\sin^2\rangle = 1/2$ and the identity $\int_{-h}^{0}\cosh^2 k(z+h)\,dz = \frac{1}{2}\left(h + \frac{\sinh 2kh}{2k}\right)$), the result simplifies to:
$$\langle KE \rangle = \frac{1}{4}\rho g a^2$$
Step 4. This establishes the equipartition of wave energy: kinetic and potential energies are equal. The total energy per unit surface area is:
$$\boxed{E = \langle KE \rangle + \langle PE \rangle = \frac{1}{4}\rho g a^2 + \frac{1}{4}\rho g a^2 = \frac{1}{2}\rho g a^2 = \frac{1}{8}\rho g H^2}$$
where $a = H/2$. This equipartition result is a fundamental property of linear surface gravity waves and holds for all water depths.
The energy flux (power per unit crest width) is the product of energy density and group velocity:
$$P = E \cdot c_g = \frac{1}{8}\rho g H^2 \cdot c_g$$
For a typical 2 m swell with 10 s period in deep water, $c_g \approx 7.8$ m/s, yielding P โ 39 kW/m. This is why wave energy is a promising renewable resource. The global wave power resource is estimated at approximately 2 TW.
Conservation of Wave Energy Flux
As waves enter shallow water and c_g changes, conservation of energy flux $E \cdot c_g = \text{const}$ along wave rays causes wave height to increase (shoaling).
Stokes Drift
Although linear wave theory predicts closed circular particle orbits, second-order analysis reveals a net forward mass transport known as Stokes drift. Particles do not return to their starting positions after one wave period:
$$u_S(z) = \frac{\omega k H^2}{8} \frac{\cosh 2k(z+h)}{\sinh^2(kh)}$$
In deep water this simplifies to $u_S = \frac{\omega k H^2}{8} e^{2kz}$
Derivation: Stokes Drift
Step 1. Consider the difference between the Lagrangian (following a fluid particle) and Eulerian (fixed point) descriptions. Let a particle's position be $\vec{x}(t) = \vec{X} + \vec{\xi}(t)$, where $\vec{X}$ is the mean position and $\vec{\xi}$ is the displacement. The Lagrangian velocity is:
$$\vec{u}_L(t) = \vec{u}(\vec{X} + \vec{\xi}, t)$$
Step 2. Taylor-expand the Eulerian velocity about the mean position to second order:
$$\vec{u}_L \approx \vec{u}(\vec{X},t) + (\vec{\xi} \cdot \nabla)\vec{u}(\vec{X},t) + \cdots$$
Step 3. At first order in wave amplitude, the particle displacement is $\vec{\xi}^{(1)} = \int \vec{u}^{(1)} dt$, and the Eulerian velocity $\vec{u}^{(1)}$ has zero time-average. The Stokes drift is the time-averaged second-order correction:
$$\vec{u}_S = \overline{(\vec{\xi}^{(1)} \cdot \nabla)\vec{u}^{(1)}}$$
Step 4. For a monochromatic wave with $u^{(1)} = a\omega\frac{\cosh k(z+h)}{\sinh kh}\cos(kx - \omega t)$, the first-order horizontal displacement is:
$$\xi^{(1)} = -\frac{a\cosh k(z+h)}{\sinh kh}\sin(kx - \omega t)$$
Step 5. Compute the correlation $u_S = \overline{\xi^{(1)}\frac{\partial u^{(1)}}{\partial x} + \zeta^{(1)}\frac{\partial u^{(1)}}{\partial z}}$, where $\zeta^{(1)}$ is the vertical displacement. The dominant contribution comes from:
$$u_S = \overline{\xi^{(1)}\frac{\partial u^{(1)}}{\partial x}} + \overline{\zeta^{(1)}\frac{\partial u^{(1)}}{\partial z}}$$
Using $\overline{\sin^2\theta} = \overline{\cos^2\theta} = 1/2$ and combining the $\cosh^2$ and $\sinh^2$ terms via $\cosh^2 + \sinh^2 = \cosh(2\cdot)$:
$$\boxed{u_S(z) = \frac{\omega k a^2 \cosh 2k(z+h)}{2\sinh^2(kh)} = \frac{\omega k H^2}{8}\frac{\cosh 2k(z+h)}{\sinh^2(kh)}}$$
The drift is always in the direction of wave propagation and is largest at the surface. It arises because forward orbital velocity at the crest (where the particle is higher) slightly exceeds backward velocity at the trough (where the particle is lower).
Stokes drift is critical for understanding surface material transport (oil spills, microplastics, larvae), and the Stokes-Coriolis force drives Langmuir circulations. The vertically-integrated Stokes transport equals $M_S = E/c$.
Langmuir Circulations
Counter-rotating helical vortices aligned with the wind direction, driven by the interaction of Stokes drift and wind-driven shear (Craik-Leibovich mechanism). Visible as windrows of foam on the surface. Cell spacing is 10โ300 m. The turbulent Langmuir number $La_t = \sqrt{u_*/u_{S0}}$ characterizes their intensity.
Wave-Mean Flow Interaction
The vortex force $\vec{u}_S \times \vec{\omega}$ (where $\vec{\omega}$ is the mean flow vorticity) couples wave effects to the Eulerian mean circulation. This is included in wave-current coupled models through the Craik-Leibovich equations and is important for modeling upper ocean mixing and mixed layer deepening.
Wave Spectra
Real ocean waves are irregular and best described statistically using spectral analysis. The wave energy spectrum S(f) describes how wave energy is distributed across frequencies:
$$\text{Pierson-Moskowitz: } S(f) = \frac{\alpha g^2}{(2\pi)^4 f^5}\exp\left[-\frac{5}{4}\left(\frac{f_p}{f}\right)^4\right]$$
$$\text{JONSWAP: } S_J(f) = S_{PM}(f) \cdot \gamma^{\exp\left[-\frac{(f-f_p)^2}{2\sigma^2 f_p^2}\right]}$$
$\gamma \approx 3.3$ is the peak enhancement factor; $\sigma = 0.07$ ($f \leq f_p$) or $0.09$ ($f > f_p$)
The significant wave height is related to the zeroth spectral moment:
$$H_s = H_{m0} = 4\sqrt{m_0}, \quad m_n = \int_0^\infty f^n S(f)\,df$$
Wave Generation by Wind
Two complementary mechanisms explain how wind generates waves:
Phillips (1957) โ Resonance Mechanism
Turbulent pressure fluctuations in the atmospheric boundary layer excite waves at frequencies that match the convection speed of the pressure pattern. This linear resonance drives initial wave growth from a flat surface. Growth is linear in time: $E \propto t$.
Miles (1957) โ Shear Flow Instability
Once waves reach finite amplitude, the wave-induced pressure perturbation at the critical layer (where wind speed equals phase speed) extracts energy from the mean wind profile. Growth becomes exponential: $E \propto e^{\beta t}$, where $\beta$ depends on the wind profile curvature at the critical height.
Fully Developed Sea
When energy input from wind balances dissipation (whitecapping) and nonlinear transfer, the spectrum reaches equilibrium. The PM spectrum describes this state. For 20 m/s wind: H_s โ 5 m, T_p โ 12.5 s.
Python: Dispersion, Spectra & Stokes Drift
Python: Dispersion, Spectra & Stokes Drift
Python!/usr/bin/env python3
Click Run to execute the Python code
Code will be executed with Python 3 on the server
Fortran: Dispersion Solver & Wave Spectrum
Fortran: Dispersion Solver & Wave Spectrum
FortranNumerical dispersion relation solver (Newton-Raphson)
Click Run to execute the Fortran code
Code will be compiled with gfortran and executed on the server
Nonlinear Wave Effects
Linear (Airy) wave theory assumes small amplitude (H/ฮป ยซ 1), but real ocean waves exhibit nonlinear behavior. Stokes wave theory includes higher-order harmonics to describe finite-amplitude waves:
$$\eta(x,t) = a\cos\theta + \frac{ka^2}{2}\cos 2\theta + \frac{3k^2 a^3}{8}\cos 3\theta + \mathcal{O}(a^4)$$
where $\theta = kx - \omega t$ and higher harmonics sharpen crests and flatten troughs
Key nonlinear phenomena include:
Wave-Wave Interactions
Resonant triad and quartet interactions transfer energy across the spectrum. Phillips (1960) showed that four gravity waves satisfying $\vec{k}_1 + \vec{k}_2 = \vec{k}_3 + \vec{k}_4$ and $\omega_1 + \omega_2 = \omega_3 + \omega_4$ can exchange energy. This nonlinear transfer shapes the spectral peak and high-frequency tail.
Wave Breaking and Whitecapping
Deep-water breaking occurs when wave steepness approaches the Stokes limit: $H/\lambda \leq 1/7 \approx 0.142$. Whitecapping dissipates energy from the spectral peak and is the primary sink balancing wind input. Phillips (1985) spectral tail: $S(f) \propto f^{-5}$ in the equilibrium range.
Cnoidal and Solitary Waves
In very shallow water (Ursell number $U_r = HL^2/h^3 > 26$), the linear sinusoidal profile is replaced by cnoidal wave profiles described by Jacobian elliptic functions. In the limit of very high $U_r$, these become solitary waves (single humps) with speed $c = \sqrt{g(h+H)}$.
Infragravity Wave Generation
Nonlinear interaction between wave groups generates bound long waves (infragravity waves) at the group frequency. These are released upon wave breaking in the surf zone and can resonantly amplify in harbors (periods 30 s to 5 min).
Operational Wave Forecasting
Modern wave forecasting models solve the spectral wave action equation, which describes the evolution of the wave spectrum in space and time:
$$\frac{\partial N}{\partial t} + \nabla \cdot (\vec{c}_g N) = \frac{S_{\text{in}} + S_{\text{nl}} + S_{\text{ds}}}{\sigma}$$
$N = S/\sigma$ is the wave action density; $S_{\text{in}}$ = wind input, $S_{\text{nl}}$ = nonlinear transfer, $S_{\text{ds}}$ = dissipation
WAM (Wave Model)
Third-generation spectral wave model. ECMWF operational since 1992. Solves the full wave action equation without any a priori assumptions about the spectral shape. Includes state-of-the-art source term parameterizations.
WAVEWATCH III
NOAA operational model. Includes unstructured grids, ice interaction, and coupling with ocean circulation models. Global 0.5ยฐ and regional 1/12ยฐ nests.
SWAN
Simulating WAves Nearshore. Designed for coastal applications with depth-induced breaking, triad interactions, and diffraction. Typically used for nearshore domains.
Key Concepts Summary
Dispersion Regimes
Deep water: dispersive, c ∝ T. Shallow water: non-dispersive, c = โ(gh). Intermediate: full tanh solution required.
Energy Transport
Energy propagates at group velocity $c_g$, not phase speed $c$. In deep water $c_g = c/2$; in shallow water $c_g = c$.
Wave Spectra
PM spectrum for fully developed seas; JONSWAP for fetch-limited growth. H_s = 4โm' from spectral moments.
Stokes Drift
Net forward mass transport at O(ka)ยฒ. Exponential decay with depth. Critical for pollutant and larval transport.