Modeling Sea Surfaces

A Tutorial on Fourier Transform Techniques

Curtis D. Mobley
Sequoia Scientific, Inc.
2700 Richards Road, Suite 107
Bellevue, WA 98005

Version 1.0

September 18, 2014

.

1
Abstract

The mathematics relating wind-blown sea surfaces and wave variance spectra are well known,
but seldom well explained, and the devil is in the details. The goal of this tutorial is to fill in those
details as needed for writing computer programs to carry out various calculations so that the results
are physically correct and computationally efficient. The tutorial first shows how measurements
of wave elevation for a random sea surface lead to a variance spectrum, which shows how much
variance (or energy) is contained in waves of different spatial or temporal frequencies. It then shows
how a variance spectrum can be used to generate random realizations of a sea surface. The basic
concepts are introduced for idealized sea surfaces with just one spatial dimension. More realistic
two-dimensional surfaces are then considered. After a thorough investigation of surfaces at a single
instant of time, the last section considers the generation of time-dependent surfaces. An appendix
gives an overview of wave spectra.
Most of this material is available as web pages in the Surfaces, Level 2 section of the Ocean
Optics Web Book, beginning at
http://www.oceanopticsbook.info/view/surfaces/level_2/modeling_sea_surfaces.
The IDL computer codes used to generate many of the figures in this report can be downloaded from
the web book pages. A pdf of this report can be downloaded from the Mobley (2014b) reference in
the references section of the Ocean Optics Web Book at
http://www.oceanopticsbook.info/view/references/publications.

Acknowledgments

The first versions of the IDL computer codes used to generate the figures in this tutorial were
developed in the course of work supported by NASA under contract NNH12CD06C to Curtis
Mobley titled Radiative Transfer Modeling for Improved Ocean Color Remote Sensing. The writing
of this tutorial was supported by my hard-working wife Ann Kruse and many unfunded weekends
and evenings. Converting the hard copy tutorial to web pages for the Ocean Optics Web Book was
supported by NASA under grant NNX14AP66G.

Permission to Reproduce

This report is Copyright c by Curtis D. Mobley, 2014. However, permission is hereby granted to re-
produce this work all or in part, with appropriate acknowledgment, for non-commercial educational
and research purposes.
.

ii
Contents

List of Figures v

List of Tables vi

1 Mathematical Preliminaries 1
1.1 Sinusoidal Wave Representations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.2 Continuous Fourier Transforms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.3 Sampling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
1.4 Discrete Fourier Transforms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.5 Parseval’s Relations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10

2 One-Dimensional, Time-Independent Surfaces 12

2.1 Wave Energy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
2.2 Wave Spectra . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
2.3 Surfaces to Variance Spectra . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
2.4 Frequency Ordering . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
2.5 Variance Spectra to Surfaces . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
2.5.1 Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
2.5.2 Example: A Roundtrip Calculation . . . . . . . . . . . . . . . . . . . . . . . . 26
2.6 Numerical Resolution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28

3 Two-Dimensional, Time-Independent Surfaces 35

3.1 Surfaces to Variance Spectra . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
3.1.1 Example: A Random Sea Surface . . . . . . . . . . . . . . . . . . . . . . . . . 36
3.1.2 Example: A Sea Surface of Crossing Sinusoids . . . . . . . . . . . . . . . . . 38
3.2 Variance Spectra to Surfaces . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
3.2.1 Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
3.2.2 Example: A Two-Dimensional Sea Surface . . . . . . . . . . . . . . . . . . . . 45
3.3 Resampling the FFT Grid for Ray Tracing . . . . . . . . . . . . . . . . . . . . . . . . 47
3.4 Example: Reflectance by Different Surface Models . . . . . . . . . . . . . . . . . . . 50

4 Two-Dimensional, Time-Dependent Surfaces 52

iii
A Wave Spectra 57
A.1 Overview of Wave Spectra . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
A.2 The Pierson-Moskowitz Omnidirectional Gravity Wave Spectrum . . . . . . . . . . . 62
A.3 The Elfouhaily et al. Directional Gravity-Capillary Wave Spectrum . . . . . . . . . . 64
A.3.1 Spreading Functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68

References 70

iv
List of Figures

1.1 Illustration of sampling the harmonics of a given wave and of aliasing . . . . . . . . 7

2.1 A 1-D random sea surface and its Fourier transform and variance spectrum . . . . . 17
2.2 A 1-D surface composed of cosine waves . . . . . . . . . . . . . . . . . . . . . . . . . 20
2.3 A 1-D surface composed of sine waves . . . . . . . . . . . . . . . . . . . . . . . . . . 21
2.4 Example of a 1-D random sea surface generated from the Pierson-Moskowitz spectrum 27
2.5 Example sampling of elevation and slope spectra for N = 1024 . . . . . . . . . . . . 31
2.6 Example sampling of elevation and slope spectra for N = 65536 . . . . . . . . . . . . 32
2.7 Example 1-D surfaces and slopes created using true and adjusted variance spectra . 34

3.1 A 2-D random sea surface and its Fourier transform and two-sided variance spectrum 37
3.2 A 2-D sea surface composed of two crossing sinusoids . . . . . . . . . . . . . . . . . . 40
3.3 Example 2-D sea surface and underlying spectrum and Fourier amplitudes . . . . . . 46
3.4 Mapping a rectangular FFT grid to a hexagonal grid of triangles . . . . . . . . . . . 48
3.5 A sea surface defined by a rectangular FFT grid and by the corresponding grid of
triangular wave facets . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
3.6 Energy reflectance by different sea surface models . . . . . . . . . . . . . . . . . . . . 51

4.1 A time-dependent sequence of 2-D sea surfaces . . . . . . . . . . . . . . . . . . . . . 54

4.2 Perspective view of a 2-D sea surface . . . . . . . . . . . . . . . . . . . . . . . . . . . 56

A.1 The Pierson-Moskowitz variance spectrum as a function of wind speed, k and ω . . . 63

A.2 The Pierson-Moskowitz slope spectrum as a function of wind speed and k . . . . . . 64
A.3 The omnidirectional part of the Elfouhaily et al. (1997) elevation and slope spectra
for fully developed seas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
A.4 The omnidirectional part of the Elfouhaily et al. (1997) elevation and slope spectra
for young to fully developed seas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
A.5 Spreading functions for directional wave spectra . . . . . . . . . . . . . . . . . . . . . 69

v
List of Tables

2.1 Spatial frequencies and wavelengths corresponding to peak variance for the Pierson-
Moskowitz spectrum . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29

4.1 Comparison of Cox-Munk mean square slopes and values for an FFT-generated 2-D
surface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55

A.1 Summary of spectral quantities. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61

vi
CHAPTER 1

Mathematical Preliminaries

“I told her [his younger sister] to start at the beginning and read as far as you can get until you’re
lost. Then start again at the beginning and keep working through until you can understand the
whole book.” –Richard Feynman in Richard Feynman: A Life in Science

*****

Wind-blown, random sea surfaces and the related concept of wave variance spectra are funda-
mental to a wide range of oceanographic problems, which include reflection and transmission of
light, exchange of momentum and energy between winds and currents, loading of structures, and
ship dynamics. Various mathematical techniques for modeling sea surfaces are therefore widely
used in oceanography and ocean engineering. Unfortunately, sea surfaces arise from very compli-
cated hydrodynamics, and it should be no surprise that the mathematics needed to describe them
can also be rather complex. This tutorial seeks to explain a subset of those techniques, namely
the use of Fourier transforms to go back and forth between surface realizations and wave variance
spectra, at a level appropriate for students coming upon the material for the first time.
The material presented here is considered well known and can be found in abbreviated form
in many texts and publications. However, “well known” often means that somewhere, sometime,
someone figured something out, but the details were not given. The implication is that if you can’t
figure out the details for yourself, then you are unworthy to benefit from The Master’s Wisdom. I
do agree that figuring out this material from scratch, as I have had to do, is a worthwhile exercise
that leads to true understanding. However, it is also very time consuming (more than a month of
effort in my case for the material presented here), and time for reinventing wheels is a luxury few
scientists can afford. I would argue that a lucid explanation can also lead to in-depth understanding,
with time left over for application of the knowledge to solving new problems. My goal here is to
explain the material at the level of detail needed to write computer programs. In particular, certain
pesky special cases are considered, as is the often confusing issue of array indexing. The reader can
judge whether or not my explanations are lucid.
We begin with a discussion of a physically and mathematically complicated but very important
topic: how to describe and model wind-blown sea surfaces. Our eventual goal is twofold. The first

1
goal is to show how to describe a time-dependent sea surface with waves of all scales propagating in
any direction. The second is to show how to generate time-dependent random realizations of such
surfaces so that the generated surface is close to what would be observed in nature. Such generated
surfaces can be used, for example, as the basis for Monte Carlo simulations of the optical reflection
and transmission properties of sea surfaces, and even to simulate the visual appearance of water
surfaces for use in computer-generated movie scenes.
Real sea surfaces exhibit extremely complicated shapes involving waves of many sizes traveling
in all directions. Those waves arise from various nonlinear physical processes that transfer energy
from the wind to the water, and between waves of different sizes. It should not come as a surprise
that the mathematics needed to describe these processes and the resulting sea surfaces is also
√
complex, both in the psychological sense and in the −1 sense.
If you come to optical oceanography from a background in physics, electrical engineering, or
numerical analysis, then you are probably already familiar with Fourier transforms, variance (or
energy or power) spectra, sampling theory, and the other tools that we will need. However, if you
began life as a biological oceanographer, then even the idea of complex numbers may strike fear
in your nonmathematical heart. We will therefore begin slowly and study various idealized cases
in detail before proceeding to real ocean surfaces. We will progress from the simplest possible
physical context of a snapshot (i.e., the surface at one instant of time) of a one-dimensional (1-D)
sea surface, then a time-independent two-dimensional (2-D) surface, and finally reach the goal of
simulating time-dependent 2-D surfaces as are found in nature.

1.1 Sinusoidal Wave Representations

We describe the sea surface using a Cartesian coordinate system. x is a unit vector chosen to point
in the downwind direction, z is directed upward (from the water to the air) normal to the mean
sea surface at z = 0, and y = z × x is in the cross-wind direction. z(x, y, t) is the surface elevation
in meters at spatial location (x, y) and time t.
Fourier proved in the early 1800’s that, under very general mathematical conditions, an arbitrary
function can be written as a sum of sines and cosines of different amplitudes and wavelengths. The
sea surface elevation z(x, y, t) is usually such a function. An exception is a sea surface with waves
curling over and breaking on a beach. There can then be multiple air-water surfaces for a given
(x, y, t): the point on the surface that sees the full sky, the point above a surfer’s head as the wave
curls over her, and the point on the water that supports the surfboard. However, for a sea surface
without breaking waves, there is only one air-water surface for a given location and time, which
can be described as a sum of sinusoids.
We first define the terms and notation needed to describe propagating sinusoidal waves:

• A is the amplitude (in units of meters) of the wave. This is one-half the distance from the
trough to the crest of a sinusoidal wave. Oceanographers often talk about the wave height, h,
which is the distance from a wave trough to a wave crest. Wave height is what most people
visualize when they talk about how “big” a wave is.

• T is the temporal period (seconds), usually called just the period. This is the time it takes for
one wave crest to pass by a fixed point, i.e. for the argument of the sinusoid to go through

2
2π radians.

• f = 1/T is the temporal frequency (1/seconds), usually called just the frequency. This is how
many waves pass by a fixed point per second.

• ω = 2π/T = 2πf is the angular frequency (radians/second). This is the number of waves
passing a fixed point in 2π seconds.

• Λ is the wavelength, or spatial period, (meters). This is the distance from one wave crest to
the next.

• ν = 1/Λ is the wavenumber or spatial frequency (1/meters). This is the number of wavelengths
per meter.

• k = 2π/Λ = 2πν is the angular wavenumber or angular spatial frequency (radians/meter).

This is the number of wavelengths in 2π meters.

• φ is the phase (nondimensional). This fixes the location of a wave crest relative to the origin
of a coordinate system.

There is no uniformity in the literature. People often refer to “the frequency” without specifying
whether they mean f or ω, and “wavenumber” or “spatial frequency” may mean either ν or k.
People use ν, ν̃, or σ for the wavenumber, some use σ for ω, and so on. You just have to figure
it out on a case by case basis. For pedagogic purposes, we will start with wavenumber ν as the
measure of spatial frequency, and then switch to the more common angular spatial frequency k.
We will use ω as the measure of temporal frequency.
Now consider a one-dimensional sea surface with elevation z(x, t) over a spatial region of length
L, e.g. 0 ≤ x ≤ L or −L/2 ≤ x ≤ L/2. We can write z(x, t) as a sum of sinusoids, the nth one of
which is
2πnx 2πnt
zn (x, t) = An cos + φn ± , (1.1)
L T
where n = 0, 1, 2, ... is simply an index for which sinusoid is being considered. (Note that we could
just as well write a sine here.) For n = 0 this reduces to a constant offset z0 (x, t) = A0 cos(φ0 ),
which is usually taken to be mean sea level and set to a reference value of 0 via setting A0 = 0.
For the moment, we take t = 0, i.e., we have a snapshot of the ocean surface at time zero. This
sinusoid has a wavelength of Λn = L/n, i.e., the cosine returns to the same value after a distance
of x = L/n. It is a pure cosine in the chosen coordinate system if the phase is 0 or an even integer
multiple of ±π. If the phase is an odd integer multiple of ±π/2 it is a sine. If time is included, the
wave is propagating in the +x direction if the time term is −2πnt/T ; the cosine returns to its initial
value after a time of t = T /n. A wave propagating in the −x direction is described by a +2πnt/T
term. The physical angular frequency ωn ≡ 2πn/T is always positive, but mathematically we can
write just +ωn t in the equation and then view a wave propagating in the +x direction as having a
negative temporal frequency.
With this interpretation of ωn t, Eq. 1.1 is conveniently written as

zn (x, t) = An cos(kn x + φn + ωn t) , (1.2)

3
where kn = 2πn/L = 2π/Λn is the angular wavenumber of the nth wave of wavelength Λn = L/n.
Likewise, ωn = 2πn/T = 2π/Tn = is the angular frequency of the nth wave with period Tn = T /n.
It may be intuitively easier to think of wavelengths per meter than of wavelengths per 2π meters,
but the convenience of writing k rather than 2π/Λ or 2πν makes k rather than ν the spatial
frequency variable used in most publications. A similar comment holds for ω vs. 2π/T = 2πf , so
ω rather than f is the common choice for the temporal frequency variable.
If we take t = 0, or just combine the time term with the phase, expanding the cosine in Eq.
(1.2) gives an equivalent form

zn (x) = an cos(kn x) + bn sin(kn x) , (1.3)

where an = An cos φn and bn = −An sin φn , n = 1, 2, ..., a0 = A0 , and b0 = 0. This equation can be
written in terms of complex numbers:

zn (x) = c+n eikn x + c−n e−ikn x , (1.4)

√
where i = −1 and

c+n = (an − ibn )/2, c−n = (an + ibn )/2, and c0 = a0 /2 . (1.5)

Recalling that e±iθ = cos θ ± i sin θ, we see that the imaginary parts of Eq. (1.4) add to zero,
so that zn (x) is still a real variable even though it is written in complex form. Also note that
c∗+n = c−n , where c∗ denotes the complex conjugate. (Complex conjugation means replace i by −i
in all terms.) Although two complex numbers in general contain four independent real numbers,
these c+n and c−n pairs contain only two independent numbers, an and bn .
We thus have three ways to describe a sinusoidal wave: (1) The cosine of Eq. (1.2) defined by
an amplitude and a phase; (2) the cosine and sine of Eq. (1.3) defined by two amplitudes, or (3)
the complex exponentials of Eq. (1.4) defined by two complex numbers containing two independent
real numbers. These equations all give the same zn (x), so we are free to choose whichever form
is mathematically convenient for the problem at hand. For visualization of sea surfaces, the real
forms (1.1) or (1.3) are convenient, but the calculations below are most conveniently carried out
using complex exponentials.
Returning now to Fourier, linear wave theory says that we can write the sea surface elevation
as a Fourier series
∞ ∞
X a0 X
z(x) = zn (x) = + [an cos(kn x) + bn sin(kn x)] , (1.6)
2
n=0 n=1

which in conjunction with Eq. (1.4) can be rewritten as

∞
X
z(x) = cn eikn x . (1.7)
n=−∞

(As previously noted, the a0 or c0 term is usually set to 0.) This equation is the mathematical
heart of our subsequent descriptions sea surfaces.
Although the physics that leads to a given sea surface in general involves non-linear interactions
between waves of different wavelengths and periods, we can write the shape of the resulting sea

4
surface as a linear sum of sinusoids of different frequencies. Owing to the orthogonality of sines and
cosines for different n values as defined above, these wave components zn (x) are independent of each
R 2π
other. (Orthogonality means that 0 cos(mx) cos(nx)dx = 0 if m 6= n, etc.) Thus a description
of a time-dependent surface based on an expansion like that of Eq. (1.6) (including the time-
dependent terms as seen in Eq. (1.2)) cannot be used to predict the nonlinear evolution of a sea
surface from a given initial condition. To do that, you must return to the world of hydrodynamics
in all of its nonlinear glory. I lived there as a graduate student and, trust me, it isn’t a pretty
place. For our present purposes, we only want to model the shapes of sea surfaces, not predict their
hydrodynamic development from an initial state, so linear wave representations are adequate.

1.2 Continuous Fourier Transforms

The Fourier transform is the fundamental mathematical tool for our subsequent study of sea sur-
faces. There are many texts on Fourier transforms. Bracewell (1986) is a standard reference, and
excellent sets of lecture notes and videos of lectures can be found on the web. We will therefore
state without proof the results needed for our problem.
Given a real function z(x) of a continuous variable x, the Fourier transform ẑ(ν) of z(x) is
defined as Z ∞
ẑ(ν) ≡ F{z(x)} ≡ z(x)e−i2πνx dx . (1.8)
−∞
The inverse Fourier transform is given by
Z ∞
z(x) = F −1 {ẑ(ν)} ≡ ẑ(ν)e+i2πνx dν . (1.9)
−∞

It can be shown that if we insert the ẑ(ν) integral of Eq. (1.8) into Eq. (1.9), we recover z(x). (This
is known as Fourier’s integral theorem.) Equations (1.8) and (1.9) are termed a Fourier transform
pair.
Understanding the units of a Fourier transform is important. In our case, z(x) is sea surface
elevation, and both z and x have units of meters. Equation (1.8) shows that ẑ(ν) thus has units
of m2 . The variance of z also has units of m2 , which gives the first hint at a profound connection
between the variance of a physical quantity and its Fourier transform. The units of m2 in the
Fourier transform also can be rewritten as m/(1/m), which is units of z divided by units of spatial
frequency ν. The transform ẑ(ν) can thus be interpreted as showing “how much of z there is
per unit frequency interval.” The inverse transform then has units of (z over frequency) times
frequency, which returns the original units of z.
A Fourier transform is a spectral density function. A spectral density function is a function whose
integral over a given frequency interval gives the variance in the physical quantity contributed by
the frequencies in the integration interval. Density functions are rather peculiar mathematical crea-
tures compared to point functions, which simply give the value of something at a given value of the
independent variable (e.g. the temperature as a function of location x and time t). The blackbody
radiation function is another example of a spectral density function. The blackbody function shows
how much energy is emitted (at a given temperature) per unit frequency interval of the emitted
electromagnetic radiation. The blackbody function is discussed on the Ocean Optics Web Book at
http://www.oceanopticsbook.info/view/light_and_radiometry/level_2/a_common_misconception.

5
If you are not familiar with the distinction between point and density functions, especially regarding
how to change variables in density functions, you may wish to take a look at that discussion.
The Fourier transfer definitions above with the 2π in the exponents are those of the “Stanford
school” of Bracewell (1986) and Goodman (1996). You will see others in the literature. If we use
k = 2πν as the frequency variable, then Eqs. (1.8) and (1.9) become
Z ∞
ẑ(k) ≡ F{z(x)} ≡ z(x)e−ikx dx (1.10)
−∞

and Z ∞
1
z(x) = F −1 {ẑ(k)} ≡ ẑ(k)e+ikx dk . (1.11)
2π −∞
This reappearance of the 2π in the second equation is required so that the inverse transform of
the transform gets you back to where you started. In the e±ikx version, some people put the 1/2π
√
in front of the other integral, and some put a 1/ 2π in front of each integral. Some people (e.g.
Press et al., 1992) put the +i in Eq. (1.8) and the −i in Eq. (1.9). The choice of which sign
to use on the i and where to put the 2π is almost a religion—most people stay with what they
first learned, are convinced of the superiority of their definition, and are willing to die rather than
change. Fortunately, it doesn’t matter which definitions you use, so long as you are consistent in
how the transform pair is defined so that you get back to where you started if you inverse transform
a transform, or vice versa.
In our work, z(x) is the sea surface elevation, which is a real number. However, even though
z(x) is real, ẑ(ν) (or ẑ(k)) is complex. Expanding the complex exponential in Eq. (1.8) as the sum
of a cosine and a sine, it is easy to see that ẑ ∗ (ν) = ẑ(−ν). Such functions are called Hermitian.
A real function has a Hermitian Fourier transform. Conversely, if a function is Hermitian, it
has a real inverse Fourier transform. This will be an important constraint below when we wish
to generate random realizations of a sea surface by computing the inverse Fourier transform of a
complex function: we will have to conjure up a Hermitian function so that we end up with a real
sea surface.

1.3 Sampling
We next consider the implications of sampling the sea surface at a finite number of discrete locations
or times. It is always the case when we want to measure the sea surface elevation that we either
make measurements at discrete locations xr , r = 0, 1, ..., N − 1 at a given time, or at discrete times
tr , r = 0, 1, ..., N − 1 at a given location.
For a specific example, consider a region of sea surface L = 10 m long where we take N = 8
evenly spaced samples. This gives a measurement every ∆x = L/N = 0.125 m at a given time. We
are going to describe this surface as a sum of sinusoids. Figure 1.1 illustrates this, with the first
few sinusoids shown as cosines. The cosines in the figure all have the same amplitude and each is
offset vertically for ease of viewing. The dots show the sampled values for each cosine component.
The longest wave that can fit into a region of length L has a wavelength of Λ = L. This is
called the fundamental wave or first harmonic. This wave is shown by the blue n = 1 curve in the
figure. Note that the n = 4 cosine in the figure has a period of 2∆x. The smallest wave than can be
captured correctly in a sampling scheme has a period of twice the sampling interval. This called the

6
Figure 1.1: Illustration of sampling the harmonics of a given wave and of aliasing. The vertical
dotted lines show the sampling locations, and the dots show the sampled values for the various
cosine components.

two-point wave or two-point oscillation. Note that a sine term with period 2∆x would be sampled
at arguments of 0, π, 2π, ... where the sine is zero and is thus not detectable. The spatial frequency
1/(2∆x) (or temporal frequency 1/(2∆t) if sampling in time at a given location) is called the spatial
(or temporal) Nyquist frequency. Waves with higher frequencies than the Nyquist frequency are
still sampled, as illustrated by the dots on the wave with n = 10. Note, however, that the pattern
of sampled points for the n = 10 wave is exactly the same as for the n = 2 wave. If all we have
are the measured points, we cannot distinguish between the n = 10 and the n = 2 waves in this
example. This illustrates the phenomenon of aliasing. In general, a wave with a frequency greater
than the Nyquist frequency (i.e., a wavelength or period less that twice the sampling interval) is
still sampled, but it appears as though the high frequency wave is a wave with a frequency less than
the Nyquist frequency. The information about the high frequency (short wavelength) wave is added
to or aliased into a lower frequency (longer wavelength) wave, thereby giving incorrect information
about the lower frequency wave. Aliasing places a severe constraint on any sampling scheme. We
can correctly sample only waves with wavelengths that lie between the size of the spatial region,
L, and twice the size of the sampling interval. (In the temporal setting, the limits are the length
of sampling time and twice the temporal sampling interval.)
The relation between a sampled frequency fs greater than the Nyquist frequency, the sampling
rate or Nyquist frequency fNy , and the frequency fa receiving the aliased signal is

fa = |mfNy − fs |

where m is the closest integer multiple of the sampling rate to the signal being sampled. In the
example of Fig. 1.1, fs = fn=10 = 10 waves per 10 meters = 1 m−1 , fNy = 0.4, and either m = 2
or m = 3 gives fa = |2 × 0.4 − 1| = 0.2 (or fa = |3 × 0.4 − 1| = 0.2). Thus the n = 10 wave is

7
aliased into the n = 2 wave, just as seen in the figure. Note that many other frequencies are also
aliased into this frequency. For example, a wave with fs = 1.4, corresponding to n = 14, also gives
fa = 0.2 with m = 4, and so on.

1.4 Discrete Fourier Transforms

Now suppose that we have sampled the sea surface z(x) at a set of N evenly spaced points xr , r =
0, 1, ..., N − 1, in a region of size L; xr = r∆x = rL/N . We want to describe this sampled sea
surface z(xr ) as a sum of sinusoids. In general, these N values can be represented as a sum of a
constant term, N/2 cosine terms, and N/2 − 1 sine terms (recall that there is no two-point sine
term):
N/2
a0 X
z(xr ) = + [au cos(ku xr ) + bu sin(ku xr )] , (1.12)
2
u=1

with bN/2 ≡ 0. Note that the sum runs from the fundamental frequency (n = 1, k1 = 2π/L)
through the Nyquist frequency (n = N/2, kN/2 = 2π/(2∆x)), with only a cosine term for the
two-point wave. This sum is equivalent to
N/2
X
z(xr ) = cu eiku xr , (1.13)
u=−N/2−1

which also contains a total of N independent real and imaginary parts of the cn coefficients. (Recall
from Eq. (1.5) that c−k = c∗k , so these coefficients are not independent for k and −k pairs.) These
equations are the discrete-variable forms of Eqs. (1.6) and (1.7).
Comment on notation: It is common to use the letters i, j, k for dummy summation indices.
√
However, I’ve already used i for −1 and k for angular wavenumber, so the preceding equations
would be hopelessly confusing if I reused i and k for indices. I will therefore use r and s for indices
on spatial variables, e.g., (xr , ys ), and u and v for indices on frequency variables, e.g., ku or νv . n
and m also will be used as needed for dummy indices.
We now have a finite number N of discrete samples z(xr ) of the sea surface, so we need a
discrete form of the Fourier transform. The discrete Fourier transform of z(xr ) is defined as
N −1
1 X
ẑ(νu ) ≡ z(xr )e−i2πνu xr .
N
r=0

Recalling that νu = u/L = u/(N ∆x) and xr = r∆x = rL/N gives νu xr = ru/N . It is also common
to write z(xr ) = z(r) and ẑ(νu ) = ẑ(u), in which case the previous equation becomes
N −1
1 X
ẑ(u) = z(r)e−i2πru/N . (1.14)
N
r=0

The corresponding inverse discrete Fourier transform is given by

N
X −1
z(r) = ẑ(u)e+i2πru/N . (1.15)
u=0

8
As always, there are competing definitions. Equations (1.14) and (1.15) are used in Bracewell
and the IDL computer language. Numerical Recipes (Press et al., 1992) interchanges the i and −i.
Matlab puts the 1/N factor on the inverse transform. Also, Matlab does not support array indices
of 0, so the summation indices are shifted from 0 to N −1 to 1 to N , with a corresponding shift from
ru to (r − 1)(u − 1) in the exponentials. As I said at the start, the devil is in the details, and details
like where to put the 1/N factor and differences in array indexing in different computer languages
can cause great misery when it comes time to actually write computer programs or compare results
computed by different canned subroutines.
The sums in the last two equations require computing complex exponentials (i.e., sines and
cosines), multiplying by the corresponding values of z(r) or ẑ(u) and adding up the results. These
equations can be evaluated for any value of N . The number of computations required to do this is
of order N 2 . The computation time thus increases very rapidly for large N .
However, a classic paper by Cooley and Tukey (1965) showed how these sums can be computed
in order N log2 N computations, if N is a power of 2. Their technique is now called the Fast Fourier
Transform or FFT. The difference in computer time becomes enormous for large N . For example,
if N = 212 = 4096, then N log2 N = 4096 × 12, and the difference in computation time is a factor of
N 2 /(N log2 N ) ≈ 341. Thus in the case of N = 4096, a roughly 6 minute computer run becomes a
1 second run. The development (or, perhaps, reinvention, since the basic idea goes all the way back
to Gauss) of the FFT was a major advance in numerical analysis, which enables the computations
on the following pages to be performed extremely efficiently. Subroutines for computing FFTs and
inverse FFTs are available in all computer languages commonly used in science (Fortran, C, IDL,
Matlab, etc). Fortunately we do not need to concern ourselves here with the details of how the
FFT algorithm actually works, any more than we need to worry about how a canned subroutine
actually computes the cosine of a number. If you are interested in how the FFT works, a web
search will yield many detailed explanations.
The one-dimensional (1-D) equations seen above are easily extended to two or more dimensions.
For two dimensions (x, y), we can sample a region of size Lx by Ly meters over Nx points in the x
direction and Ny points in the y direction, with Nx and Ny both powers of 2 so we can use FFTs.
Equations (1.12) and (1.13) then become

Nx /2 Ny /2
a0 X X
z(xr , ys ) = + [au,v cos(ku xr + kv ys ) + bu,v sin(ku xr + kv ys )]
2
u=1 v=1
Nx /2 Ny /2
X X
= cu,v ei(ku xr +kv ys ) . (1.16)
u=−Nx /2−1 v=−Ny /2−1

The corresponding 2-D FFT pair is

x −1 N
NX y −1
1 X
ẑ(u, v) = z(r, s)e−i2π(ru/Nx +sv/Ny ) (1.17)
Nx Ny
r=0 s=0

and
x −1 N
NX y −1
X
z(r, s) = ẑ(u, v)e+i2π(ru/Nx +sv/Ny ) . (1.18)
u=0 v=0

9
We emphasized above that the continuous function ẑ(ν) defined by the continuous Fourier
transform (1.8) is a density function with units of z(x) per spatial frequency interval, e.g., m/(1/m)
if z is sea surface height in meters. However, the discrete functions ẑ(u) and ẑ(u, v) defined by the
discrete Fourier transforms (1.14) and (1.17) have the same units as z(r) and z(r, s). The discrete
Fourier transform is a point function that shows how much of z(r) is contained in a finite frequency
interval ∆ν = 1/L centered at frequency νu = u/L. Discrete Fourier transforms therefore convert
point functions z(r) and z(r, s) to point functions ẑ(u) and ẑ(u, v).
It will be important below to keep notational track of continuous vs. discrete versions of various
functions. Thus S(k) will denote a continuous function of k, S(k = ku ) will denote the continuous
function evaluated at the discrete value ku , and S(ku ) = S(u) will denote a discrete point function.
Keep in mind that S(k = ku ) and S(ku ) have different units.
The differences in units between continuous and discrete Fourier amplitudes sometimes makes
it tricky to make the transition between discrete and continuous versions of the same quantity. In
particular, it will be necessary to explicitly include the frequency intervals in some of the later cal-
culations that involve both continuous and discrete variables. For example, if we have a continuous
density function and we need to convert to a corresponding discrete function, we must multiple the
continuous function by the finite frequency interval, e.g. ẑ(u) = ẑ(ν = νu )∆ν. Conversely, if we
have discrete amplitudes ẑ(u) and we wish to estimate the continuous spectral density ẑ(ν), then
we must divide by the frequency interval: ẑ(ν) = ẑ(u)/∆ν. (If you are an optical oceanographer
familiar with the scattering phase function, you can find an analogous situation in the estimation of
the scattering phase function from measurements of scattered light. The scattering phase function
is a measure of how much light energy is scattered from an incident direction into a particular
direction, per unit solid angle; it therefore has units of 1/steradian. If you measure the scattered
light using an instrument with a finite solid angle ∆Ω, then you get the total amount of energy
scattered into the solid angle ∆Ω. To estimate the phase function from this measurement, you
must divide the measured value by the solid angle of the instrument; this gets you back to units of
1/steradian.)

1.5 Parseval’s Relations

The physical and spectral variables of a continuous Fourier transform pair satisfy
Z ∞ Z ∞
2
|z(x)| dx = |ẑ(ν)|2 dν , (1.19)
−∞ −∞

which is known as Parseval’s relation. The corresponding equation for the discrete Fourier transform
pair defined by Eqs (1.14) and (1.15) is
N
X −1 N
X −1
2
|z(r)| = N |ẑ(u)|2 . (1.20)
r=0 u=0

For complex amplitudes, |ẑ|2 = ẑ ẑ ∗ . The extension to the two-dimensional case is straightforward:

x −1 N
NX y −1 x −1 N
NX y −1
X X
|z(r, s)|2 = Nx Ny |ẑ(u, v)|2 . (1.21)
r=0 s=0 u=0 v=0

10
The discrete forms of Parseval’s relations provide important checks on numerical calculations.
For example, it is easy to misplace factors of N, which appear in different places depending on the
exact form used for the definition of the discrete transforms.
You can take a class in Fourier transforms and prove many more pretty theorems about their
properties. However, we have now assembled the mathematical tools needed to describe wind-blown
sea surfaces, and so we return to oceanography.

11
CHAPTER 2

One-Dimensional, Time-Independent Surfaces

We now begin our study of sea surfaces and their Fourier transforms. For simplicity, we will
introduce the basic ideas in a one-dimensional (1-D) geometry. This corresponds to making mea-
surements of the sea surface along a line of points at a given time, which yields a set of surface
elevations z(xr ) = z(r), r = 0, 1, ..., N . The Fourier transform of z(xr ) will be in terms of the spatial
frequency νu or the angular spatial frequency ku . Conversely, we could make measurements at a
single spatial point over a set of times. This yields a set of elevations z(tr ). The frequency variable
for the Fourier transform of z(tr ) would be either the temporal frequency fu or angular temporal
frequency ωu . The math is the same in either case. Our primary interest is generating sea surfaces
at a given time, so the discussion (until Chapter 4) will be in terms of spatial sampling at a given
time, which will be taken as t = 0 and not shown.

2.1 Wave Energy

Consider a sinusoidal wave of amplitude A and wavenumber ν. It can be shown (e.g. Preisendorfer,
1976, Vol. VI, page 72) that the total energy per unit horizontal area of sea surface of this wave,
averaged over a wavelength, is

energy 1 1
= ρgA2 + τ ν 2 A2 ,
area 2 4
where

ρ ≈ 103 kg m−3 is the density of water,

g ≈ 9.8 m s−2 is the acceleration of gravity,
τ ≈ 0.072 N m−1 (at 25 deg C) is the surface tension of water.

The 12 ρgA2 term is the sum of the kinetic and gravitational potential energy of the wave. (The
potential energy is relative to the mean sea surface at z = 0.) The 14 τ ν 2 A2 term is the energy
required to stretch the level surface into a sinusoid, working against surface tension. It is easy to
see that the units of these terms are J m−2 .

12
The variance of a sinusoidal surface is the average over one wavelength Λ of the surface elevation
squared (assuming that the mean surface is at zero):
Z Λ 2
1 2πx 1
var{z} ≡ A sin dx = A2 . (2.1)
Λ 0 Λ 2

Thus the energy of a sinusoidal wave also can be written in terms of its variance:

energy 1 2
= ρg + τ ν var{z} . (2.2)
area 2

2.2 Wave Spectra

Equation (2.1) shows that the variance of the surface elevation is proportional to the amplitude
squared. Equation (2.2) shows that the energy per unit horizontal area of a sinusoidal wave of
a given amplitude and spatial frequency is also proportional to the amplitude squared, or to the
variance of the surface elevation. For a sea surface containing waves of many different amplitudes
and spatial frequencies, the total variance of the sea surface is the sum of the variances of in
the individual waves (the variance of a sum of independent random variables is is the sum of the
variances of the individual variables). Likewise, the total energy of the waves is the sum of the
energies of the individual waves .
P −1
For a discrete sample of a zero-mean sea surface, the variance is N1 N 2
r=0 |z(r)| , which ac-
PN −1 2
cording to Parseval’s relation (1.20) equals u=0 |ẑ(u)| . Most importantly, we have seen that
the discrete Fourier transform converts the surface elevations z(r) to amplitudes ẑ(u) that show,
loosely speaking, how much of z(r) there is in each finite frequency band centered at νu . We have
just seen that wave energy is proportional to wave variance. This line of reasoning led Schus-
ter (1898) to the seminal observation that Fourier transforms can be used to decompose the total
variance contained in a signal into the variance contained in each frequency. Schuster’s original
interest was in searching for what he call “hidden periodicities” in weather phenomena (as opposed
to obvious periodicities such as daily or seasonal cycles). He soon applied the technique to looking
for periodicities in earthquakes, sunspots, and other phenomena. Today, the computational speed
of the FFT allows his method of analysis to be widely applied in science, engineering, economics
and other fields where it is desired to know the energy or power of a signal as a function of spatial
or temporal frequency. Regardless of what physical quantity is under consideration, the essence of
a spectrum is that it gives the distribution of variance in that quantity as a function of frequency.
In this tutorial we are concerned with generating sea surfaces consistent with a chosen wave
spectrum. For that application, it is the roughness of the surface—that is, the variance of the
surface elevation at different spatial scales—that is of interest, not the energy contained in the
waves. It is thus the variance spectrum that is of interest. If we wanted to model the forces on an
oil derrick or ship, then we would be interested in the energy of the waves and would work with
the energy spectrum.
The function |ẑ(u)|2 , u = 0, 1, ..., N − 1 is called the discrete variance spectrum. In our context,
it is the variance in sea surface elevations that is described by |ẑ(u)|2 . It common to call |ẑ(u)|2
the discrete energy or power spectrum. (Publications often make a cryptic comment about |ẑ(u)|2
giving energy “when suitably normalized,” but they never seem to show the variance-to-energy

13
conversion factor (ρg + 21 τ ν 2 ).) This may permissible if it is relative values that are of interest, e.g.
this frequency has ten times more energy than that frequency. In those cases, it is only the shape of
the spectrum that is important. However, we are here interested in the variance of the sea surface,
so I will correctly refer to |ẑ(u)|2 as the discrete variance spectrum, whose units are m2 (the units
of both |z(r)|2 and |ẑ(u)|2 ). If you really do want to get the energy per unit area of the sea surface
contained in the frequency interval ∆ν at νu , multiply |ẑ(u)|2 by ρg ≈ 104 kg m−2 s−2 and you’ll
have Joules per square meter. (The surface-tension term contributes a negligible part of the total
variance compared to the larger gravity waves. Surface tension is significant only for a very young
sea surface when only capillary waves are present.) Some writers (e.g. Holthuijsen, 2007) take care
to distinguish between the variance spectrum |ẑ(u)|2 , which describes the statistical properties of
the waves, and the energy spectrum ρg|ẑ(u)|2 , which describes a physical property of the waves.
As forewarned in §1.4, the transition from discrete to continuous spectra requires care. (To
complicate matters, the literature is rich with different terms and symbols used for the same
quantities, and the same symbol used for different quantities.) I’ll use S(u) to denote the discrete
variance spectrum |ẑ(u)|2 , u = 0, 1, ..., N − 1. S(u) shows how much variance is contained in a
finite frequency interval ∆ν centered at frequency νu . The continuous version of this discrete point
function is the variance spectral density S(ν), which shows how much variance there is per unit
frequency interval at frequency ν. The conceptual relation between the discrete and continuous
versions of S can be written
S(u)
S(ν) = lim . (2.3)
∆ν→0 ∆ν
In this limit process, we are thinking of taking a finer and finer resolution of the sea surface
elevation z(xr ), r = 0, 1, ...N − 1 so that N → ∞ for a fixed physical domain of length L. The
discrete z(xr ) then approaches a continuous function z(x). The discrete Fourier transform of z(xr )
then gives more and more frequencies νu with smaller and smaller frequency intervals ∆ν. An
important consequence of Eq. (2.3) that the continuous variance spectral density S(ν) has units of
m2 /frequency, which in the present example is m2 /(1/m), whereas the discrete variance spectrum
S(u) has units of m2 . S(ν) is called the omnidirectional variance spectrum. It is derived in a
slightly different way in Appendix §A.1, which discusses continuous wave spectra in both one and
two dimensions. Those results will be cited as needed. If you are new to wave spectra, take a look
at that discussion before proceeding.
Plots of the Fourier amplitude squared are also often called “power spectra” because in many
fields it is power rather than energy that is of interest. An example is the power being dissipated
in an electrical circuit. In that case, the power is proportional to the time-dependent electrical
current squared; current replaces surface elevation in our equations and the relevant frequency is
the temporal frequency. It is now common that the Fourier amplitude squared of anything is called
a power spectrum, regardless of what physical variable is under consideration.
As we know from Eq. (1.11), the 1-D surface elevation z(x) is related to the Fourier amplitude
ẑ(k) by Z ∞
1
z(x) = ẑ(k)e+ikx dk .
2π −∞
Differentiating this equation with respect to x gives the 1-D slope of the sea surface as
Z ∞
dz(x) 1
σ(x) ≡ = ẑ(k)ike+ikx dk .
dx 2π −∞

14
This leads us to identify ikẑ(k) as the Fourier amplitude corresponding to the sea surface slope.
This gives us two ways to study the slope statistics of random sea surfaces, given the Fourier
amplitude ẑ(k) (which we will learn to create from variance spectra in the following sections). The
first way is to take the inverse Fourier transform of ẑ(k) to obtain z(x), and then to differentiate
z(x) to get the slope. The second way is to take the inverse transform of ikẑ(k) to get the slope
directly, without ever creating the surface z(x) itself.
It is shown in §A.1 that the variance of the sea surface slope, σ 2 , usually called the mean square
slope or mss, is given by Z ∞
σ 2 = mss = k 2 S(k) dk ,
0
where S(k) is the omnidirectional variance spectrum introduced above. As we saw, S(k) has units
of m2 /(rad/m), so the integral above gives the mss in units of rad2 . This is nondimensional, but
the label of “rad” reminds us that we can think of the slope as an angle (from the horizontal)
measured in radians. The product k 2 S(k) is called the slope spectrum.
The extensions of these spectra to two dimensions are made in §A.1, which summarizes a
number of useful relations for one- and two-dimensional variance and slope spectra as needed for
this tutorial.

2.3 Surfaces to Variance Spectra

We begin with a contrived example. We will comment in detail on this simple example in order
to reveal the mathematical subtleties and physical characteristics of Fourier transforms and wave
variance spectra derived from surface elevations. A thorough understanding of this example takes
us much of the way to understanding the case of a real ocean surface.
An ad hoc wave profile is constructed using the formula
N/2
X
z(xr ) = z(r) = A(j) cos[2πjxr νf + φ(j)] r = 0, ..., N − 1 . (2.4)
j=0

The xr locations are given by r∆x = rL/N , where L is the length of the sea surface region being
sampled and N is the number of samples. νf = 1/L is the fundamental spatial frequency, that is,
the spatial frequency or wavenumber of the wave with a wavelength of L. The amplitude of the
wave at the j th frequency, j = 1, 2, ..., N/2, is chosen to be

A(j) = 0.1 exp(−3j/N ) ,

and A(0) = 0. The phase of the j th wave component is randomly distributed over [0, 2π) using

φ(j) = 2πU

where U is a uniform [0, 1) random number. A different random number is drawn for each j value.
The upper left panel of Fig. 2.1 shows the surface generated in this manner for L = 10 m,
N = 16, and a particular set of random phases. Note that N is a power of 2 as will be needed
for the FFT. The thin colored lines show the N/2 + 1 = 9 waves for each of the frequencies. The
blue line is the wave for the fundamental frequency νf = 1/L = 0.1 m−1 ; the thin black line is the

15
two-point wave at the Nyquist frequency νNy = 1/(2∆x) = 0.8 m−1 ; the purple line is the constant
j = 0 wave, which is set to z = 0 for the mean sea surface. The black dots connected by the thick
black line show the sum of the individual waves. These points represent a discrete sampling of the
continuous sea surface elevation.
In this example, the sampled region of the sea surface is L = 10 m in length, but the N = 16
samples do not include the point at x = 10 m. This is because the surface elevation at x = L is
always the same as at x = 0 when using Fourier techniques. Resolving the surface as a sum of
sinusoids that are harmonics of the fundamental frequency νf = 1/L gives sinusoids that always
return to their initial value after distance L. Real sea surfaces are of course not periodic, but we
do not know the true value at L because it was not measured by the present sampling scheme.
(Likewise, we do not know the true surface elevations at points in between the sampled locations.)
When we use Fourier techniques to generate random surface realization, we are always generating
a sea surface that is a periodic tiling; the tile dimension is L. This periodicity is useful if we want
to generate a visual rendering of a large region of sea surface from a smaller computed region; the
edges of the small tiles will match perfectly and the larger surface will often look reasonable, if you
don’t look too closely.
We now take the sequence of the N = 16 real wave elevations z(r) seen in Fig. 2.1 and feed them
into an FFT routine. We soon get back 16 complex numbers, the ẑ(νu ) = ẑ(u) Fourier amplitudes,
at a set of 16 corresponding frequencies νu . The upper right panel of Fig. 2.1 plots the real part of
the ẑ(u) complex numbers, and the lower-left panel plots the imaginary part.
We note first that the FFT routine returned both negative and positive spatial frequencies: ν =
(− N2 + 1)∆ν = −0.7 m−1 , ..., −0.1, 0, 0.1, ..., N2 ∆ν = 0.8 m−1 , for a total of N = 16 discrete spatial
frequencies. Note that the frequency spacing ∆ν equals the fundamental frequency νf = 1/L. The
spatial frequency νu tells how many waves fit into 1 m; this is clearly a positive number with an
easily understood physical meaning. Therefore the negative frequencies returned by the FFT may
seem somewhat mysterious. However, these negative frequencies are simply the mathematical price
we pay for the convenience of using complex numbers. Consider, for example, the representation
of the cosine as a sum of complex exponentials for the uth frequency:
1 i2πνu xr
+ e−i2πνu xr

cos(2πνu xr ) = e
2
1
= ei2π(+νu )xr + ei2π(−νu )xr .
2
We can interpret the complex representation of the real cosine as having one term with a positive
frequency +νu and one term with a negative frequency −νu . A similar equation holds for the
complex representation of sin(2πνu xr ). The Fourier transform of a real function always contains
both negative and positive frequencies, which arise from the complex exponential in the definition
of the transform.
Note next that the real parts of the complex amplitudes ẑ(u) are even functions of frequency:
Re{ẑ(−νu )} = Re{ẑ(+νu )}. The imaginary parts are odd functions of frequency: Im{ẑ(−νu )} =
−Im{ẑ(+νu )}. In Fig. 2.1 the positive and negative frequencies are connected by a red arrow
for one particular frequency pair, νu = ±0.6 m−1 . If a complex function c(ν) = a(ν) + ib(ν) has
an even real part a(ν) and an odd imaginary part b(ν), then c∗ (−ν) = c(ν), i.e. the function is
Hermitian. Thus the plots verify that the amplitudes are Hermitian, as is always the case for the

16
Figure 2.1: A one-dimensional random sea surface and its Fourier transform and variance spectrum.
The black dots in the upper-left panel show the points z(u) of the sampled sea surface. The light
colored lines show the sinusoidal components used to create the surface. The upper-right panel
shows the real part of the complex amplitude ẑ(νu ) and the lower-left panel shows the imaginary
part. The lower right panel shows the two-sided and one-sided discrete variance spectra. Compare
with Figs. 2.2 and 2.3.

17
Fourier transform of a real function. The Hermitian character of the complex amplitudes means
that these N = 16 complex numbers contain only 16 independent real and imaginary numbers,
not 32 as would be the case for 16 arbitrary complex numbers. In general the FFT of N real
numbers (e.g., N spatial samples of a sea surface) gives back N independent numbers, so that the
“information content” of the physical and Fourier representations is the same.
The positive frequency at N2 ∆ν = 0.8 m−1 is the Nyquist frequency. There is, however, no value
for the negative of the Nyquist frequency. Note also in the lower left panel that the imaginary part
of the amplitude is identically zero at the Nyquist frequency. We will explain these values below.
The lower-right panel of the figure shows the values of |ẑ(νu )|2 . The values at the negative to
positive frequencies are connected by the black dotted line. These points constitute the two-sided
discrete variance spectrum, S2S (νu ) = |ẑ(νu )|2 . “Two-sided”, denoted by the subscript, refers to
spectra showing both the negative and positive frequencies. The variance at zero frequency is the
variance contained in the constant mean sea level. This value is zero because we have set the mean
sea level to zero.
We are usually concerned only with the variance at a given magnitude of the spatial frequency,
and not with whether the frequency is negative or positive. Nor is there any reason to plot the
point at zero frequency, which is usually zero by choice of zero for the mean sea level. It is therefore
customary to define the one-sided variance spectrum

S1S (νu ) = S2S (−νu ) + S2S (νu ) ,

for u = 1, 2, ..., N2 − 1, and S1S (νNy ) = S2S (νNy ). Then only the positive frequencies are plotted.
The points connected by the solid black line in the lower-right panel of Fig. 2.1 comprise the
one-sided variance spectrum. In the present simulation, the two-sided spectrum is symmetric for
positive and negative frequencies (except for the Nyquist frequency, which does not have a negative
counterpart and is always a special case), and the one-sided function is simply twice the value of the
two-sided function for the positive frequencies, except for the Nyquist frequency. When you read a
paper and it refers to or plots “the variance (or energy) spectrum” without further comment, it is
always the one-sided spectrum. We will have to account for the magnitude difference in one- and
two-sided spectra in the next section.
There is an important detail to note in the computation and plotting of S(u), as in the lower-
right panel of Fig. 2.1. The values of S(u) were obtained by the discrete Fourier transform of
Eq. (1.14), and S(u) gives the variance contained in a finite frequency interval ∆ν = 1/L at the
discrete frequency νu . ∆ν equals the fundamental frequency and is the frequency interval used in
the calculations and the plot. As noted in the discussion of the discrete transform, S(u) is a point
function. As already seen in Eq. (2.3), if we wish to convert the discrete S(u) to an estimate of the
continuous variance spectral density S(ν), we must divide by the discrete function by the frequency
interval: S(ν) = S(u)/∆ν. The units of S(ν) are then m2 /(1/m), as expected for a spectral density
function of spatial frequency. I will strive to distinguish between a discrete variance point function
and a continuous variance spectral density.
Let us return to Eq. (2.4) and set all of the phases φ(j) to zero. We are then adding together
cosines to create the surface wave profile, which is seen in the upper left panel of Fig. 2.2. The
FFT of this profile gives the real part of ẑ(νu ) as positive numbers except for the 0 frequency, and
the imaginary part is identically zero for all frequencies.

18
If we set all of the phases φ(j) to π/2, we are then adding together sines to create the surface
wave profile, which is seen in the upper left panel of Fig. 2.3. The FFT of this profile gives the
real part of ẑ(νu ) identically zero and the imaginary part is nonzero except for the 0 and Nyquist
frequencies.
These two figures show that the real part of the complex amplitude ẑ(νu ) tells us how much of
z(xr ) is composed of cosine waves, and the imaginary part shows how much of z(xr ) is composed
of sine waves. This explains why the imaginary part of the amplitude is zero at the Nyquist
frequency. The two-point wave at the Nyquist frequency is inherently a cosine wave because, as
noted previously, a two-point sine wave is sampled only at its zero values. The general case of a
wave component with a phase that is neither 0 (nor a multiple of 2π) nor π/2 (nor an odd integer
multiple of π/2) can be written as a sum of cosine and sine waves, as in Eq. (1.3). In that case,
both the real and imaginary parts of the amplitudes are nonzero (except for the special cases of
the 0 and Nyquist frequencies).
Note that in each of these three simulations, which differ only by the phases of the component
sinusoidal waves, the variance spectrum is exactly the same (except at the Nyquist frequency), as
seen in the lower right panels of Figs. 2.1-2.3. That is to say, the variance contained in a wave
does not depend on the reference coordinate system used to describe it, even though the Fourier
amplitudes ẑ(u) do depend on the coordinate system. The variance depends only on the amplitude
of the wave. The variance at the Nyquist frequency is largest when cosine waves are added and is
zero when sine waves are added. In the first case, we have the maximum possible amplitude of the
two-point cosine wave, and in the latter case there is no two-point wave.

2.4 Frequency Ordering

It is necessary to mention a peculiarity of all FFT subroutines. The plots above show the frequencies
ordered from negative to positive, as is convenient for plotting of two-sided spectra. Let ∆f
represent the fundamental frequency. If wavenumber ν is the frequency variable, then ∆f = 1/L.
If angular spatial frequency k is the frequency variable, then ∆f = 2π/L; for temporal angular
frequency ω, ∆f = 2π/T . In any case, the frequencies are equally spaced at intervals of ∆f and
are in the order
N N N N
− + 1, − + 2, ..., −1, 0, 1, ..., − 1, ∆f . (2.5)
2 2 2 2
I’ll call this the “math frequency order.”
However, FFT routines generally return their frequencies in the order of 0 first, then the positive
frequencies, then the negative frequencies in reverse order:

N N N N
0, 1, ..., − 1, , − + 1, − + 2, ..., −1, ∆f . (2.6)
2 2 2 2
I’ll call this the “FFT frequency order.” Given the Hermitian symmetry of the amplitudes about the
0 frequency, the FFT order is convenient for creating and plotting a one-sided variance spectrum,
in which case the last half of the amplitude array is ignored.
Either frequency order can be obtained from the other by a repeated circular shift, which moves
an array element off of one end of an array and copies it to the other end of the array, shifting all
elements right or left by one position in the process. The detail to watch is that different computer

19
Figure 2.2: A one-dimensional surface composed of cosine waves. Compare with Figs. 2.1 and 2.3.

20
Figure 2.3: A one-dimensional surface composed of sine waves. Compare with Figs. 2.1 and 2.2.

21
languages define a circular shift in different ways. For example, the IDL routine SHIFT (and the
Matlab routine CIRCSHIFT) moves the array elements to the right for a “positive” shift (a negative
shift moves elements to the left), whereas the Fortran 95 CSHIFT routine moves the array elements
to the left for a positive shift (a negative shift moves elements to the right). Thus

In IDL:

math order = SHIFT(FFT order, N/2-1)

FFT order = SHIFT(math order, -(N/2-1))

In Fortran95:

math order = CSHIFT(FFT order, -(N/2-1))

FFT order = CSHIFT(math order, N/2-1)

2.5 Variance Spectra to Surfaces

The previous section showed the major features of the Fourier analysis of a sea-surface elevation
record. We started with a sample of a sea surface and ended with the corresponding discrete variance
spectrum (or estimate of the variance spectral density after division by the frequency interval). In
this section we learn how to go in the reverse direction: start with a variance spectrum and generate
a random realization of the corresponding sea surface.
We will generate examples using two widely used variance spectra. The Pierson and Moskowitz
(1964) spectrum describes fully developed gravity waves. Fully developed waves are an idealization
of the sea state after the wind has been blowing for an infinitely long time over an infinite fetch. The
Elfouhaily et al. (1997) spectrum describes the full range of gravity to capillary waves, including
a parameter for the “age” of the waves, from young (the wind has just started blowing) to fully
developed seas. The mathematical formulations of these spectra are given for reference in the
Appendix.

2.5.1 Theory
The essence of the process is as follows:

1. Choose the domain size. To create a sea surface over a spatial domain at a given time,
we pick the length L of the region [0, L]. To generate a time series at a given point, we pick
the length of the time series.

2. Choose the number of points for the FFT. This number N is the number of frequencies
at which we will sample the variance spectrum, and equals the number of samples of the sea
surface that will be generated. N must be a power of 2 for the inverse FFT.

3. Choose the frequency variable. To generate a sea surface over a spatial domain at a
given time, we can use either wavenumber ν or angular spatial frequency k. To generate a
time series at a given point, we can use either f or ω.

22
4. Choose a variance spectrum The variance spectrum must be expressed in terms of the
chosen frequency variable. §A.2 gives the Pierson-Moskowitz one-sided spectrum as a function
of ν, k, or ω and shows how to change variables in a density function.

5. Choose the wind speed. Pick a wind speed, and perhaps other physical parameters such
as the age of the waves to be generated if required by the chosen variance spectrum.

6. Create random Hermitian amplitudes. This is the tricky part. We must create an
array of randomized discrete Hermitian Fourier amplitudes ẑ(u), starting with the chosen
continuous variance spectrum.

7. Take the inverse FFT of the random amplitudes. The inverse FFT converts the Fourier
amplitudes to the physical wave heights.

8. Extract the sea surface heights. The inverse FFT of the complex amplitudes returns a
complex array. The real part of this array is the random realization of the sea surface heights,
and the imaginary part is zero.

9. Check your results. This is extremely important during code development. For example,
take the FFT of the generated surface heights to see if you get back to the Fourier amplitudes
and variance spectrum you started with.

We now proceed through these steps and discuss them in detail for a specific example.
Steps 1 and 2: Let us generate a sea surface over the region from x = 0 to x = L = 100 m.
The longest wavelength that can be resolved is then 100 m. We use N = 1024, which gives a spatial
grid resolution of ∆x = L/N = 0.0977 m. This means that the shortest wavelength that can be
resolved, the two-point wave or Nyquist wavelength, is 2∆x = 0.1954 m.
Step 3: We used wavenumber ν in the previous section because of its easy interpretation. Now
let’s use angular spatial frequency k, which is more commonly used. The fundamental frequency
is then kf = 2π/L = 0.0628 rad m−1 . The highest frequency sampled, the Nyquist frequency, is
kNy = (N/2)kf = 32.15 rad m−1 . Note that (2π)/kNy = 0.1954 m which, as noted above, is the
wavelength of the two-point wave.
Step 4: We for this example, we use the Pierson-Moskowitz variance spectrum in terms of
angular spatial frequency k. This is given by Eq. (A.11). Note that this is a one-sided spectrum,
which is defined for positive k values.
Step 5: The wind speed at 10 m elevation is U10 = 5 m s−1 . The wind speed is the only input
to the Pierson-Moskowitz spectrum.
Step 6: We now discuss in detail how to generate a set of random discrete Fourier amplitudes
that are physically consistent with the chosen variance spectrum. These amplitudes must be defined
for both positive and negative frequencies, and the amplitudes must be Hermitian. We first define
r
1 S1S (k = ku )
ẑo (ku ) ≡ √ [ρ(ku ) + iσ(ku )] ∆k . (2.7)
2 2

Here S1S (k = ku ) denotes the continuous spectral density S1S (k) evaluated at k = ku . ∆k = kf
is the spatial frequency sampling interval. ẑo (ku ) must be defined for both positive and negative
discrete frequencies in order to create the Hermitian amplitudes for use in the inverse FFT. Recalling

23
the comments of the previous section, we convert the one-sided continuous spectrum S1S (k) into a
two-sided discrete variance function by
1. dividing its magnitude by 2, assuming that S2S (−k) = S2S (k);

2. multiplying the continuous spectral density by the fundamental frequency interval ∆k, which
gives the variance contained in a finite frequency interval at each frequency ku .
To emphasize the discrete vs continuous functions, and for brevity of notation, let us write the
frequency index u for the frequency ku . Then Eq. (2.7) becomes
1 p
ẑo (u) = √ [ρ(u) + iσ(u)] S2S (u) , (2.8)
2
where S2S (u) denotes the two-sided discrete variance spectrum at frequency ku . ẑo (u) can now
be evaluated for both positive and negative ku . The 0 and Nyquist frequencies are always special
cases: set S2S (0) = 0 and S2S (kNy ) = S1S (kNy ). ρ(ku ) ≡ ρ(u) and σ(ku ) ≡ σ(u) are independent
random numbers drawn from a normal distribution with zero mean and unit variance, denoted
ρ, σ ∼ N (0, 1). A different pair is drawn for each u value.
ẑo (u) is a random variable. Let h...i denote the expectation of the enclosed variable. The
expected value of |ẑo (u)|2 , hẑo (u)ẑo∗ (u)i, gives back whatever variance function is used for S2S (u):
D 1
1
E
∗
p p
hẑo (u)ẑo (u)i = √ [ρ(u) + iσ(u)] S2S (u) √ [ρ(u) − iσ(u)] S2S (u)
2 2
S2S (u)
hρ2 i + hσ 2 i = S2S (u)

=
2
because hρ2 i = hσ 2 i = 1 for N (0, 1) random variables. Thus ẑo (u) is consistent with the chosen
variance spectrum. However, ẑo (u) is not Hermitian, so the inverse FFT would not give a real sea
surface.
Next define ẑ(u) as
1
ẑ(u) ≡ √ [ẑo (u) + ẑo∗ (−u)] . (2.9)
2
This function is clearly Hermitian, so the inverse FFT applied to ẑ(u) will give a real-valued
z(xr ) ≡ z(r). Moreover, this ẑ(u) is consistent with the variance spectrum:

h|ẑ(u)|2 i = hẑ(u)ẑ ∗ (u)i

D 1 1 p 1 p

= √ √ [ρ(u) + iσ(u)] S2S (u) + √ [ρ(−u) − iσ(−u)] S2S (−u) ×
2 2 2
E
1 1 p 1 p
√ √ [ρ(u) − iσ(u)] S2S (u) + √ [ρ(−u) + iσ(−u)] S2S (−u)
2 2 2
1D 2
ρ (u) − iρ(u)σ(u) + iσ(u)ρ(u) + σ 2 (u) S2S (u)+

=
4 p p
[ρ(u)ρ(−u) + iρ(u)σ(−u) + iσ(u)ρ(−u) − σ(u)σ(−u)] S2S (u) S2S (−u)+
p p
[ρ(−u)ρ(u) − iρ(−u)σ(u) − iσ(−u)ρ(u) − σ(−u)σ(u)] S2S (−u) S2S (u)+
2 E
ρ (−u) + iρ(−u)σ(−u) − iσ(−u)ρ(−u) + σ 2 (−u) S2S (−u)

1
= [S2S (u) + S2S (−u)] = S2S (u) .
2

24
Here we have noted that hρ(u)ρ(−u)i = 0, etc., because the random variables are uncorrelated for
different u values.
Equations (2.7) and (2.9) are the key to generating random sea surfaces from variance spectra.
ẑ(u) defined by these equations contains random noise, which leads to a sea surface with random
amplitudes and phases for the component waves of different frequencies. Any one of these surfaces
has a variance spectrum that looks like the chosen spectrum plus random noise. However, on
average over many realizations, the the noise in these spectra will average out, leaving the variance
spectrum. Figure 2.4 below illustrates these important ideas, but first we must complete the surface
generation.
Step 7: Compute the inverse FFT of the ẑ(u) of Eq. (2.9). The result is a complex function
Z(x):
Z(xr ) ≡ Z(r) = F F T −1 {ẑ(u)} .
A crucial warning to this step is that the u = 0, ..., N − 1 elements of the ẑ(u) array must be
in the FFT frequency order given by Eq. (2.6). Z(r) is returned with xr values in the order from
x0 = 0 to xN −1 = (N − 1)∆x.
Step 8: Extract the surface. The inverse FFT returns a complex array Z(xr ) whose real part
is the surface elevations z(xr ) and whose imaginary part is 0. The surface elevations are extracted
as the real part of Z(xr ):
z(xr ) = Re{Z(xr )} .
√
Step 9: Check the results. There are many places along the way to lose a 2 or to mess up
array indexing. At the minimum, it is worthwhile to check that the mean of the generated surface
is zero, and that the imaginary part of Z(xr ) = 0 (to within a small amount of numerical roundoff
error).
When developing computer code, or when first learning this material, it is also a good idea to
take the forward FFT of Z(xr ) to make sure that the input Fourier amplitudes ẑ(u) are recovered,
and that the variance spectrum corresponding to z(xr ) is consistent with the one chosen in Step 4.
Indeed, this check was what led to these notes.
Equations (2.7) and (2.9) are, with minor changes in notation, Eqs. (42) and (43), respectively,
of Tessendorf (2004, §4.3). However, Tessendorf’s version of Eq. (2.7) appears to use a one-sided
variance spectrum (his example used the one-sided Phillips spectrum of his Eq. (40)) without the
division by 2 seen in Eq. (2.7), which is needed to convert the one-sided spectrum to a two-sided
spectrum. Nor does he show the ∆k factor needed to convert a continuous spectral density to a
√
discrete function. His version of Eq. (2.9) does not contain the overall factor of 1/ 2 seen above.
These missing factors mean that in a round-trip calculation

variance spectrum → FFT−1 → sea surface → FFT → variance spectrum ,

you do not get back to the original variance spectrum. In other words, the Tessendorf equations
do not conserve wave variance (i.e., wave energy). Even if he included the ∆k factor in his actual
√
computations, the missing factors of 1/ 2 in his versions of our Eqs. (2.7) and (2.9) give an overall
factor one-half on the amplitudes, which corresponds to a factor of four error in the variance. That
is, waves generated using the Tessendorf equations have amplitudes that are too large.
Tessendorf (2004) discusses much more than just Fourier transform techniques, and his notes
have been very influential in the computer graphics industry. In 2008 he deservedly received

25
an Academy Award for Technical Achievement for showing the movie industry how to generate
and render sea surfaces, as well as for his many other pioneering accomplishments in efficiently
computing and rendering fluid motions into visually appealing images. (The first movie to use his
techniques was Waterworld, followed by dozens of others including Titanic.) When I checked with
him (personal communication) about the missing numerical factors, he readily acknowledged that
Eqs. (2.7) and (2.9) are the correct ones, but noted that Hollywood doesn’t care about conservation
of energy. I suppose that should be no surprise, since movies seem to have no problem with rockets
going faster than the speed of light, sound propagating through the vacuum of outer space, or time
travel that violates causality. Tessendorf’s equations are widely cited (especially in the computer
graphics literature), always without comment about the missing scale factors. I am therefore
forced to conclude that many people have created sea surfaces without checking their calculations
for internal consistency and physical correctness. That may be acceptable in a fantasy world, but
such laxness is not permissible in science. I wrote this tutorial in part to make clear exactly how
surfaces should be generated so that they are physically correct and suitable for scientific work.

2.5.2 Example: A Roundtrip Calculation

Figure 2.4 shows an example of 1-D surface waves generated using the Pierson-Moskowitz spectrum
for a wind speed of 5 m/s, and the recovery of the variance spectrum from the generated surface.
The blue curve in the upper-left panel shows the Pierson-Moskowitz spectrum as defined by Eq.
(A.11). The red dots show the frequencies at which the continuous spectrum is sampled. Those
dots blur together at the higher frequencies because of the log scale, but the ku points are equally
spaced at intervals of the fundamental frequency ∆k = kf = 2π/L = 0.0628 rad/m. The last
sampled frequency is kNy = 32.17 rad/m. The bottom panel shows the sea surface elevations z(xr )
generated for a particular sequence of random numbers ρ(u), σ(u).
The red line in the upper-right panel shows the function

|ẑ(u)|2 1
P(u) ≡ = |F F T {z(r)}|2 . (2.10)
∆k ∆k
P(u) is the discrete variance function for this particular z(xr ) surface. Schuster (1898) called
P(u) a periodogram. The periodogram P(u) contains random noise because z(xr ) is a random
realization of the sea surface, which was generated by applying random noise to the the variance
spectrum. This particular z(xr ) is analogous to a particular measurement of the sea surface. Had
we drawn a difference sequence of random numbers for use in Eq. (2.8), we would have generated
a different sea surface, and a different P(u). However, we can expect that if we average together
many different P(u), corresponding to many different sets of z(xr ), the noise would average out
and we would be left with a curve close to the variance spectrum we started with, which is shown in
blue. Numerical experimentation shows that averaging 100 P(u) generated from 100 independent
sea surface realizations gives an average P(u) that is almost indistinguishable from the blue curve
at the scale of this plot. Thus P(u) ≡ P(ku ) is an approximation of the variance density spectrum
.
S(k), denoted P(ku ) = S(k). This averaging processes leads to the topic of spectrum estimation,
which considers such problems as how many sets of measurements of z(xr ) are needed to estimate
the variance spectrum to within certain error bounds. Fortunately, we need not pursue that here.
At the minimum, you should always check to see that Parseval’s relation, Eq. (1.20), is satisfied.

26
Figure 2.4: Example of a 1-D random sea surface generated from the Pierson-Moskowitz spectrum
for U10 = 5 m s−1 (the wind speed at 10 m above mean sea level). The longest resolvable wave has
wavelength L = 100 m. For N = 1024, the two-point wave has wavelength 2L/N = 0.195 m. This
is already less than the smallest wavelength (highest frequency) for which the Pierson-Moskowitz
spectrum should be used.

27
For the simulation of Fig. 2.4, we have
N
X −1 N
X −1
2
|z(r)| = N |ẑ(u)|2 = 19.395 m2 .
r=0 u=0

There are sometimes other checks that can be made. For example, the Pierson-Moskowitz
spectrum is simple enough that it can be analytically integrated over all frequencies. This gives
Z ∞
U4
2
hz i = SPM (k)dk = 3.04 · 10−3 10 , (2.11)
0 g2

where we have recalled that the variance spectral density is related to the variance of the sea surface.
The variance of the generated zero-mean sea surface can be computed from
N
1 X
var(z) = z(xr )2 (2.12)
N
r=0

and compared with the analytical expectation. For the surface seen in Fig. 2.4, Eq. (2.12) gives
var{z} = 0.0189 m2 vs. the theoretical value of hz 2 i = 0.0197 m2 from Eq. (2.11). This agreement
to within a small amount of random noise indicates that all is probably well with the calculations.
Indeed, the average var(z)± one standard deviation for 100 independent simulations is 0.020±0.007,
which agrees well with the theoretical value.
The significant wave height H1/3 is by definition the height (trough-to-crest distance) of the
highest one-third of the waves. To a good approximation, this is related to the expectation of the
variance by p
H1/3 = 4 hz 2 i .
In the present example, this formula gives H1/3 = 0.55 m. The average significant wave height for
100 simulations is 0.56 ± 0.09 (average ± one standard deviation). If you spend enough time in a
sea kayak to develop intuition about wave heights as a function of wind speed, a half-meter height
for the largest waves is about right for a 5 m s−1 or 10 knot wind speed.
To summarize this chapter: we have learned how to start with a wave variance spectral density
function and generate random discrete Fourier amplitudes. The inverse FFT of those amplitudes
gives a random realization of a sea surface that is physically consistent with the chosen variance
spectrum. That is all that we can ask of the Fourier transform technique.

2.6 Numerical Resolution

The question naturally arises: “How large a region L and how many points N must be used in
generating sea surface realizations?” Not surprisingly, the general answer is, “It depends on your
application.” It is now computationally possible to create 2-D surfaces (as described in the next
chapter) with sufficiently large N values in the x and y directions, but run times are prohibitive.
Kay et al. (2011) created 2-D surfaces with 65536 × 65536 points in each dimension, which sampled
the energy spectrum for wavelengths from 200 m gravity waves to 3 mm capillary waves. However,
it took 6 hours to create just one surface on a 3 GHz computer. Many ray tracing applications use
tens of thousands of surface realizations to get good statistical averages. Thus, the specific answer

28
for brute-force scientific calculations is usually, “A bigger N than you can run on your computer.”
This section shows how to finesse certain calculations so that large-N calculations can be avoided.
For creating visually appealing sea surfaces for rendering into a movie scene, N = 210 =
1024, 211 = 2048, or perhaps at most 212 = 4096 is usually adequate. This is because the visual
impression of a sea surface is, to first order, determined by the height of the waves, which in turn is
governed by the largest gravity waves for the given wind speed. If you are rendering an image, you
don’t need to have a finer resolution for the sea surface than maps to an image pixel. For example,
if you are going to simulate what a CCD with 1024 × 1024 pixels would record, and you are looking
down onto the sea surface from above at an elevation where one CCD pixel sees a 0.1 m × 0.1 m
patch of the sea surface, then you need to resolve the sea surface with a grid spacing of at most 0.1
m. Modeling a 100 m × 100 m patch of sea surface with a spatial resolution of 1024 × 1024 gives a
grid spacing of 0.1 m, which would be adequate for the CCD image simulation.
Unfortunately, for accurate scientific calculations of sea surface reflectance, N may have to be
much larger. This is because the surface reflectance depends to first order on the surface wave
slope, and the slope is strongly affected by the smallest waves with spatial scales down to capillary
size of a few millimeters.
The k value of the peak of the Pierson-Moskowitz spectrum of Eq. (A.11) is easily obtained
from setting dSPM (k)/dk = 0 and solving for the value of k. This gives
0.5
2β g
kp = 2 ,
3 U19
which corresponds to a wavelength of Λp = 2π/kp . Table 2.6 shows the values of kp and Λp for a
few wind speeds.

U19 [m/s] kp [rad/m] Λp [m/s]

5 0.276 23
10 0.069 91
15 0.031 205
20 0.017 364

Table 2.1: Spatial frequencies and wavelengths corresponding to peak variance for the Pierson-
Moskowitz spectrum SPM (k).

If we use L ≈ 2Λp for the size of the spatial domain, we will resolve the large gravity waves, which
contain most of the wave variance. (Keep in mind that density functions do not have a unique
maximum: the location of the maximum of a density function depends on which frequency variable
is chosen. See the discussion of this on the Ocean Optics Web Book page on “A Common Miscon-
ception” at
http://www.oceanopticsbook.info/view/light_and_radiometry/level_2/a_common_misconception.
√
Differentiating Eq. (A.13) and using the dispersion relation ω = gk leads to Λp = 2π(1.25/β)0.25 g/U19
2 ,

which gives 73 m at U19 = 10. However, either version of Λp gives adequate guidance for our present
purpose.)
Now consider how many points it takes to almost fully resolve, or account for, the variances
of elevation and slope spectra when they are sampled in a FFT. For a specific example, we take

29
U10 = 10 m/s and L = 200 m. We use the omnidirectional, or 1-D, part of the Elfouhaily et al.
(1997) spectrum presented in §A.3. This S(k) is given by Eq. (A.16) and the associated equations.
This variance spectrum and the corresponding slope spectrum k 2 S(k) are plotted as the green
curves in Fig. A.3, and again by the blue curves in the next 3 figures.
As we have learned (c.g. §2.2 and A.1), the expected variance of the sea surface is given by
Z ∞
hz 2 i = S(k) dk , (2.13)
0

and the mean-square slope is given by

Z ∞
2
hσ i = mss = k 2 S(k) dk . (2.14)
0

Such integrals usually must be numerically evaluated with the limits of k = 0 and ∞ being replaced
by some small value k0 and some large but finite value k∞ . In the present case, I used k0 = 0.01 and
k∞ = 104 , with 106 evaluation points in between these limits. That gives an accurate evaluation
of Eqs. (2.13) and (2.14) for the U10 = 10 m s−1 spectra. The results are hz 2 i = 0.4296 m2 and
mss = 0.06011 rad2 .
Now suppose we sample the spectrum using N = 1024 points, which leads eventually to a
sea surface realization z(xr ) with 1024 points. For L = 200 m the fundamental frequency is
kf = 2π/L = 0.0314 rad/m. This point is shown by the left-most red dot with the black center in
Fig. 2.5. For this L and N , the Nyquist frequency is kNy = 2π/(2L/N ) = 16.085 rad/m, which is
shows by the right-most black and red dot. These two end points and the N − 2 red vertical lines
in between show the discrete points where the variance spectrum is sampled.
The elevation and slope variances that are accounted for in a sampling scheme with N points
are given by
Z kNy
z 2 (N ) = S(k) dk ,
kf

and Z kNy
mss(N ) = k 2 S(k) dk .
kf

In the present example, these integrals are z 2 (N

= 1024) = 0.4219 m2 and mss(N = 1024) =
0.02584 rad2 . Thus the finite range of the sampled frequencies includes the fraction

z 2 (N ) 0.4219
fE ≡ = = 0.982
hz 2 i 0.4296
of the total variance of the sea surface elevation. However, the corresponding fraction of the
sampled slope variance is just fS = 0.02584/0.06011 = 0.430. Thus N = 1024 is sampling 98% of
the elevation variance but only 43% of the slope variance. This sampling is acceptable if we are
interested only in creating a sea surface that looks realistic to the eye. The large gravity waves will
be correctly simulated, and that is what counts for visual appearance.
However, if our interest is ray tracing to simulate the optical reflectance and transmission
properties of the sea surface, then we must also sample almost all of the slope variance. This is
because, to first order, reflectance and transmission of light are governed by the slope of the sea
surface, not by its height. As will be seen in §3.4, if we under sample the slope varaince, then

30
Figure 2.5: Example sampling of elevation and slope spectra for N = 1024. The red dots and
lines show the regions of the elevation and slope spectra sampled using 1024 points. The elevation
spectrum is adequately sampled, but the slope spectrum is under sampled.

the generated surface will be too smooth at the smallest spatial scales, and the computed optical
properties will be incorrect. The brute-force approach to sampling more of the slope spectrum is
to increase N .
Figure 2.6 shows the sampling points when N = 216 = 65536. Now the sampled points extend
into the capillary-wave spatial scale: the smallest resolvable wave is 0.006 m. (Note that the
fundamental frequency is fixed by the choice of L, which also fixes the spacing of the samples,
∆k = kf .) The elevation variance integral over the sampled ranges is unchanged (we are still
missing a small bit of the variance from the very longest waves with frequencies below kf . However,
the slope integral is now mss(N = 65536) = 0.05909 rad2 . Thus we are now sampling a fraction
fS = 0.05909/0.06011 = 0.983 of the slope variance. The resolution of both the elevation and
slope spectra are now acceptable for almost any application. Although using a very large N is
computationally feasible in one dimension, such a large N is computationally impracticable in two
dimensions, when we would need a grid of size N × N . In the present example, this would would
require 655362 = 4.3 Gbytes of storage for just one real array, as well as a corresponding increase
in run time for the FFT. Another approach to resolving the slope variance is therefore needed.
As we have seen, the unsampled frequencies greater than the Nyquist frequency can account for
a large part of the slope variance. These frequencies represent the smallest gravity and capillary
waves. The amplitudes of these waves are small, so they contribute little to the total wave variance,
but their slopes can be large. One way to account for these “missing” waves is to alias their variance
into the waves of the highest frequencies. The highest frequency waves that are sampled will then
contain too much variance, i.e., they will have amplitudes that are too large for their wavelengths,
which will increase their slopes. One way to do this is as follows.

31
Figure 2.6: Example sampling of elevation and slope spectra for N = 65536. The red dots and
lines show the sampled regions of the elevation and slope spectra. Both the elevation and the slope
spectra are adequately sampled.

We seek an adjusted or re-scaled elevation spectrum

S̃(k) ≡ [1 + δ(k)]S(k) (2.15)

such that the integral of k 2 S̃(k) over the sampled region equals the integral of the true k 2 S(k) over
the full spectral range. Then sampling the re-scaled spectrum S̃(k) will lead to the same slope
variance as would be obtained from sampling the correct spectrum S(k) over the entire range of
frequencies.
We can choose L so that the low frequencies are well sampled starting at the fundamental
frequency kf = 2π/L. There is thus no need to modify the low-frequency part of the variance
spectrum, which if done, would adversely affect the total wave variance. We want to modify only
the higher frequencies of the sampled region, which contribute the most to the wave slope but little
to the wave elevation. A simple approach is to take δ(k) to be a linear function of k between the
spectrum peak kp and the highest sampled frequency kNy , and zero elsewhere:

0 if k ≤ kp
δ(k) ≡
k−k p
(2.16)
 δNy if k > kp
kNy −kp

δNy is a parameter that depends on the spectrum (i.e., the wind speed), the size L of the spatial
domain, and the number of samples N . We must determine the value of δNy so that this δ(k) “adds
back in” the unresolved slope variance. The added variance will be zero at the peak frequency kp
and largest at the Nyquist frequency. That is, the δ(k) function will make only a small change to the
variance spectrum at the low frequencies, and the change will be largest near the highest sampled
frequencies, which is consistent with the idea that the high frequency waves have the largest slopes.

32
(There is nothing physical about the functional form of Eq. (2.16); it is simply a single-parameter
ad hoc function that works well in practice. Nonlinear functions could be used, e.g. to concentrate
more of the variance into the frequencies closest to the Nyquist frequency. However, those functions
could have more than one free parameter to determine and, in any case, the end result is the same:
a sea surface that reproduces the height and slope variances of the fully sampled spectrum.)
To determine δNy we thus have
Z ∞ Z k∞ Z kNy Z k∞
2 2 2
mss = k S(k) dk ≈ k S(k) dk = k S(k) dk + k 2 S(k) dk
0 k0 k0 kNy
kNy kNy kNy
k − kp
Z Z Z
≡ k 2 S̃(k) dk = k 2 S(k) dk + δNy k2 S(k) dk .
k0 k0 kp kNy − kp

The right-most terms of these equations give (after recalling that δ(k) = 0 for k ≤ kp )
R k∞
kNy k 2 S(k) dk
δNy = R k . (2.17)
Ny k−kp
kp k2 kNy −kp S(k) dk

Figure 2.7 shows example 1-D surface realizations created with N = 1024 samples and both
the true and scaled variance spectra. The upper two panels reproduce Fig. 2.5, except that the
sampled points of the re-scaled spectra are shown in green. It is clear that the δ(k) function has
added progressively more variance to the higher frequencies. The re-scaled variance spectrum does
of course contain somewhat more variance over the sampled region than does the true spectrum.
As the green inset fE number shows, this re-scaling has increased the fraction of sampled/true
variance from 0.982 to 1.020. However, the fS number in the upper-right panel shows that we are
now sampling 99.5% of the slope variance, as opposed to 43% for the true spectrum. Having a bit
too much total elevation variance seems to be a good tradeoff for being able to model the optically
important slope statistics.
The row 2 panel shows random realizations of the 1-D surfaces generated from these two spectra
(with the same set of random numbers). The surface elevations z(xr ) from the true spectrum (red
line) and the re-scaled spectrum (green line) are almost indistinguishable at the scale of this figure.
These significant wave heights for these two surface realizations compare well with the significant
wave height of H1/3 = 2.48 m for a Pierson-Moskowitz spectrum. The row 3 panel shows the surface
slopes computed from finite differences of the discrete surface heights. The Cox-Munk along-wind
mean square slope is given by mssx = 0.0316U12.5 , which compares well with the value of 0.032
obtained with the re-scaled spectrum. However, the value obtained from the true spectrum is only
0.022, or 70% of the Cox-Munk value. The bottom two panels replot the 0-10 m sections of the z
and dz/dx figures for better display of the details.
We thus see that when using the δ(k) correction to a variance spectrum, we can generate a
surface that reproduces, to within a few percent of the theoretical values, both the surface elevation
and slope statistics that would be obtained from the underlying true variance spectrum if it were
sampled with enough points to fully resolve the elevation and slope variances.

33
Figure 2.7: Example 1-D surfaces and slopes created using true and adjusted variance spectra. See
the text for discussion.

34
CHAPTER 3

Two-Dimensional, Time-Independent Surfaces

We have now had a detailed look at the one-dimensional case. The extension to two dimensions is
mathematically straight forward.
Let z(x, t) = z(x, y, t) be the sea surface elevation in meters at point x = (x, y) at time t. The
spatial extent of the sea surface is 0 ≤ x < Lx and 0 ≤ y < Ly . This surface is sampled on a
rectangular grid of Nx by Ny points, where both Nx and Ny are powers of 2 for the FFT. The
spatial sampling points are then
Lx
x(r) = [0, 1, 2, ..., Nx − 1] = r∆x, r = 0, ..., Nx − 1
Nx
Ly
y(s) = [0, 1, 2, ..., Ny − 1] = s∆y, s = 0, ..., Ny − 1 .
Ny

This spatial sampling frequency gives Nx and Ny spatial frequencies kx and ky (in math frequency
order)
2π
kx (u) = [−(Nx /2 − 1), ..., −1, 0, 1, ..., Nx /2] = u∆kx , u = −(Nx /2 − 1), ..., Nx /2
Lx
2π
ky (v) = [−(Ny /2 − 1), ..., −1, 0, 1, ..., Ny /2] = v∆ky , v = −(Ny /2 − 1), ..., Ny /2
Ly
.

As before, the x-dimension 0 frequency is at array element u = Nx /2 − 1 and the Nyquist frequency
is at element u = Nx − 1. Thus the positive and negative pairs of kx values are related by
kx (u) = −kx (Nx − 2 − u), u = 0, ..., Nx − 2, with the Nyquist frequency always being a special case
because there is only a positive Nyquist frequency. Corresponding relations hold for the y direction.
Let k = (kx , ky ) denote a spatial frequency vector, where kx and ky are frequencies in the x
and y directions,
q respectively. For discrete values, we write kuv = (kx (u), ky (v)). The magnitude
of k is k = kx2 + ky2 . In our (x, y) coordinate system, let the wind blow in the +x direction.
The −x direction is then upwind, and the ±y directions are the cross-wind directions. With this
choice, kx > 0 indicates frequencies of waves propagating more or less downwind, and kx < 0 for
waves propagating against the wind. kx = 0 and ky 6= 0 indicates a wave propagating at exactly a

35
cross-wind ±y direction. The angle of wave propagation relative to the downwind direction for a
wave of frequency (kx (u), ky (v)) is given by

−1 ky (v)
ϕ(kuv ) = ϕ(u, v) = tan . (3.1)
kx (u)

We can thus write the 2-D surface as z(x), its Fourier amplitude as ẑ(k), and the associated variance
spectrum as Ψ(k). The discrete variance spectrum Ψ(kuv ) = Ψ(u, v) gives the variance of the wave
with wavelength 2π/kuv propagating in direction ϕ(kuv ) relative to the downwind direction. See
§A.1 for discussion of two-dimensional variance spectra.

3.1 Surfaces to Variance Spectra

Even though the mathematical transition from one to two dimensions causes no problems, it is
again educational to take a careful look at a couple of contrived examples.

3.1.1 Example: A Random Sea Surface

For the first example, a sea surface area of size Lx × Ly = 10 × 10 m is sampled using Nx = 16 and
Ny = 8 points in the x and y directions. Sea surface elevations z(x, y) were created by drawing a
N (0, 1) random number at each (x, y) value.
The upper-left panel of Fig. 3.1 shows a contour plot of this surface for a particular sequence
of random numbers. The spatial periodicity of the Fourier representation is used to extend the
contour plot to the full (Lx , Ly ) range of the tile, which gives a good visual appearance. Thus the
elevations at x = Lx are the same as those at x = 0, those at y = Ly equal those at y = 0, and the
elevation at (Lx , Ly ) duplicates the elevation at (0, 0). The Nx × Ny grid of sample points is shown
by the solid silver dots. The points at x = Lx and y = Ly obtained by periodicity are shown by
open silver circles at the right and top of the surface plot.
The Fourier amplitudes ẑ(k) = ẑ(kx , ky ) are obtained from the 2-D FFT of z(x) = z(x, y):

ẑ(k) = F F T2D {z(x)} .

The usual warning on frequency order applies here. The 2-D FFT gives back the amplitudes ẑ(k)
with the array elements corresponding the the FFT frequency order of Eq. (2.6). Before plotting,
the ẑ(k) array elements must be shifted into the math frequency order in both the kx and ky array
directions using the appropriate shift function for the computer language used in the plotting. For
the IDL routine used to generate Fig. 3.1, the 2-D circular shift is given by the IDL command

realzhatplot = SHIFT(REAL PART(zhat), Nx /2 − 1, Ny /2 − 1)

where zhat is the complex 2-D array returned by the FFT routine, and realzhatplot is the 2-D
array plotted in the upper right panel of Fig. 3.1.
The upper right panel of Fig. 3.1 plots the real part of ẑ(k), and the lower left panel plots the
imaginary part. The Nyquist frequency kxNy = 2π/(2∆x) = 5.03 rad/m lies along the right side
of the contour plot. The Nyquist frequency kyNy = 2π/(2∆y) = 2.51 rad/m lies along the top of
the contour plot. The white space at the left and bottom highlights that there are no negative

36
Figure 3.1: A two-dimensional random sea surface and its Fourier transform and two-sided variance
spectrum. The surface elevations in the upper left panel were randomly drawn from a N (0, 1)
distribution. The black arrows highlight a particular ±k frequency pair. The blue-blue and red-red
symmetry of the real part of the Fourier amplitudes, and the red-blue symmetry of the imaginary
part, shows the Hermitian nature of the amplitudes. The gray dots show the locations of the
discrete values that were contoured to create the figures.

37
Nyquist frequencies. In each amplitude plot a particular pair of ±k values is indicated by the
black arrows; the (0, 0) frequency is shown by a black dot. Note that in the plot of the real part,
Re{ẑ(−k)} = Re{ẑ(+k)}, whereas in the plot of the imaginary part, Im{ẑ(−k)} = −Im{ẑ(+k)}.
The contouring is rather low quality for so few points, but it is easy to see in the digital output
that when the kx is the Nyquist frequency (the points along the right column of points in the
plot), the symmetries are given by Re{ẑ(kxNy , −ky)} = Re{ẑ(kxNy , +ky)} and Im{ẑ(kxNy , −ky)} =
−Im{ẑ(kxNy , +ky)}. A corresponding relation holds for ±kx when ky = kyNy (the points along
the top row of the plot). The point at the upper right of the plot corresponds to both kx and
ky being at their respective Nyquist frequencies. As always, the array elements at the Nyquist
frequencies must be treated as special cases when writing computer programs. These symmetries
show that the 2-D amplitudes are Hermitian: ẑ ∗ (−k) = ẑ(k). The discrete 2-D variance spectrum
Ψ(k) = Ψ(kx , ky ) = ẑ(kx , ky )ẑ ∗ (kx , ky ) is contoured in the lower right panel of the figure. Note the
Ψ(−k) = Ψ(k) symmetry.
For this example, the standard check on Parseval’s relation (1.21) gives

x −1 N
NX y −1 x −1 N
NX y −1 x −1 N
NX y −1
X X X
2 2
|z(r, s)| = Nx Ny |ẑ(u, v)| = Nx Ny Ψ(u, v) = 130.90 m2 .
r=0 s=0 u=0 v=0 u=0 v=0

In all of these plots, it should be remembered that the discrete values are known only at the
locations of the silver dots. The contouring routine simply interpolates between these points to
create a visually appealing figure. The continuous color in the plots does not imply that the values
are continuous and known in between the discrete points.

3.1.2 Example: A Sea Surface of Crossing Sinusoids

For a second example, let us define a surface from a pair of crossing sinusoids as follows:

z(xr , ys ) = A1 cos(kx1 xr + ky1 ys + φ1 ) + A2 cos(kx2 xr + ky2 ys + φ2 ) , (3.2)

where as usual r = 0, ..., Nx − 1; s = 0, ..., Ny − 1; and where

A1 = 1.0 m is the amplitude of the first wave

Nx1 = 2 is the number of wave lengths in the x direction in [0, Lx ] for the first wave
kx1 = 2.0πNx1 /Lx = 1.257 rad/m2 is the kx frequency of the first wave
Ny1 = 1 is the number of wave lengths in the y direction in [0, Ly ] for the first wave
ky1 = 2.0πNy1 /Ly = 0.628 rad/m2 is the ky frequency of the first wave
φ1 = 0 is the phase of the first wave; 0 gives a cosine wave
A2 = 0.5 m is the amplitude of the second wave
Nx2 = 4 is the number of wave lengths in the x direction in [0, Lx ] for the second wave
kx2 = 2.0πNx2 /Lx = 2.513 rad/m2 is the kx frequency of the second wave
Ny2 = 3 is the number of wave lengths in the y direction in [0, Ly ] for the second wave
ky2 = −2.0πNy2 /Ly = −1.885 rad/m2 is the ky frequency of the first wave

38
φ2 = π/2 is the phase of the second wave; π/2 gives a sine wave
q
The wavelength of the first wave is Λ1 = 2π/ kx12 + k 2 = 4.47 m, and that of the second wave is
y1
Λ2 = 2.00 m. The direction of propagation of the first wave relative to the +x axis is

−1 ky1
ϕ1 = tan = 26.57 deg ,
kx1

and the direction of propagation of the second wave relative to the +x axis is

−1 ky2
ϕ2 = tan = −36.87 deg .
kx2

The upper left panel of Fig. 3.2 shows this surface elevation pattern when Eq. (3.2) is sampled
with Nx = Ny = 16. The dominant red-blue pattern shows the first wave oriented with the
direction of propagation along either the +k1 direction at ϕ1 = 26.57 deg or the −k1 direction at
ϕ1 = 26.57 + 180 = 206.57 deg. We cannot of course determine the actual direction +k1 or −k1 of
propagation from sea surface elevations at a single time. The dominant wave pattern is modulated
by the second wave, which has one half the amplitude of the first wave.
The choice above of φ1 = 0 in Eq. (3.2) makes the dominate wave a cosine in our coordinate
system, and φ2 = π/2 makes the second wave a sine. As we saw in the 1-D examples, variance
associated with cosine waves appears in the real part of the Fourier amplitudes, and the variance in
sine waves appears in the imaginary parts. We see this again here for the first (cosine) and second
(sine) waves of the surface. Since each wave pattern has only one frequency, there is only one pair
of points at ±k1 in the plot of the real part, and one pair at ±k2 in the plot of the imaginary part.
The amplitudes at all other frequencies are zero. The symmetries of these points again show the
Hermitian nature of the amplitudes.
The lower right panel shows the two-sided variance function Ψ2S (k) for this surface. Note that
the first wave has four times the variance of the second wave because the amplitude of the first
wave is twice that of the second wave. Ψ2S (−k) = Ψ2S (+k). Note also that you can look at a 2-D
variance spectrum and see how much energy is propagating in a given ±k direction.
For this simple example involving just four frequencies it is easy (from the digital output) to
hand check that the right-hand side of Parseval’s relation (1.21) is
XX
|ẑ(u, v)|2 = 16 · 16 (0.5)2 + (0.5)2 + (0.25)2 + (−0.25)2 = 160 m2 .

Nx Ny
u v

2
P P
This value agrees exactly with the corresponding sums of the surface elevations, r s [z(r, s)] ,
P P
and variance values, Nx Ny u v Ψ(u, v).

39
Figure 3.2: A two-dimensional sea surface composed of two crossing sinusoids, and the resulting
Fourier amplitudes and variance.

40
3.2 Variance Spectra to Surfaces
Now, again, we face the reverse process: start with a two-dimensional, one-sided variance density
spectrum and generate a random 2-D realization of a sea surface. The first task is to conjure up a
two-dimensional variance density spectrum, which requires some effort.

3.2.1 Theory
Let Ψ(k) = Ψ(kx , ky ) denote a 2-D elevation variance spectrum. This spectrum has units of
m2 /(rad/m)2 . By definition, the integral of Ψ(kx , ky ) over all frequencies gives the elevation vari-
ance: Z ∞Z ∞
var{z} = hz 2 i = Ψ(kx , ky ) dkx dky .
−∞ −∞

(See §A.1 for the precise definition of Ψ(kx , ky ) and for further details.) As in the 1-D case of §2.5, h i
indicates expectation or ensemble average over many measurements of the sea surface. The k values
coming out of ẑ = F F T {z} point both “downwind” (positive kx values) and “upwind” (negative kx
values), i.e. there are both positive and negative frequencies represented in the “two-sided” variance
spectrum, denoted as Ψ2S [kx (u), ky (v)], u = −(Nx /2 − 1), ..., Ny /2, v = −(Ny /2 − 1), ..., Ny /2 as
above. In general, Ψ2S (−k) 6= Ψ2S (k) because more energy propagates downwind than upwind
at a given frequency. kx (Nx /2) is the Nyquist frequency for waves propagating in the x direction;
ky (Ny /2) is the Nyquist frequency for waves propagating in the y direction. Ψ2S (0, 0) is the variance
at zero frequencies, i.e., the variance of the mean sea surface; this term is normally set to 0 so that
the mean sea surface is at height 0.
The 2-D equivalents of Eqs. (2.7) and (2.8) are
r
1 Ψ1S (k = kuv )
ẑo (kuv ) ≡ √ [ρ(kuv ) + iσ(kuv )] ∆kx ∆ky
2 2
1 p
= √ [ρ(kuv ) + iσ(kuv )] Ψ2S (kuv ) .) (3.3)
2
Here, as before, ρ(kuv ) and σ(kuv ) are independent N (0, 1) random variables, with a different
random variable drawn for each kuv value. As in the 1-D case, the expected value of ẑo (kuv )
gives back whatever variance spectrum is used for Ψ2S (kuv ), but is not Hermitian. As before, the
notation of Eq. (3.3) distinguishes between the value of the continuous spectral density function
evaluated at a discrete frequency value, Ψ1S (k = kuv ), and the discrete variance point function,
Ψ2S (kuv ).
We must define random Hermitian amplitudes for use in the inverse Fourier transform. Looking
ahead to the next chapter, which extends the results of this chapter to time-dependent surfaces,
define the time-dependent spectral amplitude
1
ẑ(kuv , t) ≡ √ ẑo (kuv )e−iωuv t + ẑo∗ (−kuv )e+iωuv t .

(3.4)
2
This function is clearly Hermitian, so the inverse FFT applied to ẑ(kuv , t) will give a real-valued
z(xrs , t). Recall from Eq. (1.1) that a function of the form cos(kx − ωt) gives a wave propagating in
the +x direction. The corresponding ẑo (kuv )e−iωt in this equation gives a wave propagating in the

41
+k direction, which is in the downwind half plane of all directions. In general the temporal angular
frequency ωuv is a function of the spatial frequency kuv . For example, for deep-water gravity waves,
2 = gk .
ωuv uv
For simplicity of notation, let us momentarily drop the rs and uv subscripts on the discrete
variables. The ẑ(k, t) of Eq. (3.4) is also consistent with the variance spectrum:

h|ẑ(k, t)|2 i = hẑ(k, t)ẑ ∗ (k, t)i

D 1 1 1

−iωt
p iωt
p
= √ √ [ρ(k) + iσ(k)] Ψ2S (k)e + √ [ρ(−k) − iσ(−k)] Ψ2S (−k)e ×
2 2 2
E
1 1 1
√ √ [ρ(k) − iσ(k)] Ψ2S (k)e−iωt + √ [ρ(−k) + iσ(−k)] Ψ2S (−k)eiωt
p p
2 2 2
1 D
r2 (k) − iρ(k)σ(k) + iσ(k)ρ(k) + s2 (k) Ψ2S (k)+

=
4 p p
[ρ(k)ρ(−k) + iρ(k)σ(−k) + iσ(k)ρ(−k) − σ(k)σ(−k)] Ψ2S (k) Ψ2S (−k)ei2ωt +
[ρ(−k)ρ(k) − iρ(−k)σ(k) − iσ(−k)ρ(k) − σ(−k)σ(k)] Ψ2S (−k) Ψ2S (k)e−i2ωt +
p p

2 E
r (−k) + iρ(−k)σ(−k) − iσ(−k)ρ(−k) + s2 (−k) Ψ2S (−k)

1
= [Ψ2S (k) + Ψ2S (−k)] .
2
Here we have noted that all terms like hρ(k)ρ(−k)i are zero because of the independence of the
random variables for different k values, as are terms like hρ(k)σ(k)i. The remaining term is the
total variance associated with waves propagating in the downwind and upwind directions at the
spatial frequency of magnitude k. It should be noted that this term is independent of time even
though the waves z(x, y, t) depend on time. This is because the total variance (or energy) of the
wave field is the same at all times, even though the exact shape of the sea surface varies with time.
If only a ‘snapshot’ of the sea surface at one time is available, it is not possible to resolve how
much of the total variance is associated with waves propagating in direction k compared to the
opposite direction −k. The forward FFT, ẑ(k, t) = F F T {z(x, y, t)}, and ẑ(k, t)ẑ ∗ (k, t) then gives
Ψ2S (−k) = Ψ2S (k), in which case the last equation reduces to

hẑ(k, t)ẑ ∗ (k, t)i = Ψ2S (k) .

In any case, the amplitudes defined by Eq. (3.4) are Hermitian, so that the real part of the
inverse Fourier transform Z(x, t) = F F T −1 {ẑ(k, t)} gives a real-valued sea surface z(x, t). That
sea surface is consistent with the variance spectrum Ψ2S (k) at every time t. Wave variance (energy)
is thus conserved in a “round-trip” calculation from variance spectrum to sea surface and back to
variance spectrum. In the time-dependent case, if Ψ2S (k) > Ψ2S (−k), then more variance will be
contained in waves propagating in the +k direction than in the −k direction. This is all that we
can ask from Fourier transform techniques.
Although Eqs. (3.3) and (3.4) are compact representations of the random spectral amplitudes,
the actual evaluation of these equations in a computer program warrants further examination. In
particular, the Nyquist frequencies are always special cases because there is only a positive Nyquist
frequency, kxNy = N2x L2πx at array element Nx − 1, with a corresponding equation for the Nyquist

42
frequency in the y direction. Writing e±iωt = cos[ω(k)t] ± i sin[ω(k)t] and expanding Eq. (3.4) gives
h p p i
2 ẑ(k, t) = ρ(k) Ψ2S (k) + ρ(−k) Ψ2S (−k) cos[ω(k)t]
h p p i
− σ(k) Ψ2S (k) + σ(−k) Ψ2S (−k) sin[ω(k)t]
(
h p p i
+ i ρ(k) Ψ2S (k) − ρ(−k) Ψ2S (−k) sin[ω(k)t]
)
h p p i
+ σ(k) Ψ2S (k) − σ(−k) Ψ2S (−k) cos[ω(k)t] . (3.5)

p
These terms are all Nx × Ny arrays, but note that terms like ρ(k) Ψ2S (k) represent element-by-
element multiplications, not matrix multiplications.
For a particular array element ẑ(kx (u), ky (v), t) = ẑ(u, v, t), and using the previously noted
indexing relation k(u) = −k(N − 2 − u), u = 0, ..., N − 2 for frequencies written in math order, we
can write

2 ẑ(u, v, t)
q
p
= ρ(u, v) Ψ2S (u, v) + ρ(Nx − 2− u, Ny − 2− v) Ψ2S (Nx − 2− u, Ny − 2− v) cos[ω(k(u, v))t]
q
p
− σ(u, v) Ψ2S (u, v) + σ(Nx − 2− u, Ny − 2− v) Ψ2S (Nx − 2− u, Ny − 2− v) sin[ω(k(u, v))t]
(
p q
+ i ρ(u, v) Ψ2S (u, v) − ρ(Nx − 2− u, Ny − 2− v) Ψ2S (Nx − 2− u, Ny − 2− v) sin[ω(k(u, v))t]
)
p q
+ σ(u, v) Ψ2S (u, v) − σ(Nx − 2− u, Ny − 2− v) Ψ2S (Nx − 2− u, Ny − 2− v) cos[ω(k(u, v))t] .

(3.6)

This equation allows for efficient computation within loops over array elements. In particular,
the code can loop over the non-positive kx (u) values, u = 0, ..., Nx /2, and over all ky (v) values,
v = 0, ..., Ny − 2 to evaluate the amplitudes for all non-Nyquist frequencies. The positive kx (u)
values are then obtained from the negative kx (u) values by Hermitian symmetry. The Nyquist
frequencies are evaluated by an equation of the same form, but with one or the other index held
fixed (e.g., u = Nx − 1 while v = 0, ..., Ny − 2). The amplitude ẑ(kx = 0, ky = 0) at array element
(u, v) = (Nx /2 − 1, Ny /2 − 1) is usually set to zero, corresponding to the mean sea level being set
to 0. For generation of time-independent surfaces, we can set t = 0 so that the cosines are one and
the sines are zero, which cuts the number of terms to be evaluated in half.
If the frequencies are in the FFT order, then the last equation has the same general form, but
the indexing that expresses Hermitian symmetry is different: k(u) = −k(N − u), u = 1, ..., N − 1,

43
with k = 0 being a special case. The equation corresponding to Eq. (3.6) is then

2 ẑ(u, v, t)
q
p
= ρ(u, v) Ψ2S (u, v) + ρ(Nx − u, Ny − v) Ψ2S (Nx − u, Ny − v) cos[ω(k(u, v))t]
q
p
− σ(u, v) Ψ2S (u, v) + σ(Nx − u, Ny − v) Ψ2S (Nx − u, Ny − v) sin[ω(k(u, v))t]
(
p q
+ i ρ(u, v) Ψ2S (u, v) − ρ(Nx − u, Ny − v) Ψ2S (Nx − u, Ny − v) sin[ω(k(u, v))t]
)
p q
+ σ(u, v) Ψ2S (u, v) − σ(Nx − u, Ny − v) Ψ2S (Nx − u, Ny − v) cos[ω(k(u, v))t] . (3.7)

However, this equation does not explicitly show the special cases. Let zhat[u,v] be the array of
ẑ(k, t) = ẑ(kx (u), ky (v), t) values at a particular time t. r[u,v] and s[u,v] are the arrays of random
numbers ρ(u, v) and σ(u, v), respectively. Psiroot[u,v] is 21 Ψ2S (kx (u), ky (v)) (incorporating the 2
p

seen on the left-hand side of Eq. (3.7)). With other obvious definitions, the pseudo code to evaluate
Eq. (3.7) at a particular time t is as follows:
(1) Loop over the non-zero frequencies:

LOOP OVER u=1,Nx-1

LOOP OVER v=1,Ny-1
zhat[u,v] = COMPLEX(
( r[u,v]*Psiroot[u,v]
+ r[Nx-u,Ny-v]*Psiroot[Nx-u,Ny-v] )*cosomegat[u,v]
- ( s[u,v]*Psiroot[u,v]
+ s[Nx-u,Ny-v]*Psiroot[Nx-u,Ny-v] )*sinomegat[u,v] ,
( r[u,v]*Psiroot[u,v]
- r[Nx-u,Ny-v]*Psiroot[Nx-u,Ny-v] )*sinomegat[u,v]
+ ( s[u,v]*Psiroot[u,v]
- s[Nx-u,Ny-v]*Psiroot[Nx-u,Ny-v] )*cosomegat[u,v] )
END v LOOP
END u LOOP

(2) Loop over all ky values for the special case of kx = 0 at frequency array index u = 0. Note that
for kx = 0, the ±ky values are complex conjugates.

LOOP OVER v=1,Ny/2

zhat(0,v) = COMPLEX(
( r(0,v)*Psiroot[0,v]
+ r[0,Ny-v]*Psiroot[0,Ny-v] )*cosomegat[0,v]
- ( s[0,v]*Psiroot[u,v]
+ s[0,Ny-v]*Psiroot[0,Ny-v] )*sinomegat[0,v] ,
( r[0,v]*Psiroot[u,v]

44
- r[0,Ny-v]*Psiroot[0,Ny-v] )*sinomegat[0,v]
+ ( s[0,v]*Psiroot[u,v]
- s[0,Ny-v]*Psiroot[0,Ny-v] )*cosomegat[0,v] )
zhat(0,Ny-v) = CONJ( zhat(0,v) )
END u LOOP

(3) Loop over all kx values for the special cases of ky = 0 at frequency array index v = 0. Note
that for ky = 0, the ±kx values are complex conjugates.

LOOP OVER u=1,Nx/2

zhat(u,0) = COMPLEX(
( r(u,0)*Psiroot[u,0]
+ r[Nx-u,0]*Psiroot[Nx-u,0] )*cosomegat[u,0]
- ( s[u,0]*Psiroot[u,0]
+ s[Nx-u,0]*Psiroot[Nx-u,0] )*sinomegat[u,0] ,
( r[u,0]*Psiroot[u,0]
- r[Nx-u,0]*Psiroot[Nx-u,0] )*sinomegat[u,0]
+ ( s[u,0]*Psiroot[u,0]
- s[Nx-u,0]*Psiroot[Nx-u,0] )*cosomegat[u,0] )
zhat(Nx-u,0) = CONJ( zhat(u,0) )
END v LOOP

(4) Finally, set the (kx , ky ) = (0, 0) value to 0, which sets the mean sea level to zero:

zhat[0,0] = COMPLEX(0.0, 0.0)

Array zhat[u,v] = ẑ(kx (u), ky (v)) as just defined is Hermitian and has the frequencies in FFT order
ready for input to an FFT routine.

3.2.2 Example: A Two-Dimensional Sea Surface

As an example, let us use a course grid sampling of Nx × Ny = 64 × 64 points, which makes
it easier to see certain features in the associated plots. The physical region to be simulated is
Lx × Ly = 100 × 100 m. The two-dimensional, one-sided variance spectrum of Elfouhaily et al.
(1997) is used; the defining equations are given in §A.3.
The upper left plot of Fig. 3.3 shows this spectrum for a wind speed of 5 m s−1 and a fully
developed sea state. The contours are of Ψ1S (kx , ky ) evaluated at the discrete grid points; the
1
2 ∆kx ∆ky factor seen in Eq. (3.3) has not yet been applied to create a discrete two-sided spectrum.
The line along ky = 0 corresponds to the spectrum plotted in Fig. A.3. We see in both plots that
for a 5 m s−1 wind, the spectrum peaks near 0.2 − 0.3 rad/m.
The plots of the real and imaginary parts of ẑ(kx (u), ky (v)) show that most of the variance is at
low frequencies and that the amplitudes have the Hermitian symmetry illustrated previously in Figs.
3.1 and 3.2. The lower right panel of the plot shows a contour plot of the sea surface generated from
the inverse FFT of the amplitudes. The significant wave height for this surface realization is 0.60
m, in good agreement with the expected value given above for the Pierson-Moskowitz spectrum,
which is similar to the Elfouhaily et al. (1997) spectrum in the gravity wave region.

45
Figure 3.3: Example two-dimensional sea surface generated from a 2-D, one-sided variance spec-
trum. The resolution is only 64 × 64 grid points, so as to make the features of the underlying
spectrum and the Fourier amplitudes easier to see.

46
3.3 Resampling the FFT Grid for Ray Tracing
One of the most important uses of FFT-generated sea surfaces is Monte-Carlo ray tracing for
the computation of sea surface optical reflectance and transmittance. Such calculations require
repeatedly finding the point of intersection of a light ray with the random sea surface and finding
the angle of the ray from the normal to the surface at the point of intersection. Those calculations
are most easily performed if the sea surface is constructed as an array of triangular wave facets. The
three vertices of each wave facet define a plane, so that the ray tracing then amounts to finding the
intersection of a line with a plane. The direction of the normal to the surface is constant within each
triangular wave facet and is easily determined from the elevations at the vertices of the triangles.
Sea surfaces generated using FFT techniques are defined by a rectangular grid of surface ele-
vations. Ray tracing with such a surface would be greatly complicated by the fact that the four
points of each grid cell or wave facet do not in general define a plane, but rather a curved surface
whose form depends on the surface elevations at the four points. Moreover, the local normal to the
curved sea surface varies from point to point within each rectangular cell. It is therefore desirable
to “resample” the rectangular FFT grid to obtain a grid of triangles. This is easily done as shown
in Fig. 3.4. The blue lines of this figure show the rectangular FFT grid for the case of an 8 × 8 m
region of water surface with the origin of the x − y coordinate system placed at the center of the
region. The FFT grid has Nx = 16 and Ny = 8 points. The red dots show a hexagonal grid of
triangles defined using every other x grid point and every y point of the FFT grid; one of these
triangles is colored in red. It is clear that this mapping requires that Ny = Nx /2. The number of
N N
triangle vertices defined in this way is Nvertices = 3 2y ( 2y + 1) + 1, and the number of triangular
wave facets generated from these vertices is Nfacets = 32 Ny2 .
It may seem wasteful to compute the surface elevation at each of the FFT grid points (the points
where the blue lines intersect) and then discard roughly half of these points when resampling the
FFT grid to obtain a hexagonal grid of triangular wave facets. However the extremely fast run time
of the FFT algorithm makes this option much more efficient than the regular Fourier technique
given in Mobley (1994, his Eq. 4.77). That technique computes the surface elevations only at
the needed triangle vertices, but requires much more time to evaluate the 2-D discrete Fourier
transforms than the FFT technique described here.
Figure 3.5 shows an example of a sea surface defined on a rectangular FFT grid and the corre-
sponding hexagonal patch of sea surface defined by a grid of triangular wave facets derived from a
rectangular grid as in Fig. 3.4. The origin of the coordinate system has been placed at the center
of the Lx × Ly = 100 × 100 m region. The hexagonal-grid sea surface is ready for use in ray tracing
computations with existing algorithms. Although those computations are conceptually simple, the
actual calculations are quite tedious. Mobley (2014, Appendix B) gives the details of a ray tracing
algorithm for a sea surface defined by a hexagonal grid of triangular wave facets. (That report also
gives the details of how ray tracing is performed for polarized light.)

47
Figure 3.4: Example mapping of a rectangular computational grid used in FFT simulations of a
random sea surface (blue lines) to a hexagonal grid of triangles as used in ray tracing (red dots).
The FFT grid points at the right and top, shown by the dotted blue lines, are obtained by the
inherent spatial periodicity of the surface as determined by FFT techniques. Thus the surface
elevation of point 35 is the same as that of point 27, point 57 is the same as point 1, etc. One of
the triangles generated from the FFT grid is shaded in red.

48
Figure 3.5: The blue lines in the top panel show a sea surface specified on a rectangular FFT grid.
The red dots show the points used to define an hexagonal grid of triangles as in Fig. 3.4. The
bottom figure shows the surface specified by the hexagonal grid of triangular wave facets shown by
the red dots in the top figure.

49
3.4 Example: Reflectance by Different Surface Models
We have now seen how to correctly generate 2-D surface realizations. The previous chapter made
two important points that underlie surface generation. First, in §2.5.1 we pointed out two errors
in the widely used equations of Tessendorf (2004). His version of our Eq. (2.7) does not divide
the one-sided energy spectrum by 2 to convert it to a two-sided spectrum, and his version of Eq.
(2.9) does not have the overall scale factor of √12 . Omission of these factors in either the 1-D
equations (2.7) and (2.9), or in the 2-D equivalents (3.3) and 3.4), makes the Fourier amplitudes ẑ
a factor-of-two too large. Consequently the wave amplitudes (surface elevations) will be too large
and energy, which is proportional to wave amplitude squared, will not be conserved by a factor of
four in a round-trip calculation like that of §2.5.2. Second, in §2.6 we pointed out the importance
of resolving the sea surface slopes when the generated surfaces are used for optical calculations, and
we presented a simple way to rescale variance spectra so that the generated surfaces adequately
reproduce both the elevation and slope variances of real sea surfaces, even though the elevation
variance spectrum is sampled over a limited range of spatial frequencies.
Figure 3.6 shows the consequences for optical calculations of either omitting the scale factors
on the amplitude equations or failing to resolve the slope variance. Polarized ray tracing was
performed as described in Mobley (2014) to compute the percentage of the total incident energy
reflected into all directions for the sun at a given zenith angle in an otherwise black sky. The wind
speed was 10 m s−1 , and the surface was modeled using Lx = Ly = 200 m and Nx = 1024, Ny = 512
as required for the resampling of the FFT grid as described in §3.3. Rays were traced for 50,000
independent surface realizations to obtain average reflectances with negligible Monte Carlo noise.
The black curve shows the energy reflectance (without regard for the state of polarization) for a level
sea surface; this reflectance is given by the well known Fresnel reflectance equations (e.g. Mobley,
1994, Eq. 4.14). The solid green curve shows the reflectance for FFT-generated surfaces when the
correct amplitude equations (3.3) and (3.4) are used, and when the slope variance is fully resolved
via the δ(k) function of Eqs. (2.15) and (2.16) for the given values of Lx , Ly , Nx and Ny . The
dashed green curve shows the reflectance for surfaces generated using the well known Cox-Munk
wind speed–wave slope statistics as described in Mobley (2014, §4.3). Cox-Munk surfaces model
the surface slopes, but not the surface wave elevations. We see that the FFT-generated surface is
in excellent agreement with the Cox-Munk surfaces. The agreement is good because both surfaces
resolve the slope statistics. (Calculations not shown here indicate that Cox-Munk and FFT surfaces
do differ by a few percent at large incident angles for other wind speeds.)
The red curve shows the reflectance computed when the correct amplitude equations (3.3) and
3.4) are used, but the 2-D elevation variance spectrum is sampled only over the spatial frequency
ranges determined by the values of Lx , Ly , Nx and Ny . That is, the slope variance is not fully
resolved. Because the undersampled slope variance is less than the true slope variance, the generated
surfaces are too smooth for the high spatial frequencies (short wavelengths), and the reflectance
is too large. The difference is significant at sun angles near the horizon. The blue curve shows
the reflectance when the slope variance is fully resolved, but the Tessendorf equations are used
to generate the Fourier amplitudes. The surface is now too rough (the amplitudes are too large),
which reduces the reflectance is the same manner as a higher wind speed does. (The reflectance at
large incident angles decreases with increasing wind speed because of increased surface roughness.)

50
Figure 3.6: Energy reflectance by different sea surface models for a wind speed of 10 m s−1 .

Although ray tracing and optics are beyond the scope of this tutorial, Fig. 3.6 is included here to
emphasize the importance of properly modeling the sea surface when doing calculations of surface
optical properties.

51
CHAPTER 4

Two-Dimensional, Time-Dependent Surfaces

One final step remains: the addition of time dependence to a sequence of the sea surface real-
izations. Many scientific applications do not need time dependence. This is the case when many
independent random realizations of sea surfaces are used for Monte Carlo ray tracing to compute
the average optical reflectance and transmittance properties of wind-blown sea surfaces. However,
time dependence is crucial for applications such as computer-generated imagery for movies.
We already have the needed theory in hand. The fundamental Eqs. (3.3) and (3.4), and their
expansions such as Eq. (3.7) include time dependence. Most importantly, these equations made no
simplifying assumptions about the ±k symmetry of the two-sided variance spectrum Ψ2S (k).
As we have learned, the Fourier transform of a snapshot of the sea surface gives a two-sided
variance spectrum with identical values for −k and +k. This represents equal amounts of energy
propagating in the −k and +k directions, i.e., equal amounts of energy in waves propagating upwind
and downwind. Such a situation in nature gives standing waves. Here also, if Ψ2S (−k) = Ψ2S (k),
the time-dependent surface will be standing waves composed of waves of equal energy propagating
the +k and −k directions. To obtain waves propagating downwind, as is the case in a real ocean,
we must use an asymmetric two-sided spectrum with Ψ2S (−k) << Ψ2S (+k), so that almost all
of the energy is propagating downwind. Note, howerver, that we cannot simply set Ψ2S (−k) = 0,
which represents no energy at all propagating upwind, because Ψ2S (−k) = 0 destroys the Hermitian
property of Eq. (3.4). Thus we must conjure up an asymmetric variance spectrum that allows a
nonzero (but perhaps negligible) amount of energy to propagate upwind.
An asymmetric two-sided variance spectrum can be constructed using an asymmetric spread-
ing function Φ(k, ϕ) as described in §A.3.1. The spreading functions of the form Φ(k, ϕ) =
Cs cos2s(k) (ϕ/2) described there allow some energy to propagate in −k directions, i.e. when
|ϕ| > π/2. With this choice of a spreading function, we can let the magnitude of Ψ2S (+k) equal
the magnitude of the one-sided variance function Ψ1S (+k), which gives the total variance, because
we assume that a negligible amount of the total energy propagates in −k directions. With this
assumption, there is no division of Ψ1S (+k) by 2 in the first line of Eq. (3.3).
Once an asymmetric Ψ2S (±k) has been defined, we can evaluate Eq. (3.3) for a set of random
numbers ρ(kuv ) and σ(kuv ). This is done only once. Then to generate a sequence of sea surface
realizations for times t = 0, ∆t, 2∆t, ..., we multiply the time-independent amplitudes by the time-

52
dependent cosines and sines as shown in Eq. (3.7). We thus obtain the amplitudes ẑ(kuv , t) at the
current time t. The inverse Fourier transform of this ẑ(kuv , t) gives the sea surface z(xrs ) at time
t.
The physics (or lack thereof) of this process is simple. We start with a realization of the sea
surface at time zero. This surface contains waves of many discrete frequencies kuv traveling in all
directions ϕuv . Then to get the surface at any later time t, we simply propagate the sinusoidal waves
of each frequency kuv in their original direction of travel through a phase change determined by the
time step and the dispersion relation ω(kuv ). For deep-water gravity waves, the dispersion relation
√
is ω(kuv ) = gkuv . It thus visually appears that the waves are propagating and the sea surface
profile is changing with time. However, this Fourier transform technique is really just moving a
collection of independent, non-interacting sinusoids through frequency-dependent phase angles to
create a new surface realization from the sum of the phase-shifted sinusoids. In a real ocean, waves
of difference frequencies can interact with each other (redistribute energy among frequencies) in
highly complex and nonlinear ways, so that the sea surface statistics may not be time-independent.
This is, in particular, how little waves grow to big waves when the wind begins to blow over a
level surface. Modeling the nonlinear evolution of a sea surface is beyond the abilities of Fourier
transform techniques which are, at heart, just a mathematical artifice based on a time-independent
directional variance density spectrum.
Figure 4.1 shows a sequence of surface realizations for a 10 m/s sec wind speed and a spatial
grid of size Lx × Ly = 100 × 100 m, and a grid resolution of Nx × Ny = 256 × 256. A time series of
images made with such a course grid would be useful for many purposes.
We have seen that the spatial pattern of a sea surface generated by a Fourier transform is
periodic. This is convenient for tiling a small computed region into a visually acceptable larger
region. When time dependence is included and a finite-length time series of images is created, the
sequence of images is not periodic in time because the various sinusoids comprising the surface do
not have a common period. As pointed out by Tessendorf (2004, §4.2), this can be remedied by
“quantizing” the temporal frequency ω(kuv ) as follows.
Let Tr be the length of time over which the time sequence of surface realizations should repeat.
The number of frames Nt in the video loop determines the time step between frames, ∆t = Tr /Nt .
For a smoothly moving sea surface, Nt must be large enough that ∆t is less than 0.1 s. Define
√
ωo ≡ 2π/Tr . For deep-water gravity waves, the true temporal frequency ω(kuv ) = gkuv is replaced
by $√ %
gkuv
ω̃(kuv ) = ωo , (4.1)
ωo
where bxc converts a real number x into its integer part (e.g., 15.23 is converted to 15). This
operation slightly alters the temporal frequency of each wave component so that each component
returns to exactly its initial shape after time Tr . A video loop can then be created from the sequence
of surfaces, and the loop will match perfectly when the frame for time Tr − ∆t goes to the frame
for time Tr , which is the same surface as time 0, after which the surfaces repeat. This adjustment
to ω is greatest for the lowest frequencies, but the adjustment becomes smaller and smaller as the
repeat time becomes larger and larger. It is thus easy to create a time-dependent small area of sea
surface that can be tiled in both space and time to create an image of a larger sea surface over a
longer time. This is good enough to fool movie-goers who have not studied this tutorial.

53
Figure 4.1: A time-dependent sequence of 2-D sea surfaces for a 10 m s−1 wind speed. The physical
domain is 100 × 100 m; the sampling is 256 × 256 points. The vertical scale of the plots is -3 to +3
m. The wind is blowing in the +x direction, which is from the upper left to the lower right of each
figure.

54
Finally, we note that in order to employ a re-scaled variance spectrum as described in §2.6, we
determine the value of the δNy parameter using the omnidirectional variance spectrum, as seen in
Eq. (2.17). We then apply the δ(k) correction to the directional spectrum Ψ(kx , ky ) with k = kuv
for each (kx (u), ky (v)) point of the rectangular grid.
An example of a 20-second (simulated time) video loop created using all of these tricks can be
seen at
http://www.oceanopticsbook.info/view/surfaces/level_2/twodimensional_timedependent_
surfaces. Figure 4.2 shows the first frame of that video. The wind speed is 10 m s−1 and the
omnidirectional variance spectrum is that of Elfouhaily et al. as described in §A.3, combined with a
cosine-2s spreading function as described there. The simulated region is 100×100 m using 512×512
grid points. Note in this figure that the mean square slopes compare well with the corresponding
Cox-Munk values shown in Table 4.1. The mss values for the generated surface are computed from
finite differences, e.g.
z(r + 1, s) − z(r, s)
mssx (r, s) =
x(r + 1) − x(r)
averaged over all points of the 2-D surface grid. The hθx i and hθy i values shown in the figure are the
average angles of the surface from the normal to a level surface (the +z axis). These are computed
from
θx (r, s) = tan(mssx (r, s)) ,
etc., averaged over all points of the surface. The mssx values of this 2-D image can be compared
with those of the 1-D simulation of Fig. 2.7.

slope variable FFT value Cox-Munk formula Cox-Munk value

mssx 0.032 0.0316U 0.032
mssy 0.021 0.0192U 0.019
mss 0.053 0.001(3 + 5.12U ) 0.054

Table 4.1: Comparison of Cox-Munk mean square slopes and values for the first frame of the video.

55
Figure 4.2: The first frame of the video loop.

56
APPENDIX A

Wave Spectra

The entire business of wave spectra can be exceedingly confusing. Just as for the FFT techniques
to generate sea surfaces, this is well known material. However, journal articles always assume that
the reader already understands the underlying ideas and mathematics. It is therefore worthwhile to
review the results needed for this tutorial. The best-written reference I’ve found for the development
of wave spectra is Holthuijsen (2007), who is very careful in his notation and terminology. The
notation used below is chosen to conform to that used in Elfouhaily et al. (1997), which will often
be referred to as the ECKV spectrum, after the last letters of the authors’ names.
The figures of the previous chapters were illustrated using two commonly used wave variance
spectra. The one-dimensional surfaces of §2.5 used the 1-D Pierson-Moskowitz spectrum, and the
two-dimensional surfaces of §3.2.2 and Chapter 4 used the 2-D spectrum of Elfouhaily et al. (1997).
After the introductory overview in the next section, the remaining two sections present these two
specific wave spectra for easy reference.

A.1 Overview of Wave Spectra

As we saw in §2.1 and 2.2, there is a close connection between the variance of the sea surface waves
and the energy per unit area contained in the waves. For a discretely sampled sea surface with
elevations z(r), the discrete Fourier transform (1.14) gives the corresponding amplitudes ẑ(u). The
discrete variance spectrum is then the function S(u) = |ẑ(u)|2 , u = 0, ..., N − 1, which is often (but
not here!) loosely referred to as the discrete energy or power spectrum.
Care is required to formulate wave spectra for continuous variables. As in Eqs. (1.2) and (1.6),
we can write the continuous surface as a sum of sinusoids, e.g.
∞
X ∞
X
z(x) = zn (x) = An cos(kn x + φn ) . (A.1)
n=0 n=0

As we saw in §2.1, the variance of the sinusoid with frequency kn is 21 A2n . The waves of different
frequencies are independent, so the total variance of the surface can be written as the sum of the

57
variances of the individual waves:
∞
X 1
var{z} = A2n . (A.2)
2
n=0

Now let ∆kn be a frequency interval centered on frequency kn , whose sinusoid has amplitude An .
We then define
1 2
A
S(k = kn ) ≡ lim 2 n . (A.3)
∆kn →0 ∆kn

In this definition, we keep in mind that each An is associated with a particular frequency kn , and
that the limit operation holds for each value of n. We are thus defining a function of the spatial
frequency, which becomes a continuous function of k as the bandwidth ∆kn goes to zero.
The quantity S(k) is called the omnidirectional variance spectrum. “Omnidirectional” means
that there is no reference direction (e.g., a direction relative to the wind direction) included in the
quantity. This is the case if a wave record is made at a single point as a function of time: the waves
go past and their elevations are recorded, but no information is obtained on the direction the waves
are traveling. S(k) is also called the omnidirectional elevation spectrum for obvious reasons. As is
to be expected, there is no uniformity of notation for this spectrum, but S seems to be the most
common symbol—and what is used in both Pierson and Moskowitz (1964) and Elfouhaily et al.
(1997) in the next sections—so that what I will use here. (E seems to be the second-most-common
symbol and is used in Holthuijsen (2007) and several other papers on my desk.) The units of
S(k) are clearly m2 /(rad/m). Equations (A.2) and (A.3) show that integrating the omnidirectional
variance spectrum over all frequencies gives the total variance:
Z ∞
var{z} = hz 2 i = S(k)dk .
0

(The equations above are written in terms of spatial angular frequency k, as used for surface
generation, but the reasoning and functional form of Eq. (A.3) are the same for any other spatial
or temporal frequency variable.)
We can now repeat this process for two dimensions, starting with
∞ X
X ∞
z(x, y) = An,m cos(kn x + km y + φn ) ,
n=0 m=0

and ending up with a function

1 2
2 An,m
Ψ(kx , ky ) ≡ lim lim . (A.4)
∆kn →0 ∆km →0 ∆kn ∆km
Here the x and y subscripts on k remind us that there are independent spatial frequencies associated
with the x and y directions. Ψ(kx , ky ) is the directional variance spectrum in Cartesian coordinates.
Its units are m2 /(rad/m)2 . This spectrum is often called the “two-dimensional wavenumber spec-
trum,” and its arguments are often written in vector form, Ψ(k), where k = (kx , ky ) denotes the
location of the frequency point in the 2-D frequency plane. This spectrum depends on direction,
i.e., on the direction of the (kx , ky ) point in a two-dimensional frequency plane. Usually, the +x
direction is chosen to be pointing downwind. In this case, the angle ϕ = tan−1 (ky /kx ) gives the

58
direction relative to the wind direction. At Eq. (A.4) shows, the integral of Ψ(kx , ky ) over all
frequencies gives the variance of the two-dimensional surface:
Z ∞Z ∞
var{z} = hz 2 i = Ψ(kx , ky ) dkx dky .
−∞ −∞

It is also common to define a directional spectrum in terms of polar coordinates given by the
magnitude k and direction ϕ of the vector k. These are are related to kx , ky by
q
k = kx2 + ky2

−1 ky
ϕ = tan
kx

and inversely by

kx = k cos ϕ
ky = k sin ϕ .

In this case we define

1 2
2 An,m
Ψ̃(k, ϕ) ≡ lim lim . (A.5)
∆k→0 ∆ϕ→0 ∆k∆ϕ
This spectrum has units of m2 /[(rad/m) rad]. (The tilde notation is used here to distinguish this
spectrum from the Ψ(k, ϕ) spectrum of Elfouhaily, which is defined below. Some authors reserve
the name “directional spectrum” for Ψ̃(k, ϕ) and refer to Ψ(kx , ky ) as the wavenumber spectrum.)
As before, definition (A.5) shows that integrating Ψ̃(k, ϕ) over k and ϕ gives the variance:
Z ∞ Z 2π
var{z} = Ψ̃(k, ϕ) dk dϕ .
0 0

The directional spectrum given in §A.3 is specified in terms of polar coordinates k, ϕ. However,
we need a spectrum in terms of Cartesian coordinates kx , ky for use in a rectangular FFT grid. The
change of variables from polar to Cartesian coordinates is effected by the Jacobian

∂(k, ϕ)
Ψ(kx , ky ) = Ψ̃(k, ϕ)

∂(kx , ky )

∂k ∂k
∂kx ∂ky
= Ψ̃(k, ϕ) ∂ϕ ∂ϕ
∂kx ∂ky
1
= Ψ̃(k, ϕ) .
k

Note that the 1/k factor converts the units of Ψ̃(k, ϕ) into the units of Ψ(kx , ky ).
In Eq. (A.15) of §A.3 this last equation is partitioned as
1
Ψ(kx , ky ) = S(k)Φ(k, ϕ) ≡ Ψ(k, ϕ) , [ECKV 45] (A.6)
k
where S(k) is an omnidirectional spectrum and Φ(k, ϕ) is a nondimensional spreading function,
which shows how waves of different frequencies propagate (or “spread out”) relative to the downwind

59
direction at ϕ = 0. Labels such as [ECKV 45] refer to the corresponding equations in Elfouhaily
et al. (1997). The spreading function by definition satisfies
Z 2π
Φ(k, ϕ) dϕ = 1 (A.7)
0

for all k.
Equation (A.6) shows that to obtain the ECKV variance spectrum in Cartesian coordinates we
need only evaluate the ECKV Ψ(k, ϕ) spectrum for the corresponding values of k and ϕ, i.e.
q
Ψ(kx , ky ) = Ψ k = kx2 + ky2 , ϕ = tan−1 (ky /kx ) . (A.8)

Note in particular that there is no “extra” k factor involved in the conversion of Ψ(k, ϕ) to Ψ(kx , ky );
both quantities have the same units.
Integration of Eq. (A.6) over the respective (kx , ky ) and (k, ϕ) frequency planes gives the
variance as
Z ∞Z ∞
2
hz i = Ψ(kx , ky ) dkx dky
−∞ −∞
Z ∞ Z 2π
1
= S(k) Φ(k, ϕ) k dk dϕ
k
Z0 ∞ 0

= S(k) dk , [ECKV A2] (A.9)

0

after noting the normalization of Eq. (A.7). Thus the variance of the sea surface is still contained
in the omnidirectional spectrum, even in the two-dimensional case. The omnidirectional spectrum
S(k) is obtained from Ψ(k, ϕ) via
Z π
S(k) = Ψ(k, ϕ) k dϕ . [ECKV A3]
−π

Now return to Eq. (A.1) and take the derivative to get the slope of the sea surface for the nth
wave:
dzn (x)
= −An kn sin(kn x + φn )
dx
As in Eq. (2.1), the variance of this slope is
Z Λn
dzn 1 1
var ≡ [An kn sin(kn x + φn )]2 dx = A2n kn2 .
dx Λn 0 2
A limit operation corresponding to Eq. (A.3) gives
1 2 2
2 An kn
lim = k 2 S(k) . (A.10)
∆kn →0 ∆kn
The quantity k 2 S(k) is the omnidirectional slope variance spectrum, usually called just the slope
spectrum. The units of k 2 S(k) are m rad. Integrating the slope spectrum over all frequencies gives
the total variance σ 2 of the sea surface slope:
* 2 + Z ∞
2 dz dz
σ ≡ var = = k 2 S(k)dk .
dx dx 0

60
This variance is usually called the mean square slope or mss. The units of mss are rad2 . This is
a nondimensional number, but the label of rad2 reminds us that we can think of the slope as an
angle from the horizontal measured in radians.
The corresponding relations for two dimensions are derived in the same fashion and lead to
similar results. Assuming that the wind is blowing in the +x direction, the mean-square slope in
the along-wind direction is given by either of
D ∂z(x, y) 2 E Z ∞ Z ∞
≡ σx2 ≡ mssx = kx2 Ψ(kx , ky ) dkx dky
∂x
Z−∞ −∞
∞ Z π
= k 2 cos2 ϕ Ψ(k, ϕ) k dk dϕ . [ECKV A4]
−∞ −π

The corresponding equation for the cross-wind direction is

D ∂z(x, y) 2 E Z ∞ Z ∞
≡ σy2 ≡ mssy = ky2 Ψ(kx , ky ) dkx dky
∂y
Z−∞ −∞
∞ Z π
= k 2 sin2 ϕ Ψ(k, ϕ) k dk dϕ . [ECKV A5]
−∞ −π

Recalling that variances of random variables add to get the total variance due to all variables gives
the total mean square slope
Z ∞Z ∞
mss = mssx + mssy = (kx2 + ky2 ) Ψ(kx , ky ) dkx dky
−∞ −∞
Z ∞Z π
= k 2 Ψ(k, ϕ) k dk dϕ
Z−∞∞
−π

= k 2 S(k) dk . [ECKV A6]

−∞

Thus, even in the 2-D case, the total slope variance can be obtained from the 1-D omnidirectional
slope spectrum.
Table A.1 summarizes the spectral quantities used in this tutorial.

Spectrum Name Symbols Units

1-D or omnidirectional
variance or elevation S(k) m2 /(rad/m)
slope k 2 S(k) m rad
2-D or directional
variance or elevation Ψ(kx , ky ), Ψ(k, ϕ) m2 /(rad/m)2
alongwind slope kx2 Ψ(kx , ky ), k 2 cos2 ϕ Ψ(k, ϕ) rad2
crosswind slope ky2 Ψ(kx , ky ), k 2 sin2 ϕ Ψ(k, ϕ) rad2
total slope (kx2 + ky2 ) Ψ(kx , ky ), k 2 Ψ(k, ϕ) rad2

Table A.1: Summary of spectral quantities.

61
A.2 The Pierson-Moskowitz Omnidirectional Gravity Wave
Spectrum
The Pierson-Moskowitz spectrum (Pierson and Moskowitz, 1964) describes gravity waves in a “fully
developed” sea. A fully developed sea is an idealization of the statistically steady-state wave field
resulting from a steady wind blowing for an infinitely long time over an infinite fetch. (In practice,
the duration and fetch required to achieve something close to a fully developed sea depend on
the wind speed. A steady wind of 5 m s−1 blowing for 10 hours over a fetch of 60 km might be
adequate; for hurricane winds of 35 m s−1 , a fetch of a few thousand kilometers and a duration of
several days would be required. Thus it is much easier to approach a fully developed sea at low
wind speeds than at high.)
The one-sided, omnidirectional Pierson-Moskowitz spectrum, formulated in terms of angular
spatial frequency k, is
g 2 1
α
SPM (k) = 3 exp −β 4 [m2 /(rad/m)] , (A.11)
2k k U19

where

α = 0.0081,
β = 0.74,
g = 9.82 m s−2 is the acceleration of gravity,
U19 is the wind speed in m s−1 at 19.5 m above the sea surface, and
k is the angular spatial frequency in rad m−1 .

The wind speed at 19.5 m can be obtained from the more commonly used wind at 10 m above the
sea surface by the approximate formula

U19 ≈ 1.026U10 .

As has already been noted, it is often of interest to express a variance spectrum in terms of other
variables, such as the wavenumber ν or the temporal angular frequency ω. To change variables
in a spectral density function, the key is to recall that a variance density function by definition
expresses the variance per unit frequency interval. The variance contained in some interval dk of
the spatial angular frequency equals the variance contained in the corresponding interval dν of the
wavenumber or the interval dω of the temporal frequency. Thus we have

SPM (k)dk = SPM (ν)dν = SPM (ω)dω .

To change the variable from k = 2πν to ν, the previous equation gives

dk
SPM (ν) = SPM (k) = SPM (k = 2πν)2π ,
dν
which results in g 2 1
α
SPM (ν) = 2 3 exp −β 4 [m2 /(1/m)] . (A.12)
8π ν 2πν U19

62
Figure A.1: The Pierson-Moskowitz variance spectrum as functions of k and ω for wind speeds of
U10 = 5, 10, and 15 m s−1 . The vertical dotted lines at k = 370 rad/m and ω = 60.3 rad/s show
the boundary between gravity and capillary waves.

To change variables from spatial angular frequency k to temporal angular frequency ω, we use the
dispersion relation for gravity waves in deep water,

ω 2 = gk ,

to evaluate
dk 2ω
= ,
dω g
which leads to " 4 #
αg 2

g
SPM (ω) = 5 exp −β [m2 /(rad/s)] . (A.13)
ω ωU19
All of these formulas have units of meters squared over the appropriate frequency. (The quantities
dk/dν and dk/dω seen in the conversions are the Jacobians for the one-dimensional changes of
variables.) Figure A.1 shows the Pierson-Moskowitz spectrum as functions of k and ω for wind
speeds of U10 = 5, 10, and 15 m s−1 .
This spectrum has withstood the test of time quite well (e.g., Alves and Banner, 2003) as a
description of gravity waves in fully developed seas. However, it should not be used for high-
frequency (short-wavelength) gravity waves, and certainly not for capillary waves. Likewise, it does
not describe young seas, which have not had the time or fetch needed to approach the state of a
well developed sea.
Figure A.2 shows the Pierson-Moskowitz slope spectra for three wind speeds. Note that the
slope spectrum falls off much more slowly for high frequencies than does the elevation spectrum.
That means that the higher frequencies contribute much more to the total slope variance than they
do to the total elevation variance.

63
Figure A.2: The Pierson-Moskowitz slope spectrum for wind speeds of U10 = 5, 10, 15 m s−1 . The
plot uses the same ordinate scale as used for the elevation spectrum in the left panel of Fig. A.1 in
order to highlight the slow falloff of the slope spectrum compared to the elevation spectrum. The
vertical dotted line is the boundary between gravity and capillary waves.

A.3 The Elfouhaily et al. Directional Gravity-Capillary

Wave Spectrum
In order to model two-dimensional sea surfaces z(x, y), we need a wave variance spectrum that
describes the distribution of wave variance for waves propagating in difference directions.
As the examples in §3.1 showed, if we have measurements of the two-dimensional sea surface
elevation at a given time, z(xrs ) = z(xr , ys ) = z(r, s), then the two-dimensional discrete Fourier
transform of Eq. (1.17) gives the two-dimensional amplitudes

ẑ(kuv ) = ẑ[kx (u), ky (v)] = ẑ(u, v) = F{z(r, s)} = F{z(xrs )} .

Dividing by the discrete frequency bandwidths gives an estimate (a 2-D periodigram) of the two-
dimensional variance spectral density

|ẑ(kuv )|2 Ψ(kuv )

P (kuv ) ≡ = , (A.14)
∆ku ∆kv ∆ku ∆kv
just as in the 1-D case. Here I have written z(xrs ), ẑ(kuv ), and Ψ(kuv ) to emphasize that these
are discrete functions, whereas Ψ(k) denotes a spectral density function of the continuous variable
k. Just as in the 1-D case, the 2-D periodogram is an approximation of the 2-D variance spectral
.
density, P (kuv ) = Ψ(k). As Eq. (1.17) shows, Ψ(kuv ) has units of m2 , whereas P (kuv ) and Ψ(k)
have units of m2 /(rad/m)2 .
Unfortunately, making measurements of 2-D sea surfaces is extremely difficult. There are very
few such data sets obtained, for example, by laser reflectance measurements (e.g., Huang et al.,

64
2000), and these do not cover the full range of spatial scales. Given the paucity of empirical 2-D
wave data from which to develop 2-D variance spectra, the common procedure is to start with a
1-D or omnidirectional spectrum and add an angular spreading function to distribute the wave
energy over different directions relative to the downwind direction. In nature, most waves travel
more or less downwind, a small amount of energy is contained in waves propagating in nearly
cross-wind directions, and almost no energy is contained in waves that by some chance (such as
the breaking of a larger wave generating smaller waves in all directions) might be propagating in
upwind directions. The spreading function must capture this behavior. Although omnidirectional
wave spectra are well grounded in observations, the choice of a spreading function is something of
a black art.
Elfouhaily et al. (1997) presented an omnidirectional variance spectrum that explicitly includes
both the gravity and capillary wave scales. For brevity, I will often denote this paper and their model
by “ECKV,” taken from the initials of the authors’ last names. The boundary between gravity and
p
capillary waves is taken to be k = ρg/τ = 370 rad/m, the angular spatial frequency at which
the restoring forces (which tend to bring a wave surface back to an unperturbed level surface) of
gravity and surface tension are equal. The corresponding wavelength is Λ = 2π/370 = 0.017 m.
They then combine their omnidirectional spectrum with a spreading function to obtain their one-
sided, directional variance spectrum. Using their notation, the one-sided, 2-D ECKV spectrum has
the form
1
Ψ(k, ϕ) = S(k)Φ(k, ϕ) [ECKV 45] . (A.15)
k
Here S(k) is the 1-D omnidirectional spectrum with units of m2 /(rad/m), and Φ(k, ϕ) is a non-
dimensional spreading function. Ψ(k, ϕ) thus has units of m2 /(rad/m)2 . Equation labels such as
[ECKV 45] give for reference the corresponding equation in the ECKV paper.
The ECKV omnidirectional spectrum is
Bl + Bh
S(k) = [ECKV 30] , (A.16)
k3
where Bl is the low-frequency (long gravity wave) contribution to the variance, and Bh is the high-
frequency (short gravity wave to capillary wave) contribution. (The quantity k 3 S(k) is called the
curvature or saturation spectrum and is of interest in physical oceanography because it is related
to the rate of variance dissipation of the waves. Thus ECKV refer to Bl and Bh as the low and
high frequency curvature spectra.) The components of the omnidirectional spectrum are given by

LP M = exp[−1.25(kp /k)2 ] [ECKV 2]

1
q
Γ = exp{− 2 [( k/kp − 1)2 ]} [below ECKV 3]
2σ
Jp = γ Γ [ECKV 3]
q
Fp = LP M Jp exp{−0.3162Ωc ( k/kp − 1)} [ECKV 32]
Fm = LP M Jp exp[−0.25(k/km − 1)2 )] [ECKV 41]
(Note : A typo in ECKV Eq. 41 omitted the LP M Jp factor in Fm )
Bl = 0.5αp (cp /c)Fp [ECKV 31]
Bh = 0.5αm (cm /c)Fm [ECKV 40]

65
where

α = 0.0081,
β = 1.25,
g = 9.82 m s−2 is the acceleration of gravity,
U10 is the wind speed in m s−1 at 10 m above the sea surface
k is the angular spatial frequency in rad m−1
Ωc is defines the age of the waves for the given wind speed:

= 0.84 for a fully developed sea (corresponds to Pierson-Moskowitz)

= 1 for a ”mature” sea [used in ECKV Fig 8a]
= 2 to 5 for a “young” sea; the maximum allowed value is 5

Cd10N = 0.00144 is a drag coefficient [value deduced from ECKV Fig 11]
√
u∗ = Cd10N U10 is the friction velocity [using ECKV 61[
ao = 0.1733 (ECKV 59)
ap = 4.0
km = 370.0 rad/m
cm = 0.23 m/s is the phase speed of the wave with spatial frequency km
am = 0.13u∗ /cm [ECKV 59]
γ = 1.7 if Ωc ≤ 1 else γ = 1.7 + 6 log10 (Ωc )
σ = 0.08(1 + 4Ω−3
c )

αp = 0.006Ωc0.55 [ECKV Eq. 34]

αm = 0.01[1 + ln(u∗ /cm )] if u∗ ≤ cm else αm = 0.01[1 + 3 ln(u∗ /cm )] [ECKV 44]
2
ko = g/U10
kp = ko Ω2c is the spatial frequency of the maximum of the spectrum
p
cp = g/kp is the phase speed of the wave with spatial frequency kp
p
c = (g/k)(1 + (k/km)2 ) is the phase speed of the wave

At the lower frequencies, the ECKV spectrum is essentially the Pierson-Moskowitz spectrum
(the LP M term above) with an enhancement (the Jp term) that adds more energy to the lower fre-
quencies. The highest frequencies have a cutoff due to viscosity damping out the smallest capillary
waves. The ECKV omnidirection elevation and slope spectra are illustrated in Fig. A.3 for the
case of a fully developed sea and three wind speeds. Figure A.4 shows the spectra as a function of
wave age for a wind speed of U10 = 10 m s−1 .

66
Figure A.3: The omnidirectional part S of the Elfouhaily et al. (1997) elevation variance spectrum
(left panel) and slope spectrum k 2 S (right panel) for fully developed seas and wind speeds of
U10 = 5, 10, 15 m s−1 . The gray lines show the corresponding Pierson-Moskowitz spectra from Fig.
A.2.

Figure A.4: The omnidirectional part of the Elfouhaily et al. (1997) elevation spectrum (left panel)
and slope spectrum (right panel) for a wind speed of U10 = 10 m s−1 and wave ages from very
young (Ωc = 5) to mature (Ωc = 1) to fully developed (Ωc = 0.84). Compare with Fig. A.3.

67
A.3.1 Spreading Functions
The ECKV spreading function is given by
1
Φ(k, ϕ) = [1 + ∆(k) cos(2ϕ)]
2π
1
1 + tanh ao + ap (c/cp )2.5 + am (cm /c)2.5 cos(2ϕ)

= (A.17)
2π
Note that this function is symmetric about ϕ = π/2; i.e., the function has symmetric spreading
downwind and upwind. This is consistent with a symmetric variance spectrum Ψ(−k) = Ψ(k) as
would be obtained from a the Fourier transform of a snapshot of a sea surface.
A commonly used family of alternate spreading functions is given by the “cosine-2s” functions
(Longuet-Higgins et al., 1963), which have the form

Φ(k, ϕ) = Cs cos2s (ϕ/2) , (A.18)

where the normalizing coefficient is

1 Γ(s + 1)
Cs = √ ,
2 π Γ(s + 1/2)

and s is a spreading parameter that in general depends on k, U10 , and wave age. In this equation
R∞
Γ is the customary gamma function defined by Γ(p) ≡ 0 xp−1 e−x dx where p > 0. The cosine-2s
functions are asymmetric, with much stronger downwind than upwind propagation.
The ECKV and cosine-2s spreading functions are illustrated in Fig. A.5. Both of these functions
satisfy the normalization condition (A.7). Both spreading functions transition from strongly forward
peaked at low spatial frequencies (long gravity waves; the red curves) to curves with significant
propagation at right angles to the wind at high frequencies (capillary waves; the blue curves). The
cosine-2s curves are asymmetric in ±k and have at least a small amount of upwind propagation
at all frequencies (except at exactly upwind, ϕ = π). Not surprisingly, the real ocean is more
complicated than either of these models. In particular, observations of long-wave gravity waves tend
to show a bimodal spreading about the downwind direction, which transitions to a more isotropic,
unimodal spreading at shorter wavelengths (e.g. Heron et al., 2006). However, the simple models
of Eqs. (A.17) and (A.18) are adequate for the present purpose of illustrating surface-generation
techniques.

68
Figure A.5: Example spreading functions according to the ECKV model (left) and the cosine-2s
model (right) for a wind speed of 10 m s−1 . Small s values correspond to large spatial frequencies
k. Downwind is to the right, upwind is to the left.

69
References

J. H. G. M. Alves and M. L. Banner, 2003. Revisiting the Pierson-Moskowitz asymptotic limits for
fully developed wind waves. J. Phys. Oceanogr., 33:1301–1323.

R. R. Bracewell, 1986. The Fourier Transform and Its Applications, Second Edition, Revised.
McGraw-Hill.

J. W. Cooley and J. W. Tukey, 1965. An algorithm for the machine calculation of complex Fourier
series. Math. Comput., 19:297–301.

T. Elfouhaily, B. Chapron, K. Katsaros, and D. Vandemark, 1997. A unified directional spectrum

for long and short wind-driven waves. J. Geophys. Res., 102:15781–15796.

J. W. Goodman, 1996. Introduction to Fourier Optics, Second Edition. McGraw-Hill.

M. L. Heron, W. J. Skirving, and K. J. Michael, 2006. Short-wave ocean wave slope models for use
in remote sensing data analysis. IEEE Trans. Geosci. Rem. Sens., 44:1962–1973.

L. H. Holthuijsen, 2007. Waves in Oceanic and Coastal Waters. Cambridge University Press.

P. A. Huang, D. W. Wang, E. J. Walsh, W. B. Krabill, and R. N. Swift, 2000. Airborne measure-

ments of the wavenumber spectra of ocean surface waves. Part II: Directional distribution. J.
Phys. Oceanogr., 30:2768–2787.

S. Kay, J. Hedley, S. Lavender, and A. Nimmo-Smith, 2011. Light transfer at the ocean surface
modeled using high resolution sea surface realizations. Optics Express, 19:6493–6504.

M. S. Longuet-Higgins, D. E. Cartwright, and N. D. Smith, 1963. Observations of the directional

spectrum of sea waves using the motions of a flotation buoy. In Ocean Wave Spectra, pages
111–136. Prentice Hall.

C. D. Mobley, 2014. HydroPol Mathematical Documentation: Invariant Imbedding Theory for the
Vector Radiative Transfer Equation. Technical report, Sequoia Scientific, Inc., Bellevue, WA
98075, 2014. URL http://www.oceanopticsbook.info/view/references/publications.

C. D. Mobley, 1994. Light and Water: Radiative Transfer in Natural Waters. Academic
Press, 1994. URL www.oceanopticsbook.info/view/introduction/level_2/text_books_
relevant_to_ocean_optics.

70
W. J. Pierson and L. Moskowitz, 1964. A proposed spectral form for fully developed wind seas
based on the similarity theory of S. A. Kitaigorodskii. J. Geophys. Res., 69:5181–5190.

R. W. Preisendorfer, 1976. Hydrologic Optics. U.S. Department of Commerce, Pacific Marine

Environmental Laboratory. URL www.oceanopticsbook.info/view/introduction/level_2/
text_books_relevant_to_ocean_optics. In 6 volumes; Volume VI treats surfaces.

W. H. Press, S. A. Teukolksy, W. T Vetterling, and B. P. Flanney, 1992. Numerical Recipes in

Fortran: The Art of Scientific Computing, Second Edition. Cambridge Univ. Press.

A. Schuster, 1898. On the investigation of hidden periodicities with application to a supposed 26

day period of meteorological phenomena. Terrestrial Magnetism, III:13–41.

J. Tessendorf, 2004. Simulating Ocean Water. Technical report, SIGGRAPH Course Notes. URL
http://people.clemson.edu/~jtessen/reports.html.

71