# No section

## Realistic simulation of ocean surface using wave spectra

*First presented at the International Conference on
Computer Graphics Theory and Applications (GRAPP)
2006, extended and revised for JVRB*

urn:nbn:de:0009-6-12294

**Abstract**

We present a method to simulate ocean surfaces away from the coast, with correct statistical wave height and direction distributions. By using classical oceanographic parametric wave spectra, our results fit real world measurements, without depending on them. Since wave spectra are independent of the ocean model, Gerstner parametric equations and Fourier transform method can be used with them. Moreover, since they are simple to use and need very few parameters, they allow easy production of ocean surface animations usable in movies and games. We explain how to accurately sample them, to achieve oceanic waves in deep water according to given wind parameters, in a realistic way.

Keywords: Natural Phenomena, Realistic Ocean Waves, Procedural Animation, Parametric Energy Spectra

Subjects: computer simulation, computer animation, computer graphics, realistic computer graphics

The reproductions of many natural phenomena is still an open problem in Computer Graphics. In particular, the realistic simulation of oceanic surfaces continues to be of great interest. Many fields of activity rely on it: virtual reality, movies and games, but also sailing [ CGG01 ] and oceanographic simulations, among others.

While beautiful pictures and animations such as wave refraction [ GlS00 ] [ GM02 ] are produced, there is still a need for oceanic scenes with realistic wave features. Models based on fluid mechanics equations [ EMF02 ] [ MMS04 ] should give the more realistic results, but are inadequate for large oceanic surfaces.

For twenty years, procedural models for wave representation in Computer Graphics have been developed and improved [ Igl04 ]. The Gerstner parametric equations and their Fourier transform rewriting are well known by the oceanographic community. Based on the linear wave theory, they allow representation of natural wave shapes and moves in deep water above the capillary size. As they not describe any net mass transport, they are limited to non-breaking waves and to scenes without violent storms.

Wave height measurement data allow to simulate a given sea state. In the more general case, oceanic surface statistics are taken in account by the use of parametric wave spectra. For given wind parameters, they indicate the distribution of wave energy as a function of frequency. However, none of the existing methods are able to correctly reproduce sea states described by these spectra.

In this paper, we explain how wave spectra are related to oceanic surface energy and we detail a method to deduce wave heights from it. We propose a solution for optimal spectrum sampling, suitable for parametric and Fourier transform models. This method enables the accurate simulation of oceanic surface. Users can easily produce realistic wave animations, with respect to supplied weather conditions.

The structure of this paper is as follows. Section 2 summarizes ocean models for Computer Graphics. In section 3 , we detail our method and give classical parametric spectrum examples. We handle wave statistics in section 4 . We discuss our results in section 5 , andconclude in section 6 .

The *surface elevation*, η, is the vertical displacement
from the water level of the oceanic surface
at rest at
given point and time. The *wave amplitude*, *a*, is the
maximum surface elevation due to this wave. The *wave
height*, *H*, is the vertical distance between a trough and
an adjacent crest. For a single wave, *H = 2a*. The
*wavelength*, λ, is the distance between two successive
crests. The *period*, *T*, is the time interval between the
passage of successive crests through a fixed point.
The *frequency*, *f = 1/T*, is the number of crests that
pass a fixed point in 1 second. The *wave speed* or
*phase speed*, *c = λ/T*, is the speed of wave crests or
troughs. It can also be expressed as *c = ω/κ*, where
*ω = 2πf* is the *angular velocity* or *angular frequency*, and
*κ = 2π/λ* is the *wavenumber*. The *water depth*, *d*, is the
vertical distance between the floor and the water
level at rest.

The concept of *deep water* is wave dependent and is
related to the ratio of the water depth to the wavelength of
the wave. This comes from the transcendental expression

where *g ≈ 9.807 m⋅s ^{-2}
* is the

*standard acceleration of gravity*. Since

*tanh(x) ≈ 1*when

*x > π*, and

*tanh(x) ≈ x*when

*x < π/10*, water is considered deep if

*d/λ > 1/2*and shallow if

*d/λ < 1/20*. This leads to some simplifications in the expressions of the terms given above (see table 1 ).

The base model for ocean waves simulation comes from
the *linear* (or *small-amplitude*) *wave theory* by [
Air45
],
widely used in oceanic engineering, and in Computer
Graphics by [
Pea86
]. It describes waves with a sinusoidal
shape, which corresponds to calm weather conditions.
Considering a location *x _{0}
* lying on a 1D surface at
rest, the elevation

*z = η*at time

*t*due to a wave with wavenumber

*κ*, amplitude

*a*and angular velocity

*ω =*is

When the steepness of the wave increases, its crests
become sharper and its troughs flatter. Therefore, [
FR86
]
used a more realistic description based on *trochoids*, made
by [
vG04
] and by [
Ran63
]. The surface equation is now

where *z _{0}
* is the surface elevation at rest (typically 0).
Note that the important point here is the phase term
difference of

*-π/2*between

*x(x*and

_{0},t)*z(x*, not the use of a particular sine or cosine function. This Lagrangian model describes the trajectory in a vertical plane of a particle

_{0},t)*(x,z)*around its position at rest

*(x*. Summing several waves and extending equation 3 to a 2D surface gives

_{0},z_{0})
where
*= (x,y)* is the horizontal particle position at
time *t*,
*
_{0} = (x_{0},y_{0})* its position at rest,

*= (κ*a wave vector, i.e. a vector with magnitude

_{x},κ_{y})*||*

*|| = κ*and direction of the considered wave propagation and

*=*

*/||*

*||*is the unit vector of . Because of the presence of uniformly distributed random phase terms,

*() ∈ [0,2π[*, this model is also known as

*random waves*. Any set of

*can be taken, which allows the use of adaptive mesh, for optimal surface sampling [ HNC02 ].*

_{0}
Another method for computing equation 4 is the 2D
*inverse discrete Fourier transform*, used by Mastin,
Watterberg and Mareda [
MWM87
], Premože and
Ashikhmin [
PA01
], and Tessendorf [
Tes01
]. The set of
*N x M* complex numbers *F _{n,m}
* is transformed into a set
of complex numbers

*f*by

_{p,q}

where *p ∈ {1,..,N}* and *q ∈{1,..,M} *. This is of
particular interest because the *fast Fourier transform*
(FFT) algorithm is an efficient way to compute these
sum. However, in comparison with the previous
method, its usage is more restrictive. The set of
*
_{0}
*
is defined as a regular grid of

*N x M*points with dimensions

*L*and origin at the center, such that

_{x}x L_{y}*, where*

_{0}= (n_{x0 }L_{x}/N,m_{y0 }L_{y}/M)*n*and

_{x0 }∈ [-N/2,n/2[*m*are integers. We rewrite equation 4 as

_{x0 }∈ [-M/2,M/2[where () and () correspond, respectively, to the imaginary part and real part of the expression, and

Rather than the random phase term , authors commonly use a normally distributed random amplitude (see section 4 ). To comply with equation 5 , the set of wave vectors must be

where *n _{κ}
* and

*m*are defined as

_{κ}*n*and

_{x0 }*m*above. The resulting surface is periodic in all directions and can be used as a tile [ Tes01 ].

_{x0 }To decrease the FFT computation time by a factor two, real numbers can be used instead of complex ones. This is achieved by adding two waves with same amplitude and travelling in opposite directions [ Tes01 ]. Equation 7 becomes

This results in a single stationary wave with time
dependent amplitude. This phenomenon is known as
*standing wave*. The drawback of this method is that the
generated surface is only half the size of the expected one.
The second half is obtained with a rotation of the first one
by 180 degrees around the origin, leading to pattern
duplications.

To be able to use the models discussed so far, we need a
function *a()*. For the parametric model, we also need to
know which set of wave vectors characterizes the surface
we want to generate. The more natural method is to use
records of surface elevation measurement data used by
oceanic engineering [
TG02
]. Measurements are made withaccelerometers mounted on buoys, satellite altimeters,
high frequency radars and others techniques. A spectral
analysis of the collected data is done using the Fourier
transform from time domain to frequency domain.
This decomposition of the surface elevation gives a
histogram of wave frequencies *ω*, with corresponding
amplitudes and possibly directions, known as *frequency
amplitude spectrum*. The same analysis method can
be used with oceanic surface photographs [
Tes01
]
[
TG02
].

An amplitude spectrum can be used to generate a surface with the same characteristics as the related one. However, to get several oceanic surfaces with different sea states, as many spectra are needed. Instead, we rely on parametric wave spectra, which can fit any environment conditions, like wind speed and fetch.

Oceanographers are mostly concerned by wave energy
rather than amplitude. For this purpose, they use another
type of spectrum, *S(ω)*, called *wave frequency energy
spectrum* or *wave spectrum*, which gives information
about wave energy distribution as a continuous function of
frequency. This is a smoothed version of a modified
amplitude spectrum (see below). Directional information
are taken into account with a *directional wave spectrum*:

where *θ* is the direction of the wave and *D(ω,θ)* is
a *directional spreading function*, defined such that

for all values of *ω*. The Pierson-Moskowitz parametric
wave spectrum (see section 3.2 ) has been regularly used
in the field of Computer Graphics since [
MWM87
] as a
direct evaluation of wave amplitudes. But this method
cannot reflect statistical characteristics of the oceanic
surface. First, this is a continuous function, so it needs to
be integrated to be used with a discrete model. Second,
wave spectrum does not give direct information about
wave amplitude but energy.

The mean energy per unit surface area is *E = ρg var(η)*,
where *ρ* is the water density and *var(η)* the *variance* of
the surface elevation. The wave spectrum is defined to be
related to the surface energy by

For a zero mean process, like the sinusoidal wave modelfrom equation 2 , the variance is

where *x _{rms}
* is the

*root mean square*of the values. For a general process, this becomes

*x*

^{2}_{rms}-*. The mean amplitude of the trochoidal wave model from equation 3 is*

^{2}*= -(κa*, which can be reasonably neglected in most cases for practical purposes. So, we apply equation 13 to this model as well for simplicity. The rms of a sinusoidal or trochoidal wave with amplitude

^{2}/2*a*is

As we see the oceanic surface as the sum of many
waves with independent random phases, the total variance
of the surface elevation is the sum of the variance of
each wave. So, taking a finite range of frequencies
*[ω _{min},ω_{max}]* such that it is representative of the whole
energy, i.e.

we divide it in *N* samples of width *Δω _{n}
* and tag the
mean frequency of each of them as

*ω*to get the relation

_{n}

Regarding equations 12 and 15 , this is rewritten as

and so we found the amplitude of each wave with

To take into account the direction of the wave,
we sample the whole directional spectrum *E(ω,θ)*:

So, for each wave with angular frequency *ω* and
direction *θ*, we can now compute the corresponding
amplitude *a(ω,θ)*.

Parametric spectrum models are empirical expressions with user defined parameters that have been found to fit ocean surface elevation measurements. Most common parameters include wind speed and direction and fetch. They are useful to simulate realistic ocean waves without measurement data, but can also be altered to fit specific data.

The base of modern parametric wave spectrum was found by [ Phi57 ] and has the form

where *α = 0.0081* is called the *Phillips constant*.
Taking into account wind speed, [
PM64
] found an
expression that describes a theoretical fully developed sea
with infinite fetch (see figure 1 ). The wind has blown
with no change over a large area for several hours,
and waves growth is almost null. The spectrum is

where *ω _{p}
* is the frequency of the peak of the spectrum,
i.e.

*dS*for

_{PM}/dω = 0*ω*. Usual value for

_{p}*ω*is

_{p}

where *U _{10}
* is the wind speed at a height of 10 meters
above the sea surface.

For a fetch limited sea, where waves continue to grow, the JONSWAP spectrum has been found by [ HBB73 ] to be more appropriate (figure 2 ). This is the previous spectrum with an additional term:

Here, *γ* is a peak enhancement factor, i.e. when same
values are used for *α* and *ω _{p}
* in both spectra, we get

*γ = S*. The parameter

_{J}(ω_{p}/S_{PM}(ω_{p})*σ*controls the width of the peak. Usual parameter values for this spectrum are

where *F* is the fetch in meters,

*γ = 3.3*, but can vary between 1 and 7, and

As the fetch is an important element of the sea state description, we use the JONSWAP spectrum, rather than the Pierson-Moskowitz one.

**Figure 2. JONSWAP frequency spectrum, S_{J}(ω), with
a wind speed of U_{10} = 15m⋅s^{-1}
, and different fetch
values F.**

Other parametric spectrum models can be found in
oceanographic literature, including double peaked
spectrum for swell coming from distant storm. These
spectra are often expressed as a function of frequency
*f = ω/(2π)*. To preserve integral equality, conversion
between spectra is done using substitution:

For the dispersion relation, we use the expression found by [ LH62 ]:

(24)

where *θ _{w}
* is the main direction of the spectrum, usually
the direction of the wind, and

*N(s(ω))*is a normalization factor defined as

and *Γ()* is the *gamma function* (figure 3 ). The
function *s(ω)* controls the sharpness of the directional
spreading and has been defined by [
MTS75
] as

**Figure 3. 3D plot of the JONSWAP spectrum,
E_{J}(ω,θ), with U_{10} = 25m⋅s^{-1}
, F = 100 km and θ_{w}
= 0 rad.**

For a simple use, only wind speed, direction and fetch need to be given to the spectrum. In order to make the spectrum curve fit particular data, like oceanic measurements, the default parameter values can be changed, until the desired curve shape is reached (figure 4 ).

**Figure 4. JONSWAP frequency spectrum,
S_{J}(ω), showing how different parameter values alter
the curve shape (U_{10}=15m⋅s^{-1},F=50km).**

We now look for a selection of wave frequencies and
directions. For the FFT based model, the size and
sample number of the rendered surface determine the set
of wave vectors to be used (see equation 8 ). For the parametric one, there is not such restriction.
We sample a spectrum *E(**)* in the Cartesian wavenumber domain (figure 5). To preserve integral equality, spectrum conversion is done according to substitution rule. First, we convert the frequency spectrum *S _{ω}(ω)* to a wavenumber spectrum

*S*:

_{κ}(κ)
**Figure 5. The same spectrum as in figure 3 , but in
the wave vector domain κ_{x} x κ_{y}
(i.e. E_{J}(
)).**

Second, we convert the wavenumber directional
spectrum *E _{κ,θ}(κ,θ) = S_{κ}(κ)D(ω =*

*θ)*from polar coordinates to Cartesian ones:

we select waves taht are the most representative of the spectrum energy dispersion. For this purpose, we sample the wave spectrum over a limited domain *[κ _{xmin}, κ_{xmax}] x [κ_{ymin}, κ_{ymax}]*. We want this domain to be the most representative of the whole spectrum, i.e. we want it to be as dense as possible. So, by iteration, we found one of the smallest domains that contains roughly 95% of the total energy.

For simplicity, a regular sampling could be used by taking a constant size

**Figure 6. The same spectrum as in figure 3 , but in
the wave vector domain κ_{x} x κ_{y}
(i.e. E_{J}(
)).**

Measurements of oceanic wave statistics have shown
that, to a reasonable accuracy, the surface elevation *η*
follows a *normal distribution*. Attributes of oceanic waves
we simulate are consistent with these observations. If we
look at the oceanic surface as an infinite sum of waves
with infinitely small amplitude, we see that the wave
spectrum is a probability density function of wave
frequencies and directions. Furthermore, the variance of
the surface elevation is finite (equation 12 ). So, by
virtue of the *central limit theorem*, we know that the more
waves we sum, the closer to the normal distribution is the
simulated surface elevation.

Another observed characteristic of oceanic waves is that
their heights follow a *Rayleigh distribution*, as noted by
[
FR86
]. This distribution is used with a parameter
*σ = H _{s}/2*, where

*H*is known as the

_{s}≈ 4*significant height*. From this, it can be shown that the probability that a wave has a larger height than

*H*is

_{s}*exp(-2) ≈ 0.1353*.

Since they use defined amplitudes, the ocean models
presented in section 2 are known as *deterministic*
methods. [
TCC84
] showed that uniformly distributed
random phase terms should be replaced with random
amplitudes. An expression like *a cos(κx _{0} - ωt + )*
becomes

where *r _{1}
* and

*r*are random numbers from the

_{2}*standard normal distribution*, i.e. with a mean of zero and a standard deviation of one,

*R =*follows a Rayleigh distribution and

However, since the values of *x _{0}
* and

*t*depend on arbitrary origins, we can use the uniformly distributed term in place of .

This *non-deterministic* method has better statistical
properties than the previous one. In order to use it with
parametric and FFT methods, we rewrite equations 4
and 6 according to equation 28 .

We made a simple real-time implementation of ocean models. We focused on wave shapes and animations, without considering effects like Fresnel reflectivity and transmissivity, or foam and spray (figure 7 ). Of course, rendering could have been achieved with classical high quality techniques as well. Interactions of objects with ocean surfaces are not specially handled by the Lagrangian model we use, and are beyond the scope of this paper.

Simulation of the interactions of objects with the ocean surface is limited, due to the Langrangian nature of the model. Velocity vectors can easily be computed at each point of the surface grid, and can be used to simulate small object motion. More complex interactions will require additional information and resort to fluid dynamics.

Since our method preserves the main spectrum energy, the global structure of the rendered surface is independent of the number of waves we sum. For both ocean models, the sea state we get is always consistent with the provided wind parameters, which cannot be achieved with other existing methods. Since these parameters are the only ones needed to have full control of the method, there is no need for the user to adjust the resulting surface by trial and error.

As previously noted, the FFT model can be particularly
tedious to use. To catch a reasonable part of the spectrum,
the grid length and the number of samples have to be
carefully chosen, whatever are the desired surface
characteristics. For example, taking 512 samples up to a
frequency of 25 rad⋅s^{-1} requires a grid length of no more
than 25 m. Furthermore, as this leads to a lowest
frequency of about 1.5 rad⋅s^{-1}, the spectrum peak may be
under-sample.

We have tested our implementation with a 3 GHz Pentium 4 PC and a Radeon 9200 graphic board. With the FFT model, we got about 65 fps and 13 fps with, respectively, a grid of 128 x 128 and 256 x 256 samples. With the parametric equations, we got 8 fps with a regular grid of 128 x 128 samples and 50 waves. And taking less than 100 waves leads to poor detailed results. Clearly, only the FFT can reach interactive rate. Although we do not implement an adaptive surface mesh, it seems obvious this could not compete with the FFT speed. Parametric equations should be kept for non-interactive rendering, since they are easiest to use.

**Figure 7.
FFT with a grid of 512 x 512 samples.
Top: grid length is 50 m and wind speed is 5 m⋅s ^{-1}.
Middle: grid length is 450 m and wind speed is 15 m⋅s^{-1}.
Bottom: grid length is 1500 m and wind speed is 25 m⋅s^{-1}.
White square is 1m wide.**

We have presented a method for accurate wave energy spectrum sampling that allows realistic ocean surface synthesis and animation. For given wind parameters, the wave heights and directions are computed such that statistical properties of the resulting surface are correct. Since it does not rely on any ocean model, this method is suitable for Gerstner equations and Fourier transforms.

[Air45]
*Tides and waves*,
Encyclopaedia Metropolitana,
1845,
,
no. 192,
241—396,
B. Fellowes,
London.

[CGG01]
*A new efficient wave model for maritime training simulator*,
Proceedings of the 17th Spring conference on Computer graphics,
p. 202,
2001,
isbn 0-7695-1215-1.

[EMF02]
*Animation and rendering of complex water surfaces*,
ACM Transactions on Graphics,
(2002),
no. 3,
736—744,
ACM Press,
issn 0730-0301.

[FR86]
*A simple model of ocean waves*,
SIGGRAPH Computer Graphics,
(1986),
no. 4,
75—84,
issn 0097-8930.

[GlS00]
*On modeling and rendering ocean scenes*,
Journal of visualisation and computer simulation,
(2000),
no. 1,
27—37,
January,
issn 1049-8907.

[GM02]
*An accurate model of wave refraction over shallow water*,
Computers & Graphics,
(2002),
no. 2,
291—307,
April,
issn 0097-8493.

[HBB73]
*Measurements of wind-wave growth and swell decay during the Joint North
Sea Project (JONSWAP)*,
Ergänzungsheft zur Deutschen Hydrographischen Zeitschrift,
1973,
no. 12,
A8.

[HNC02]
*Interactive Animation of Ocean Waves*,
Symposium on Computer Animation,
pp. 161-166,
2002,
isbn 1-58113-573-4.

[Igl04]
*Computer graphics for water modeling and rendering: a survey*,
Future generation computer systems,
(2004),
no. 8,
1355—1374,
November,
issn 0167-739X.

[LH62]
*The distribution of intervals between zeros of a stationary random
function*,
Philosophical Transactions for the Royal Society of London,
(1962),
no. 1047,
557—599,
Series A, Mathematical and Physical Sciences,
issn 0080-4614.

[MMS04]
*Animation and control of breaking waves*,
Symposium on Computer Animation,
2004,
pp. 315—324,
isbn 3-905673-14-2.

[MTS75]
*Observations of the directional spectrum of ocean waves using a
cloverleaf buoy*,
Jour. of physical oceanography,
(1975),
no. 4,
750—760,
issn 0022-3670.

[MWM87]
*Fourier synthesis of ocean scenes*,
IEEE Computer Graphics and Applications,
(1987),
no. 3,
16—23,
March,
issn 0272-1716.

[PA01]
*Rendering Natural Waters*,
Computer Graphics Forum,
(2001),
no. 4,
189—199,
issn 0167-7055.

[Pea86]
*Modeling waves and surf*,
SIGGRAPH Computer Graphics,
(1986),
no. 4,
65—74,
August,
issn 0097-8930.

[Phi57]
*On the generation of waves by turbulent wind*,
Journal of fluid mechanics,
(1957),
417—445,
issn 0022-1120.

[PM64]
*A Proposed Spectral Form for Fully Developed Wind Seas Based on the Similarity Theory of S. A. Kitaigorodskii*,
Journal of geophysical research,
(1964),
5181—5203,
issn 0022-1406.

[Ran63]
*On the exact form of waves near the surface of deep water*,
Philosophical transactions of the Royal society of London,
(1863),
127—138,
issn 1471-2946.

[TCC84]
*Numerical simulation of a random sea*,
Applied ocean research,
(1984),
no. 2,
118—122,
issn 0141-1187.

[Tes01]
*Simulating ocean water, *
ACM SIGGRAPH course notes,
http://graphics.ucsd.edu/courses/rendering/2005/jdewall/tessendorf.pdf,
2001,
Updated in 2004.

[TG02]
*Ocean waves synthesis and animation using real world information*,
Computers and Graphics,
(2002),
no. 1,
99—108,
February,
issn 0097-8493.

[vG04]
*Theorie der Wellen*,
Abhandlungen der Königlichen Böhmischen Gesellschaft der Wissenschaften,
1804,
1—65.

## License ¶

Any party may pass on this Work by electronic means and make it available for download under the terms and conditions of the Digital Peer Publishing License. The text of the license may be accessed and retrieved at http://www.dipp.nrw.de/lizenzen/dppl/dppl/DPPL_v2_en_06-2004.html.

## Recommended citation ¶

Jocelyn Fréchot, *Realistic simulation of ocean surface using wave spectra*. JVRB - Journal of Virtual Reality and Broadcasting, **4(2007)**, no. 11. (urn:nbn:de:0009-6-12294)

Please provide the exact URL and date of your last visit when citing this article.