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

script.py88 lines

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

Fortran

Numerical dispersion relation solver (Newton-Raphson)

program.f9073 lines

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.