Overview

Stochastic sea surface modeling plays a crucial role in ocean engineering, coastal structure design, and wave prediction. This article demonstrates how to generate realistic random ocean surfaces by combining Monte Carlo simulation with two-dimensional wave spectrum models.

Step 1: Define the Wave Energy Spectrum

The wave anergy spectrum describes how wave energy distributes across different frequencies and directions. Common spectrum models include Pierson-Moskowitz (PM) spectrum and JONSWAP spectrum. The following implementation uses the PM spectrum:

function S = compute_pm_spectrum(wavenumber, gravity, wind_speed, peak_factor)
    % Pierson-Moskowitz spectrum implementation
    if peak_factor == 0
        peak_factor = 3.3;
    end
    
    % Calculate spectral density based on PM formula
    denominator = wavenumber.^5 * wind_speed^4;
    numerator = gravity^2;
    exponential_term = exp(-(gravity / (wind_speed^2 * wavenumber)).^4);
    
    S = (numerator ./ denominator) .* exponential_term;
end

Step 2: Generate Complex Amplitudes in 2D Wavenumber Domain

Using the Monte Carlo method, complex amplitudes are generated based on the wave spectrum. Random phases are introduced to simulate the stochastic nature of ocean waves:

function [amplitude, kx, ky] = generate_wave_amplitudes(length_x, length_y, ...
    grid_x, grid_y, gravity, wind_speed, peak_factor, random_seed)
    % Create wavenumber grids
    kx = linspace(0, 2*pi/length_x, grid_x);
    ky = linspace(0, 2*pi/length_y, grid_y);
    [KX, KY] = meshgrid(kx, ky);
    wavenumber_magnitude = sqrt(KX.^2 + KY.^2);
    
    % Initialize random number generator with specified seed
    rng(random_seed, 'twister');
    
    % Generate uniformly distributed random phases
    phase_angle = 2*pi * rand(grid_y, grid_x);
    
    % Compute spectral density values
    spectral_density = compute_pm_spectrum(wavenumber_magnitude, gravity, ...
        wind_speed, peak_factor);
    
    % Construct complex amplitudes from spectral density and random phases
    amplitude = sqrt(spectral_density ./ 2) .* exp(1i * phase_angle);
end

Step 3: Apply 2D Inverse Fourier Transform

Transform the complex amplitudes from wavenumber domain to spatial domain using inverse Fourier transform:

function surface_elevation = compute_spatial_surface(amplitude, length_x, ...
    length_y, grid_x, grid_y)
    % Perform 2D inverse FFT with frequency shifting
    spatial_data = ifft2(ifftshift(amplitude));
    
    % Extract real component and normalize
    surface_elevation = real(spatial_data) / (length_x * length_y * grid_x * grid_y);
end

Step 4: Complete 2D Random Sea Surface Simulation

Integrate all components to generate the stochastic ocean surface:

% Simulation parameters
length_x = 1000;      % Surface length in meters
length_y = 1000;      % Surface width in meters
grid_x = 256;         % Grid points in x direction
grid_y = 256;         % Grid points in y direction
gravity = 9.81;       % Gravitational acceleration (m/s^2)
wind_speed = 10;      % Wind speed at 10m height (m/s)
peak_factor = 3.3;    % Peak enhancement factor for PM spectrum
random_seed = 42;     % Seed for reproducible results

% Generate complex amplitudes in wavenumber domain
[amplitude, kx, ky] = generate_wave_amplitudes(length_x, length_y, ...
    grid_x, grid_y, gravity, wind_speed, peak_factor, random_seed);

% Transform to spatial domain
elevation = compute_spatial_surface(amplitude, length_x, length_y, ...
    grid_x, grid_y);

% Visualization
figure;
imagesc(elevation);
colormap('ocean');
colorbar;
caxis([-1 1] * max(abs(elevation(:))));
title('2D Random Ocean Surface Elevation');
xlabel('X Position (m)');
ylabel('Y Position (m)');
axis equal;

Parameter Considerations

Several key parameters affect the simulation accuracy:

• Grid resolution: Higher grid counts (Nx, Ny) provide better spatial detail but increase computational cost
• Domain size: Larger physical dimensions capture lower-frequency wave components
• Wind speed: Directly influences the energy distribution in the spectrum
• Spectrum model selection: JONSWAP spectrum may be more appropriate for fetch-limited conditions

This methodology provides a foundation for realistic ocean surface generation. Advanced applications may incorporate directional spreading functions, nonlinear wave interactions, and more sophisticated spectrum models such as JONSWAP for better representation of specific sea states.