Wave Theory

Wave Theory SME Skill

Comprehensive ocean wave theory expertise including wave mechanics, spectral analysis, wave statistics, and irregular sea modeling for offshore engineering applications.

When to Use This Skill

Use wave theory knowledge when:

• Wave spectra - JONSWAP, Pierson-Moskowitz, scatter diagrams
• Wave statistics - Significant wave height, spectral parameters
• Irregular seas - Generate time series from spectra
• Wave kinematics - Particle velocities and accelerations
• Wave transformation - Shoaling, refraction, diffraction
• Extreme values - Design wave estimation

Core Knowledge Areas

1. Regular Wave Theory

Linear (Airy) Wave Theory:

import numpy as np

def airy_wave_properties(
    H: float,
    T: float,
    d: float,
    g: float = 9.81
) -> dict:
    """
    Calculate Airy wave properties.

    Valid for: H/L < 0.14, d/L > 0.5 (deep water) or d/L < 0.05 (shallow)

    Args:
        H: Wave height (m)
        T: Wave period (s)
        d: Water depth (m)
        g: Gravity (m/s²)

    Returns:
        Wave properties dictionary
    """
    # Wave frequency
    omega = 2 * np.pi / T

    # Dispersion relation: ω² = gk·tanh(kd)
    # Solve iteratively for wave number k
    from scipy.optimize import fsolve

    def dispersion(k):
        return omega**2 - g * k * np.tanh(k * d)

    k0 = omega**2 / g  # Deep water approximation
    k = fsolve(dispersion, k0)[0]

    # Wave length
    L = 2 * np.pi / k

    # Wave celerity (phase speed)
    C = omega / k

    # Group velocity
    n = 0.5 * (1 + 2*k*d / np.sinh(2*k*d))  # Shoaling coefficient
    Cg = n * C

    # Deep water classification
    if d / L > 0.5:
        regime = "Deep water"
    elif d / L < 0.05:
        regime = "Shallow water"
    else:
        regime = "Intermediate"

    return {
        'height_m': H,
        'period_s': T,
        'depth_m': d,
        'wavelength_m': L,
        'wave_number': k,
        'frequency_rad_s': omega,
        'frequency_hz': omega / (2*np.pi),
        'celerity_m_s': C,
        'group_velocity_m_s': Cg,
        'regime': regime,
        'd_over_L': d / L,
        'H_over_L': H / L,
        'steepness': H / L
    }

# Example
wave = airy_wave_properties(H=8, T=12, d=1500)

print(f"Wave Properties (H={wave['height_m']}m, T={wave['period_s']}s):")
print(f"  Wavelength: {wave['wavelength_m']:.1f} m")
print(f"  Regime: {wave['regime']} (d/L = {wave['d_over_L']:.3f})")
print(f"  Celerity: {wave['celerity_m_s']:.2f} m/s")
print(f"  Steepness: {wave['steepness']:.4f}")

Wave Kinematics:

def wave_particle_kinematics(
    z: float,
    H: float,
    T: float,
    d: float,
    t: float = 0,
    x: float = 0,
    g: float = 9.81
) -> dict:
    """
    Calculate wave particle velocities and accelerations.

    Args:
        z: Vertical position (0 at SWL, negative below)
        H: Wave height (m)
        T: Wave period (s)
        d: Water depth (m)
        t: Time (s)
        x: Horizontal position (m)
        g: Gravity (m/s²)

    Returns:
        Particle kinematics
    """
    # Wave properties
    wave = airy_wave_properties(H, T, d, g)
    k = wave['wave_number']
    omega = wave['frequency_rad_s']

    # Amplitude
    a = H / 2

    # Hyperbolic functions
    cosh_kz_d = np.cosh(k * (z + d))
    sinh_kz_d = np.sinh(k * (z + d))
    cosh_kd = np.cosh(k * d)
    sinh_kd = np.sinh(k * d)

    # Wave phase
    phase = k * x - omega * t

    # Horizontal velocity
    u = (omega * a * cosh_kz_d / sinh_kd) * np.cos(phase)

    # Vertical velocity
    w = (omega * a * sinh_kz_d / sinh_kd) * np.sin(phase)

    # Horizontal acceleration
    ax = -(omega**2 * a * cosh_kz_d / sinh_kd) * np.sin(phase)

    # Vertical acceleration
    az = (omega**2 * a * sinh_kz_d / sinh_kd) * np.cos(phase)

    # Dynamic pressure
    p_dynamic = g * a * (cosh_kz_d / cosh_kd) * np.cos(phase)

    return {
        'horizontal_velocity': u,
        'vertical_velocity': w,
        'horizontal_acceleration': ax,
        'vertical_acceleration': az,
        'dynamic_pressure': p_dynamic,
        'total_velocity': np.sqrt(u**2 + w**2),
        'total_acceleration': np.sqrt(ax**2 + az**2)
    }

# Example: Surface velocity (z=0)
kinematics = wave_particle_kinematics(z=0, H=8, T=12, d=1500, t=0, x=0)

print(f"Surface Particle Kinematics:")
print(f"  Horizontal velocity: {kinematics['horizontal_velocity']:.2f} m/s")
print(f"  Vertical velocity: {kinematics['vertical_velocity']:.2f} m/s")
print(f"  Total velocity: {kinematics['total_velocity']:.2f} m/s")

2. Wave Spectra

JONSWAP Spectrum:

def jonswap_spectrum(
    frequencies: np.ndarray,
    Hs: float,
    Tp: float,
    gamma: float = 3.3,
    alpha: float = None
) -> np.ndarray:
    """
    Calculate JONSWAP wave spectrum.

    S(f) = α g² (2π)^-4 f^-5 exp[-5/4(f/fp)^-4] γ^exp[-(f-fp)²/(2σ²fp²)]

    Args:
        frequencies: Frequency array (Hz)
        Hs: Significant wave height (m)
        Tp: Peak period (s)
        gamma: Peak enhancement factor (3.3 for North Sea)
        alpha: Phillips constant (calculated if None)

    Returns:
        Spectral density S(f) (m²/Hz)
    """
    g = 9.81
    fp = 1 / Tp  # Peak frequency (Hz)

    # Calculate alpha if not provided
    if alpha is None:
        # Relationship: Hs = 4*sqrt(m0)
        # For JONSWAP: alpha ≈ 5.061 * Hs² / Tp⁴ * (1 - 0.287*ln(γ))
        alpha = 5.061 * Hs**2 / Tp**4 * (1 - 0.287 * np.log(gamma))

    # Sigma parameter
    sigma = np.where(frequencies <= fp, 0.07, 0.09)

    # Pierson-Moskowitz spectrum
    S_PM = alpha * g**2 * (2*np.pi)**(-4) * frequencies**(-5) * \
           np.exp(-1.25 * (frequencies / fp)**(-4))

    # Peak enhancement
    r = np.exp(-(frequencies - fp)**2 / (2 * sigma**2 * fp**2))
    gamma_factor = gamma ** r

    # JONSWAP spectrum
    S = S_PM * gamma_factor

    return S

# Example: Generate JONSWAP spectrum
freq = np.linspace(0.01, 0.5, 500)
S = jonswap_spectrum(freq, Hs=8.5, Tp=12.0, gamma=3.3)

# Verify Hs
m0 = np.trapz(S, freq)
Hs_calc = 4 * np.sqrt(m0)

print(f"JONSWAP Spectrum:")
print(f"  Input Hs: 8.5 m")
print(f"  Calculated Hs: {Hs_calc:.2f} m")
print(f"  Peak frequency: {1/12:.4f} Hz")

Pierson-Moskowitz Spectrum:

def pierson_moskowitz_spectrum(
    frequencies: np.ndarray,
    Hs: float,
    Tp: float = None,
    U19_5: float = None
) -> np.ndarray:
    """
    Calculate Pierson-Moskowitz spectrum (fully developed sea).

    Args:
        frequencies: Frequency array (Hz)
        Hs: Significant wave height (m)
        Tp: Peak period (s) - optional
        U19_5: Wind speed at 19.5m height (m/s) - optional

    Returns:
        Spectral density S(f) (m²/Hz)
    """
    g = 9.81

    if Tp is not None:
        # Use peak period
        fp = 1 / Tp
    elif U19_5 is not None:
        # Calculate from wind speed
        fp = 0.877 * g / (2 * np.pi * U19_5)
    else:
        raise ValueError("Must provide either Tp or U19_5")

    # Phillips constant
    alpha = 0.0081  # For fully developed seas

    # P-M spectrum
    S = alpha * g**2 * (2*np.pi)**(-4) * frequencies**(-5) * \
        np.exp(-1.25 * (frequencies / fp)**(-4))

    return S

# Example
S_PM = pierson_moskowitz_spectrum(freq, Hs=8.5, Tp=12.0)

m0_PM = np.trapz(S_PM, freq)
Hs_PM = 4 * np.sqrt(m0_PM)

print(f"P-M Spectrum Hs: {Hs_PM:.2f} m")

3. Wave Statistics

Spectral Parameters:

def calculate_spectral_parameters(
    S: np.ndarray,
    frequencies: np.ndarray
) -> dict:
    """
    Calculate spectral wave parameters.

    Args:
        S: Wave spectrum (m²/Hz)
        frequencies: Frequency array (Hz)

    Returns:
        Spectral parameters
    """
    # Spectral moments
    m0 = np.trapz(S, frequencies)
    m1 = np.trapz(S * frequencies, frequencies)
    m2 = np.trapz(S * frequencies**2, frequencies)
    m4 = np.trapz(S * frequencies**4, frequencies)

    # Significant wave height
    Hs = 4 * np.sqrt(m0)

    # Mean period
    Tm01 = m0 / m1

    # Zero-crossing period
    Tz = np.sqrt(m0 / m2)

    # Peak period (from spectrum maximum)
    peak_idx = np.argmax(S)
    Tp = 1 / frequencies[peak_idx]

    # Spectral width
    epsilon = np.sqrt(1 - m2**2 / (m0 * m4))

    # Wave steepness
    k_mean = 2 * np.pi / (9.81 * Tz**2 / (2*np.pi))  # Deep water approx
    steepness = k_mean * Hs / 2

    return {
        'm0': m0,
        'm1': m1,
        'm2': m2,
        'm4': m4,
        'Hs': Hs,
        'Tp': Tp,
        'Tz': Tz,
        'Tm01': Tm01,
        'spectral_width': epsilon,
        'steepness': steepness
    }

# Example
params = calculate_spectral_parameters(S, freq)

print(f"Spectral Parameters:")
print(f"  Hs: {params['Hs']:.2f} m")
print(f"  Tp: {params['Tp']:.2f} s")
print(f"  Tz: {params['Tz']:.2f} s")
print(f"  Spectral width: {params['spectral_width']:.3f}")

Wave Height Distribution:

def rayleigh_distribution(
    H: np.ndarray,
    Hs: float
) -> np.ndarray:
    """
    Rayleigh distribution for wave heights in irregular seas.

    P(H) = probability that wave height exceeds H

    Args:
        H: Wave height array (m)
        Hs: Significant wave height (m)

    Returns:
        Exceedance probability
    """
    # Rayleigh parameter
    H_rms = Hs / np.sqrt(2)

    # Exceedance probability
    P = np.exp(-(H / H_rms)**2)

    return P

def significant_wave_statistics(Hs: float) -> dict:
    """
    Calculate wave statistics from Hs using Rayleigh distribution.

    Args:
        Hs: Significant wave height (m)

    Returns:
        Wave statistics
    """
    H_rms = Hs / np.sqrt(2)

    # Various statistical wave heights
    H_mean = H_rms * np.sqrt(np.pi / 2)
    H_1_10 = H_rms * np.sqrt(2 * np.log(10))  # Average of highest 1/10
    H_1_100 = H_rms * np.sqrt(2 * np.log(100))  # Average of highest 1/100
    H_max_1000 = H_rms * np.sqrt(2 * np.log(1000))  # Most probable max in 1000 waves

    return {
        'Hs': Hs,
        'H_mean': H_mean,
        'H_rms': H_rms,
        'H_1_10': H_1_10,
        'H_1_100': H_1_100,
        'H_max_1000': H_max_1000
    }

# Example
stats = significant_wave_statistics(Hs=8.5)

print(f"Wave Statistics (Hs = {stats['Hs']} m):")
print(f"  Mean height: {stats['H_mean']:.2f} m")
print(f"  H_1/10: {stats['H_1_10']:.2f} m")
print(f"  H_1/100: {stats['H_1_100']:.2f} m")
print(f"  H_max (in 1000 waves): {stats['H_max_1000']:.2f} m")

4. Time Series Generation

Generate Irregular Wave Time Series:

def generate_irregular_wave_time_series(
    S: np.ndarray,
    frequencies: np.ndarray,
    duration: float,
    dt: float,
    random_seed: int = None
) -> tuple[np.ndarray, np.ndarray]:
    """
    Generate irregular wave elevation time series from spectrum.

    Args:
        S: Wave spectrum (m²/Hz)
        frequencies: Frequency array (Hz)
        duration: Duration (s)
        dt: Time step (s)
        random_seed: Random seed for reproducibility

    Returns:
        (time, elevation) arrays
    """
    if random_seed is not None:
        np.random.seed(random_seed)

    # Time array
    time = np.arange(0, duration, dt)

    # Initialize elevation
    eta = np.zeros_like(time)

    # Frequency resolution
    df = frequencies[1] - frequencies[0]

    # Generate wave components
    for i, f in enumerate(frequencies):
        if S[i] > 0:
            # Amplitude from spectrum
            amplitude = np.sqrt(2 * S[i] * df)

            # Random phase
            phase = np.random.uniform(0, 2*np.pi)

            # Wave component
            omega = 2 * np.pi * f
            eta += amplitude * np.cos(omega * time + phase)

    return time, eta

# Example: Generate 1 hour of wave data
t, elevation = generate_irregular_wave_time_series(
    S, freq,
    duration=3600,
    dt=0.1,
    random_seed=42
)

# Verify statistics
Hs_timeseries = 4 * np.std(elevation)

print(f"Time Series Statistics:")
print(f"  Target Hs: {params['Hs']:.2f} m")
print(f"  Generated Hs: {Hs_timeseries:.2f} m")
print(f"  Duration: {len(t) * 0.1 / 3600:.2f} hours")

5. Wave Scatter Diagrams

Create Wave Scatter Diagram:

def create_wave_scatter_diagram(
    Hs_bins: np.ndarray,
    Tp_bins: np.ndarray,
    location_data: dict
) -> np.ndarray:
    """
    Create wave scatter diagram (probability table).

    Args:
        Hs_bins: Hs bin edges (m)
        Tp_bins: Tp bin edges (s)
        location_data: Historical wave data or hindcast

    Returns:
        Probability matrix (sum = 1.0)
    """
    # This is typically based on hindcast data
    # Simplified example using lognormal distribution

    n_Hs = len(Hs_bins) - 1
    n_Tp = len(Tp_bins) - 1

    scatter = np.zeros((n_Hs, n_Tp))

    # Simplified: Tp roughly proportional to sqrt(Hs)
    for i in range(n_Hs):
        Hs_mid = (Hs_bins[i] + Hs_bins[i+1]) / 2

        # Expected Tp for this Hs (empirical: Tp ≈ 3.6*sqrt(Hs))
        Tp_expected = 3.6 * np.sqrt(Hs_mid)

        # Hs occurrence (Weibull distribution)
        from scipy.stats import weibull_min
        p_Hs = weibull_min.pdf(Hs_mid, c=2, scale=2.5)

        # Tp distribution given Hs (normal around expected)
        from scipy.stats import norm
        for j in range(n_Tp):
            Tp_mid = (Tp_bins[j] + Tp_bins[j+1]) / 2
            p_Tp_given_Hs = norm.pdf(Tp_mid, loc=Tp_expected, scale=1.5)

            scatter[i, j] = p_Hs * p_Tp_given_Hs

    # Normalize to probabilities
    scatter /= scatter.sum()

    return scatter

# Example
Hs_bins = np.array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10])
Tp_bins = np.array([0, 4, 6, 8, 10, 12, 14, 16])

scatter = create_wave_scatter_diagram(Hs_bins, Tp_bins, {})

# Annual hours
annual_hours = scatter * 8760

print(f"Wave Scatter Diagram:")
print(f"  Total probability: {scatter.sum():.4f}")
print(f"  Most probable sea state: Hs={Hs_bins[np.unravel_index(scatter.argmax(), scatter.shape)[0]]}-{Hs_bins[np.unravel_index(scatter.argmax(), scatter.shape)[0]+1]}m")

6. Extreme Value Analysis

Design Wave from Return Period:

def calculate_extreme_wave_height(
    return_period_years: float,
    Hs_annual_max: np.ndarray = None,
    distribution: str = 'weibull'
) -> dict:
    """
    Calculate design wave height for given return period.

    Args:
        return_period_years: Return period (years)
        Hs_annual_max: Array of annual maximum Hs values
        distribution: 'weibull' or 'gumbel'

    Returns:
        Extreme wave height statistics
    """
    from scipy.stats import weibull_min, gumbel_r

    if Hs_annual_max is None:
        # Example data: 25 years of annual maxima
        np.random.seed(42)
        Hs_annual_max = weibull_min.rvs(c=2, scale=10, size=25)

    # Fit distribution
    if distribution == 'weibull':
        params = weibull_min.fit(Hs_annual_max)
        c, loc, scale = params
        dist = weibull_min(c, loc, scale)
    elif distribution == 'gumbel':
        loc, scale = gumbel_r.fit(Hs_annual_max)
        dist = gumbel_r(loc, scale)
    else:
        raise ValueError("Unknown distribution")

    # Exceedance probability for return period
    exceedance_prob = 1 / return_period_years

    # Extreme value
    Hs_extreme = dist.ppf(1 - exceedance_prob)

    # Confidence intervals (simplified)
    Hs_lower = dist.ppf(1 - exceedance_prob - 0.1)
    Hs_upper = dist.ppf(1 - exceedance_prob + 0.1)

    return {
        'return_period_years': return_period_years,
        'Hs_extreme': Hs_extreme,
        'Hs_lower_bound': Hs_lower,
        'Hs_upper_bound': Hs_upper,
        'distribution': distribution,
        'exceedance_probability': exceedance_prob
    }

# Example: 100-year return period
extreme_100yr = calculate_extreme_wave_height(
    return_period_years=100,
    distribution='weibull'
)

print(f"100-Year Wave:")
print(f"  Hs: {extreme_100yr['Hs_extreme']:.2f} m")
print(f"  Range: {extreme_100yr['Hs_lower_bound']:.2f} - {extreme_100yr['Hs_upper_bound']:.2f} m")

Complete Examples

Example 1: Complete Wave Analysis

def complete_wave_analysis(
    Hs: float,
    Tp: float,
    depth: float,
    duration: float = 3600,
    output_dir: str = 'reports/wave_analysis'
) -> dict:
    """
    Complete wave analysis: spectrum, time series, statistics.

    Args:
        Hs: Significant wave height (m)
        Tp: Peak period (s)
        depth: Water depth (m)
        duration: Time series duration (s)
        output_dir: Output directory

    Returns:
        Complete analysis results
    """
    import plotly.graph_objects as go
    from plotly.subplots import make_subplots
    from pathlib import Path

    output_path = Path(output_dir)
    output_path.mkdir(parents=True, exist_ok=True)

    # 1. Generate spectrum
    freq = np.linspace(0.01, 0.5, 500)
    S = jonswap_spectrum(freq, Hs, Tp)

    # 2. Calculate spectral parameters
    params = calculate_spectral_parameters(S, freq)

    # 3. Generate time series
    t, eta = generate_irregular_wave_time_series(S, freq, duration, dt=0.1)

    # 4. Wave statistics
    wave_stats = significant_wave_statistics(Hs)

    # 5. Regular wave properties (using Tp)
    regular_wave = airy_wave_properties(Hs, Tp, depth)

    # 6. Create visualizations
    fig = make_subplots(
        rows=2, cols=2,
        subplot_titles=(
            'JONSWAP Spectrum',
            'Wave Elevation Time Series',
            'Wave Height Distribution',
            'Wave Steepness'
        )
    )

    # Plot 1: Spectrum
    fig.add_trace(
        go.Scatter(x=freq, y=S, name='S(f)', line=dict(color='blue')),
        row=1, col=1
    )

    # Plot 2: Time series (first 10 minutes)
    t_plot = t[:6000]
    eta_plot = eta[:6000]
    fig.add_trace(
        go.Scatter(x=t_plot, y=eta_plot, name='η(t)', line=dict(width=1)),
        row=1, col=2
    )

    # Plot 3: Wave height distribution
    H_array = np.linspace(0, Hs*2, 100)
    P_exceedance = rayleigh_distribution(H_array, Hs)

    fig.add_trace(
        go.Scatter(
            x=H_array, y=P_exceedance,
            name='Rayleigh',
            line=dict(color='red')
        ),
        row=2, col=1
    )

    # Plot 4: Steepness vs frequency
    steepness_freq = (2*np.pi*freq)**2 / 9.81 * np.sqrt(S)

    fig.add_trace(
        go.Scatter(x=freq, y=steepness_freq, name='Steepness'),
        row=2, col=2
    )

    fig.update_layout(height=800, showlegend=True, title_text=f'Wave Analysis (Hs={Hs}m, Tp={Tp}s)')
    fig.write_html(output_path / 'wave_analysis.html')

    # Export summary
    summary = {
        'input': {
            'Hs': Hs,
            'Tp': Tp,
            'depth': depth
        },
        'spectral_params': params,
        'statistics': wave_stats,
        'regular_wave': regular_wave,
        'time_series': {
            'duration_s': duration,
            'timestep_s': 0.1,
            'points': len(t)
        }
    }

    import json
    with open(output_path / 'wave_summary.json', 'w') as f:
        json.dump(summary, f, indent=2, default=str)

    print(f"✓ Wave analysis complete")
    print(f"  Output: {output_dir}")

    return summary

# Example
analysis = complete_wave_analysis(
    Hs=8.5,
    Tp=12.0,
    depth=1500,
    duration=3600
)

Resources

• Shore Protection Manual: US Army Corps of Engineers
• Ocean Waves and Oscillating Systems: J. Falnes
• Water Wave Mechanics for Engineers and Scientists: R.G. Dean & R.A. Dalrymple
• DNV-RP-C205: Environmental Conditions and Environmental Loads
• ISO 19901-1: Metocean design and operating considerations

Use this skill for all wave analysis in DigitalModel!