Simulating Ocean water

4.1 Gerstner Waves

The physical model is to describe the surfacce in terms of the motion of individual points on the surface. To a good approximation, points on the surface go through a circular motion as a wave passes by. If a point on the undisturbed surface is labelled ${\mathbf{x_0} = (x_0, y_0)}$ and the undistured height is ${y_0 = 0}$, then as a singl wave with amplitude A passes by the point on the surface is displaced at time t to

${\mathbf{x} = \mathbf{x_0} - (\mathbf{k}/k)A\sin(\mathbf{k \cdot x_0} - \omega t)}$ (26)

${y = A\cos(\mathbf{k \cdot x_0 - \omega t})}$ (27)

In these expressions, the vector ${\mathbf{k}}$, called the wavevector, is a horizontal vector that points in The direction of travel of the wave, and has magnitude ${k}$ related to the length of athe wave ${(\lambda)}$ by

${k = 2 \pi / \lambda}$

Gerstner waves can be generalized to a more complex profile by summing a set of sin waves. One picks a set of wavevectors ${\mathbf{k_i}}$ , amplitudes ${A_i}$, frequencies ${\omega{_i}}$ and phase ${\phi{_i}}$, for ${i = 1,...,N}$ to get the expression

${\mathbf{x} = \mathbf{x_o} - \sum_{i=1}^N (\mathbf{k_i}/k_i) A_i \sin(\mathbf{K_i \cdot x_0} - \omega{_i} t + \phi{_i}) }$ (29)

${y = \sum_{i=1}^N A\cos(\mathbf{k_i \cdot x_0} - \omega{_i} t + \phi{_i})}$ (30)

4.2 Animating Waves: The Dispersion Relation

The animated behavior of Gerstner waves is determined by the set of frquencies ${\omega{_i}}$ chosen for each component. For water waves, there is a well-known relationship between these frequencies and the magnitude of their corresponding wavevectors ${\mathbf{k}_i}$. In deep water, where the bottom may be ignore, that relationship is

${\omega^2(k) = gk}$ (31)

The parameter ${g}$ is the gravitational constant, nominally ${9.8m/sec^2}$. This dispersion relationship holds for Gerstner waves, and also for FFT-based waves. There are several conditions in which the dispersion relationship is modified. When bottom is relatively shallow compared to the length of the waves, the bottom has a retarding effect on the waves. For a bottom depth ${D}$ below the mean water level, the dispersion relation is

${\omega^2(k) = gk\tanh(kD)}$ (32)

A second situation which modifies the dispersion relation is surface tension. Very small waves, with a wavelength of about 1 cm or less, have an additional term

${\omega^2(k) = gk(1 + k^2L^2)}$ (33)

and the parameter ${L}$ has units of length. Its magnitude is the scale for the surface tension to have effect.

Using these dispersion relationships, it is very difficult to create a sequence of frames of water surface which work for a continuous loop In order to have the sequence repeat after a certain amout of time ${T} for example, it is necessary that all frequencies be multiples of the basic frequency

${\omega_0 } \equiv \frac{2 \pi}{T}$ (34)

However, when the wavevector ${\mathbf{k}}$ are distributed on a regular lattice, it is impossible t arrange the dispersion-generated frequencies to also be on a uniform latticewith spacing ${\omega_o}$. The solution to that is to not use the dispersion frequencies, but instead a set that is close to them For a given wavenumber ${k}$, we use the frequency

${ \bar{\omega}(k) = [[\frac{\omega (k)}{\omega_{0}}]] \omega_0 }$ (35)

4.3 Statistical Wave Model and the Fourier Transform

The FFT based representation of a wave height field expresses the wave height ${h(\mathbf{x}, t)}$ at the horizontal position ${\mathbf{x} = (x, z)}$ as the sum of sinusoids with complex, time dependent amplitudes:

${h(\mathbf{x}, t) = \sum{_k} \tilde{h}(\mathbf{k}, t) \exp{(i\mathbf{k \cdot x})}}$ (36)

where ${t}$ is the time and ${\mathbf{k}}$ is a two-dimensional vector with components ${k = (k_x, k_z)}$, ${k_x = 2 \pi n / L_x, k_z = 2 \pi m / L_z}$ and ${n}$ and ${m}$ are integers with bounds ${-N/2 \leq n \leq N/2 }$ and ${-M/2 \leq m \leq M/2 }$. The fft process generates the height field at discrete points ${\mathbf{x} = (n L_x/N, m L_z /M)}$.

The height amplitude Fourier components, ${\tilde{h} (\mathbf{k}, t)}$, determine the structure of the surface.

The slope vector can be obtained by using more ffts:

${ \epsilon (\mathbf{x}, t) = \Delta h(\mathbf{x}, t) = \sum_k i \mathbf{k} \tilde{h}(\mathbf{k}, t) \exp(i \mathbf{k} \cdot \mathbf{x}) }$ (37)

In terms of this fft representation, the finite difference approach would replace the term ${i \mathbf{k}}$ with terms proportional to

${\exp(i \mathbf{k} \cdot \mathbf{ \triangle x} ) - 1}$ (38)

Statistical analysis of a number of wave-buoy, photographic and radar measurements of the ocean surface demonstrates that the wave height amplitudes ${\tilde{h} (\mathbf{k}, t)}$ are nearly Statistically stationary, independent, gaussian fluctuations with a spatial spectrum denoted by

${P_h(\mathbf{k}) = < \vert \tilde{h^*} (\mathbf{k}, t) \vert{^2} > }$ (39)

A useful model for wind-driven waves larger than capillary waves in a fully developed sea is the phillips spectrum

${ P_h(\mathbf{k}) = A \frac{\exp(-1/(kL)^2)}{k^4} \vert \hat{\mathbf{k}} \cdot \hat{ \omega } \vert{^2}}$ (40)

where ${L = V^2 / g}$ is the largest possible waves arising from a continuous wind of speed ${V}$, ${g}$ is the gravitational constant and ${\hat{ \omega }}$ is the direction of the wind. The cosine term ${ \vert \mathbf{ \hat{k} } \cdot \hat{ \omega } \vert{^2}}$ in the Phillips spectrum eliminates waves that move perperndicular to the wind direction. This model, while relatively simple, has poor convergence properties at high values of the wavenumber ${ \vert \mathbf{k} \vert }$. A simple fix is to suppress waves smaller than a small length ${l \ll L}$, and modify the phillips spectrum by the multiplicative factor

${\exp(-k^2 l^2) }$ (41)

Building a Random Ocean Wave Height Field

Realization of water wave height fields are created from the principles elaborated up to this point: gaussian random numbers with spatial spectra of a prescribed form. This is most efficiently accomplished directly in the fourier domain. The fourier amplitudes of a wave height field can be produced as

${ \tilde{h_0} = \frac{1}{\sqrt{2}} (\xi{_r} + i \xi{_k}) \sqrt{P_h(\mathbf{k})} }$ (42)

where ${\xi{_r} ~ and ~ \xi{_i}}$ are ordinary indepenedent draws from a gaussian random number generator, with mean 0 and standard deviation 1.

Given a dispersion relation ${ \omega (k)}$, the Fourier amplitudes of the wave field realization at time ${t} are

${ \tilde{h} (\mathbf{k}, t) = \tilde{h_0}(\mathbf{k}) \exp \{ i \omega (k) t \} + \tilde{h_0^*}(-\mathbf{k}) \exp \{ -i \omega (k) t \} }$ (43)