July 31, 2017 16:58 ws-book9x6 10666-main page 1

Chapter 1
by [Link] on 11/12/25. Re-use and distribution is strictly not permitted, except for Open Access articles.

Introduction
Ocean Surface Waves Downloaded from [Link]

1.1 Waves in the Ocean and Their Significance

Atmosphere and ocean form a coupled system which continuously exchanges

heat, momentum and mass at the air-sea interface governing the dynamics
of weather, climate and global change. Due to energy flow from the atmo-
sphere to the ocean this interface presents an aerodynamically rough surface
that can comprise dynamic, unsteady, very high and steep surface waves,
which are the results of the various external forces of nature dictating what
types of waves can be induced in the ocean. A realistic description of the
physical processes occurring at the ocean-atmosphere interface is essential
for a reliable determination of the air-sea fluxes of momentum, sensible and
latent heat, CO2 and other trace gases, and aerosols. As the wave field is
involved in these exchange processes, wave models should be able to com-
pute not only the wave spectra and wave statistics, but also the processes
at the air-sea interface that govern the fluxes across the interface.
The most obvious cause of surface waves is an action of wind. The
ancient Greeks were well aware of interaction between the atmosphere and
sea surface. Aristotle (384–322 BC) realised that wind acting on the sea
surface plays a very important role in the development of waves. Pliny (AD
23–79) observed that oil poured upon waves calms them. From the time of
Aristotle to the Renaissance of science in the ‘Golden Age of Discovery’ in
the late fifteenth century, very little progress was made towards an under-
standing of the generation and growth of waves. It was not until well into
the nineteenth and twentieth centuries that more fundamental knowledge
of what causes waves and how they behave was accumulated.
At present, the study of surface wave mechanics are concentrated pre-
dominantly on three aspects. On the one hand, wave mechanics was

1
July 31, 2017 16:58 ws-book9x6 10666-main page 2

2 Introduction

treated as a pure mathematical problem, important also for numerical mod-

elling (Stokes, 1847, 1880; Lamb, 1932; Stoker, 1957; Phillips, 1966, 1977;
LeBlond and Mysak, 1978; Davidan et al., 1978; Mei, 1983; Komen et al.,
1994; Massel, 1989, 2007, 2013; Dingemans, 1997; Ochi, 1998; Kantha and
Clayson, 2000; Lavrenov, 2003; Pedlosky, 2003; Janssen, 2004; Mei et al.,
by [Link] on 11/12/25. Re-use and distribution is strictly not permitted, except for Open Access articles.

2005; Holthuijsen, 2007; Babanin, 2011 and others).

On the other, dynamics of surface waves were used to estimate the re-
sultant wave loads on offshore and coastal structures under extreme and
operating conditions. Textbooks and scientific papers provides many ex-
amples of the direct applications of wave mechanics to ocean and coastal
Ocean Surface Waves Downloaded from [Link]

engineering (Krylov et al., 1976, 1986; Sarpkaya and Isaacson, 1981; Dean
and Dalrymple, 1998; Goda, 1985, 2000; Mei et al., 2005 and others).
In recent years, an increasing attention has been directed to climate
change and the role played by the oceans. The interaction between the
ocean and atmosphere at the air-sea interface is critical to our understand-
ing of the earth’s climate as ocean surface forms a filter to the exchange
of heat, moisture, momentum and trace constituents. The energy trans-
fer from the atmosphere to the ocean enhances heat flux and mixed layer
during the circulation of the upper ocean, while the energy fluxes from the
ocean to the atmosphere affect gas and aerosol concentration, as well as
the atmospheric circulation, weather and climate. Oceans are not empty,
but support a large variety of organisms and plants. Ocean waters inter-
act, influence and constrain the ocean’s biological life. In particular, ocean
waters are a critical factor for enhancement of mixing in the upper ocean
layer and transferring of oxygen and nutrients into water column (Mann
and Lazier, 1996; Massel, 1999, 2007; Janssen, 2004; Babanin, 2011).
In general, five basic types of waves, present in the ocean can be dis-
tinguished: sound, capillary, gravity, internal and planetary waves. Sound
waves are due to water compressibility, which is, in fact, very small. Gravity
forces, acting on water particles displaced from equilibrium at the ocean
surface or at an internal geopotential surface in a stratified fluid, induce
gravity waves (surface or internal). The combination of the turbulent wind
and surface tension gives rise to short, high frequency capillary waves. On
the other hand, very slow, large-scale planetary or Rossby waves are in-
duced by the variation of the equilibrium potential vorticity due to changes
in depth or latitude. All of the above wave types can occur together, pro-
ducing more complicated patterns of oscillations.
Moreover, the gravitational attraction between the Earth, Moon and
Sun generate tides, which are very long oceanic waves with periods of the
July 31, 2017 16:58 ws-book9x6 10666-main page 3

1.1. Waves in the Ocean and Their Significance 3

order of 12 hours and more. Complexity of land contours on the Earth and
complicated bathymetry cause great difficulties in the prediction of tides at
some points on the Earth. However, the tides can still be predicted using
the harmonic analysis of a sufficiently long record of water level fluctuation.
Under specific conditions, long waves can also be generated in closed wa-
by [Link] on 11/12/25. Re-use and distribution is strictly not permitted, except for Open Access articles.

ter basins such as lakes, bays, estuaries and harbours, and they are known
as seiches. The most probable sources of seiches forcing mechanisms are
meteorological disturbances, impact of wind gusts on the water surface,
seismic oscillations of the earth during earthquakes and others. The fre-
quency of the free oscillations (natural oscillations) are determined by the
Ocean Surface Waves Downloaded from [Link]

shape of the basin and water depth. However, the meteorological conditions
may considerably change the height of predicted seiches and complicate the
prediction process.
The oscillations of water level known as the storm surges, are the dis-
turbances of sea level due to meteorological factors resulting from wind,
atmospheric pressure changes, moving storm systems and bottom configu-
ration. Storm surges may be expected to occur along coasts with relatively
shallow waters which are affected by the passage of storms.
Strong earthquakes with an epicentre in or near an ocean basin generate
large disturbances of the ocean surface known as tsunami. Tsunami are
impulsively driven water waves caused by the sudden displacement of a large
mass of water. They propagate radially from their origin on the ocean
surface and their impact is particularly severe on the ocean coasts, ranging
from erosion and over-washing of beaches, localised or regional floods to
the complete devastation of coastal villages and structures. Occasionally,
during special wave events appearing in nature, such as meteorites impact
on ocean surface or glaciers calving in polar regions, impulsive type waves
of the multi-frequency nature are generated.
In general, the frequency range associated with external forces is very
wide and ocean surface response occupies an extraordinarily broad range
of wave lengths and periods, from capillary waves, with periods of less than
a second, through wind-induced waves and swell with periods of the order
of a few seconds, to tidal oscillations with periods of the order of several
hours and storm surges with periods of the order of days. In Fig. 1.1 and in
Table 1.1, the schematic representation of energy contained in the surface
waves, and the physical mechanisms generating these waves, are listed. The
Figure gives some impression of the relative importance of various kinds
of surface oscillations, but does not necessarily reflect the actual energy
content.
July 31, 2017 16:58 ws-book9x6 10666-main page 4

4 Introduction
by [Link] on 11/12/25. Re-use and distribution is strictly not permitted, except for Open Access articles.
Ocean Surface Waves Downloaded from [Link]

Fig. 1.1 Schematic distribution of wave energy in frequencies.

Table 1.1 Waves, physical mechanisms, and periods.

Wave type Physical mechanism Periods

Capillary waves Surface tension < 10−1 s

Wind waves Wind shear, gravity < 15 s
Rogue waves Nonlinear mechanisms < 15 s
Swell Wind waves < 30 s
Surf beat Wave groups 1–5 min
Seiche Meteorological disturbances 2–40 min
Internal waves Density structure instability 10 min – 15 h
Tsunami Earthquake 10 min – 2 h
Tides Gravitational action of the moon 12–24 h
and sun, earth rotation
Storm surges Wind stresses and atmospheric 1–3 days
pressure variation

For offshore engineering operations, offshore structures and ship perfor-

mance, gravity waves are of the greatest importance, as the influence of
the loads due to gravity waves on engineering structures is the most criti-
cal and hostile. Marine structures must be designed to sustain the forces
and velocities induced by waves. An accurate assessment of the maximum
water particle velocities beneath the high wave crest is required for drag
July 31, 2017 16:58 ws-book9x6 10666-main page 5

1.1. Waves in the Ocean and Their Significance 5

force calculations, while the maximum accelerations beneath the steepest

section of the wave profile are critical for inertial force calculations. There-
fore, a thorough understanding of the interaction of waves with offshore
structures has now become a vital factor in the safe and economical de-
sign of such structures. The calculation procedures necessary to establish
by [Link] on 11/12/25. Re-use and distribution is strictly not permitted, except for Open Access articles.

the structural loading generally involve the few steps in which knowledge
on the surface waves is essential, namely establishing the wave climate
in the vicinity of the structure, estimating design wave conditions for the
structure, and selecting and applying a wave theory to determine the hy-
drodynamical loading on the structure. In particular, the exceptionally
Ocean Surface Waves Downloaded from [Link]

large waves, called rogue waves or freak waves, represent a danger to ships
and marine structures. A number of ship and offshore platform accidents
related to the action of rogue waves have been reported in the literature.
The role of waves in coastal marine environment can not be overes-
timated. Waves approaching the shoreline break and dissipate their en-
ergy on beaches. Storm and cyclone waves impose large forces on natural
coastal and man made structures. Long-shore currents transport sediments
and create areas of erosion and accumulation. The waves are very impor-
tant for transporting energy, sediment and chemical and biological matter.
Knowledge of wave motion and the sediment budget provides the key to the
proper selection of protecting structures and methods of shoreline stabilisa-
tion. Our present understanding of these processes is still not adequate to
develop effective management plans or engineering designs which are today
subjected to very stringent requirements to minimise their impact on the
environment. Being dependent on the quality of the environment, human
activities are recognised as drivers of significant pressures on marine and
coastal biodiversity, pollution, and nutrient emissions.
Other types of ocean waves which play a smaller role in ocean and
coastal engineering or have been treated in other books are not included
here. In particular, for a systematic review of planetary waves, tides and
storm surges, internal waves and coastal trapped Kelvin waves I refer the
reader to the books by Phillips (1977), LeBlond and Mysak (1978), Pedlosky
(1986), Monin (1988), Cushman-Roisin (1994), Miropolsky (2001), Dijkstra
(2005), Flór (2010), Sutherland (2010) and Massel (2015).
July 31, 2017 16:58 ws-book9x6 10666-main page 6

6 Introduction

1.2 Basic Assumptions on Seawater and Wave Motion

1.2.1 Continuous fluid and water particle concepts

Prior to introducing the basic concepts on fluids and fluid motion, we will
adopt the rectangular coordinate system O(x, y, z) or O(xi , z); i = 1, 2.
by [Link] on 11/12/25. Re-use and distribution is strictly not permitted, except for Open Access articles.

The origin of the system is at the mean sea surface. The axes x and y are
horizontal. The z axis is directed opposite to the force of gravity.
A fluid is composed of a large number of molecules which are in constant
motion and underlying collisions with each other. Therefore when studying
the mechanics of a fluid it is possible in principle to consider the motion
Ocean Surface Waves Downloaded from [Link]

of the fluid molecules themselves, as it done in kinetic theory or statistical

mechanics. However, the average quantities are usually sufficient for the
explanation of the observed macroscopic phenomena. Therefore, the dis-
crete molecular structure is replaced with a continuous distribution, known
as continuum. By continuous fluid in motion, we mean that the velocity u
is everywhere finite and continuous while its space derivatives of the first
order are finite (but not necessarily continuous). Thus, any closed surface
S which moves with the fluid, permanently and completely separates the
fluid matter inside S from that outside. The fluid can be treated as contin-
uous at the scale A (obstacle dimension, a pore diameter, a turning radius
et.) when the flow past an obstacle of the dimension A is much larger than
the average free path of the molecule l0 (for water l0 ≈ 3 · 10−10 m). The
mean free path specifies the average distance that molecule travels before it
demonstrates its presence, temperature, or momentum to other molecules.
The rate (l0 /A) is known as the Knudsen number (Kn). If
l0
Kn = < 0.01, (1.1)
A
the fluid can be treated as continuous (Kundu et al., 2016). For most
terrestial situations, above requirement is not a great restriction. However,
for flow over small objects or in very narrow passages, the molecular-kinetic
approach is necessary.
With regard to continuous fluid, we can define a fluid particle as con-
sisting of the fluid contained within an infinitesimal volume, that is to say,
a volume whose size may be considered so small that for the particular pur-
pose in hand its linear dimensions are negligible. We can then treat a fluid
particle as a geometrical point. According to continuum hypothesis, fluid
properties like temperature, density, or velocity are defined at every point
in space.
August 14, 2017 10:22 ws-book9x6 10666-main page 7

1.2. Basic Assumptions on Seawater and Wave Motion 7

1.2.2 Properties of seawater and its motion

The equations of motion for the fluid particle depends on the physical
properties of the fluid and motion itself. In order to render the subject
amenable to exact mathematical treatment, we examine some physical and
chemical properties of seawater and simplifying assumptions on the medium
by [Link] on 11/12/25. Re-use and distribution is strictly not permitted, except for Open Access articles.

and its motion, i.e.:

(a) seawater temperature, salinity and density
Seawater is not a pure water. It contains salts, dissolved organic substances,
and mineral and organic suspended matter. Water has many unique chem-
Ocean Surface Waves Downloaded from [Link]

ical and physical properties. Unlike most substances which contract when
frozen, water expands, allowing ice to float on the surface of water basins.
The principal seawater state parameter is density, or mass per unit volume.
The density of water controls many fundamental processes in the ocean, e.g.
the hydrostatic equilibrium or motion of water particles and propagation
of internal waves, turbulence and mixing in water column, heat transfer,
concentration of plankton and sediment transport, locomotion of marine
species and many others. Even small alterations in the density of seawater
result in great changes in water flow and its thermal and chemical status.
Density of water is normally symbolised by %w and in the SI system it is
expressed in kilogram per cubic meter (kg m−3 ). For example, the density
of pure water is approximately 1000 kg m−3 being 770 to 890 times that of
air, %a , at sea level. Unlike measuring salinity, temperature and pressure,
there are no practical methods of measuring the density of seawater in situ.
The density of seawater is usually determined through its dependence on
pressure, p, temperature, T , and salinity, S. Pressure, p, has an insignif-
icant effect on the density of water for most applications, unless one is
dealing with water at great depths within the ocean.1 The density is much
more dependent on temperature, generally decreasing as the temperature
increases. However, the density of freshwater is not a monotonic function
of temperature, and water density reaches a maximum at temperature of
3.98◦ C under normal atmospheric pressure.2 This anomalous dependence
is a result of the specific molecular structure of water (Dera, 1992).

1 For a depth of 1 km the density of water is only 0.5% greater than at the surface,

and at the deepest point in the earth’s oceans, the Marianas Trench, water is about
6% denser than water at the surface, assuming that the temperature and salinity are
constant (Denny, 1993).
2 Normal atmospheric pressure at sea level, called the normal atmosphere, is defined to

be 1.01325 · 105 Pa at 15◦ C (1 Pa = 1 N m−2 ).

July 31, 2017 16:58 ws-book9x6 10666-main page 8

8 Introduction

The mean temperature, T , of the surface water of oceans is about 17.8◦ C

and it is 3.6◦ C higher than the mean air temperature on the Earth’s surface.
In the Pacific Ocean, the mean water temperature is equal to 19.37◦ C, in
the Atlantic Ocean – to 17.58◦ C and in the Indian Ocean – to 17.85◦ C, while
at the Arctic Ocean depth it is about (1.0–1.5)◦ C.
by [Link] on 11/12/25. Re-use and distribution is strictly not permitted, except for Open Access articles.

The relatively constant composition of the main constituents of sea salt

has made it possible to introduce a single parameter defining the salt con-
centration in seawater. This is known as salinity (S). Away from coasts, the
salinity varies from 32 to 37 PSU (Practical Salinity Units) and the mean
salinity of the ocean waters (excluding Arctic waters) is equal to 34.71 PSU.
Ocean Surface Waves Downloaded from [Link]

There is a weak correlation between salinity and the latitude. In general,

cold waters are less saline than warmer waters. In closed or semi-closed
seas, such as the Baltic Sea, salinity is equal to 7–8 PSU, being, however,
dependent on the rivers discharge, and the difference between precipitation
and evaporation.
Practical Salinity, S, is today a globally accepted standard for reporting
the salinity of seawater and has been in use since 1978. The Practical Salin-
ity Scale, PSS–78, is based on a carefully determined salinity-conductivity
relationship. According to the formulation of PSS–78, the Practical Salin-
ity, S, of a sample of seawater, is defined in terms of the ratio K15 of
the electrical conductivity of seawater at the temperature of 15◦ C and the
pressure of one standard atmosphere, to that of a potassium chloride (KCl)
solution, in which the mass fraction is 32.4356 · 10−3 , at the same tem-
perature and pressure. A K15 value exactly equal to 1 corresponds, by
definition, to a Practical Salinity exactly equal to 35 (UNESCO, 1981).
Except the Practical Salinity Scale, the so called Absolute Salinity SA
of seawater has come under discussion. The Absolute Salinity SA is defined
as the ratio of mass of dissolved material in seawater to the total mass
of seawater. In principle, the dissolved material includes not only those
compounds in water conducting electricity. Such components are usually
present in many local sea basins. However, the SA quantity cannot be
measured directly. Millero et al. (2008) proposed a new salinity variable,
called the Reference-Composition Salinity, or Reference Salinity, to provide
a means to more accurately estimate the Absolute Salinity SA of a local
seawater sample.
Density of seawater %w (S, T, p) for temperature T (0 < T < 40◦ C),
salinity S (0 < S < 42 PSU) and pressure p can be determined from the
July 31, 2017 16:58 ws-book9x6 10666-main page 9

1.2. Basic Assumptions on Seawater and Wave Motion 9

International Equation of State for Seawater (UNESCO, 1981; Dera, 1992)

%w (S, T, 0)
%w = p , (1.2)
1−
K(S, T, p)
by [Link] on 11/12/25. Re-use and distribution is strictly not permitted, except for Open Access articles.

where K(T, S, p) is module of seawater compressibility. The expression

for %w (S, T, p) is rather complicated and it will not be given here.3 For
illustration, some values of density %w (S, T, 0) are listed below:
%w (0, 10◦ , 0) = 999.702 kg m−3 ;
%w (35, 10◦ , 0) = 1026.952 kg m−3 ;
Ocean Surface Waves Downloaded from [Link]

%w (0, 20◦ , 0) = 998.206 kg m−3 ;

%w (35, 20◦ , 0) = 1024.763 kg m−3 .
In place of the water density, the specific weight will be sometimes used in
this book. This quantity, usually symbolised by γ, is equal to the product
of water density %w and gravitational acceleration, g, so γ = %w g. The
specific weight is expressed in newtons per cubic metre (N m−3 ).
The density of seawater and its distribution in the water column deter-
mines the hydrostatic stability of water masses and influences the turbu-
lence and internal waves and sound propagation. However, the influence of
density on surface waves is negligible, except perhaps for the stage of wave
generation under wind action (see Chapter 2).
(b) water is an inviscid and incompressible fluid
Viscosity is a property which is a measure of a fluid’s resistance to ‘deforma-
tion’ during motion. Within a fluid, momentum of rapidly moving particle
is exchanged with the momentum of relatively slower particles. Those ex-
changes produce a shearing stress. In general we can write for shear stress τ
du
τ =µ , (1.3)
dn
where τ is tangential stress, u = (u, v, w) is fluid velocity vector, n denotes
vector normal to vector u, and µ is the coefficient of dynamic molecular
viscosity of the fluid. It has units of (Ns m−2 ). A fluid which shows the
direct proportionality between the applied shear stress and the resulting
rate of deformation is called a Newtonian fluid. Many biological materials,
such as blood, synovial fluid, mucus of various consistencies, can not be
treated as Newtonian fluids. Some of them have a memory of previous
3 The algorithm of determination of density for given salinity, S, temperature, T , and

pressure, p, according to UNESCO Formula is given in UNESCO, 1981 or Massel, 2015.

July 31, 2017 16:58 ws-book9x6 10666-main page 10

10 Introduction

shape and elasticity. A description of such fluids is a subject of another

branch of fluid mechanics, so-called rheology.
As force on an element of fluid varies like τ but the mass of that element
varies like %w , the acceleration and hence the velocity field is determined
by the ratio τ /%w , known as the coefficient of kinematic viscosity, ν
by [Link] on 11/12/25. Re-use and distribution is strictly not permitted, except for Open Access articles.

µ
ν= , (1.4)
%w
where ν has units of m2 s−1 . For example, the coefficient of kinematic
viscosity for seawater of salinity S = 35 PSU and of temperature T = 20◦ C
becomes 1.064 · 10−6 m2 s−1 . It is interesting to note that there is almost
Ocean Surface Waves Downloaded from [Link]

no liquid with viscosity lower than that of water, and that air viscosity
νair ≈ 15 νwater .
Both coefficients, µ and ν, are physical properties of fluid, independent
of fluid motion. When a fluid is completely at rest or is moving very slowly,
the rate of diffusion of momentum is essentially determined by molecular
motion. Eddies are small and the velocity shear is great, and the internal
resistance of the water due to molecular viscosity smooths out the gradients
in velocity. This smoothing of the flow is the way the energy in turbulence
is finally converted to heat and dissipated.
However, as soon as the fluid is stirred and the elements of fluid are much
greater than molecules mix, the coefficient µ increases considerably. In the
oceans and atmosphere, eddies and turbulent motions in the flow can be so
effective in moving particles among themselves that the effects of molecular
diffusion are overwhelmed. For such situations, the coefficient of kinematic
viscosity, ν, is replaced by the coefficient of turbulent viscosity, which is
several hundred to many thousand times larger. By analogy to molecular
exchange, the turbulent viscosity coefficients result from the hypothesis that
the turbulent momentum flux is proportional to the averaged turbulent flow
velocity. Here we only note that the estimates of the coefficient of turbulent
viscosity in the ocean vary enormously, from 10 m2 s−1 to 104 m2 s−1 in the
horizontal plane, and from 10−4 m2 s−1 to 10−1 m2 s−1 in the vertical plane.
The ratio Re = uL/ν (the Reynolds number), in which L is the differen-
tial length scale and u is the flow velocity, represents the relative magnitudes
of the inertial and viscous terms. In many oceanic motions, the Reynolds
number is very large. Thus, the viscous influence is often quite negligible
over most of the field of motion. The viscous forces are important only
in narrow regions of the flow, where the local inertial and viscous forces
are comparable. In the ocean, the interfacial layer between the air and
the water, as well as the bottom boundary layer are such regions. The
July 31, 2017 16:58 ws-book9x6 10666-main page 11

1.2. Basic Assumptions on Seawater and Wave Motion 11

1/2
thickness δ of the surface boundary layer is of the order δ ≈ (2ν/ω) ,
where ω is a wave frequency, and ν is kinematic coefficient of viscosity. For
typical ocean wave frequencies, the thickness δ ∼ 0.001 m. For the bound-
ary layer near the natural sea bottom, the eddy viscosity is much higher.
The thickness of the boundary layer is then about 0.1 m, which is still quite
by [Link] on 11/12/25. Re-use and distribution is strictly not permitted, except for Open Access articles.

small. Therefore, the boundary layer regions are but a very small fraction
of a fluid volume, and the influence of the molecular viscosity on the wave
motion can be neglected.
The compressibility of water is very small and the Young’s modulus is
of order E ≈ 3.05 · 108 N m−2 (Dera, 1992). As the typical velocity of
Ocean Surface Waves Downloaded from [Link]

seawater is much smaller than the speed of sound, the very small water
compressibility has no influence on water motion.
(c) conservation laws
The basic principles in fluid mechanics are the conservation laws for mass,
momentum and energy. In the following we present them in differential
form, applicable at a point or to a fluid particle. Therefore, conservation
of mass takes the form of equation
∂ ∂ ∂ ∂
(%w ) + (%w u) + (%w v) + (%w w) = 0. (1.5)
∂t ∂x ∂y ∂z
For constant density throughout the flow field, and for incompressible flow
from Eq. (1.5) we have
∂u ∂v ∂w
div u = ∇ · u = + + = 0, (1.6)
∂x ∂y ∂z
where
∂ ∂ ∂
∇= i+ j+ k. (1.7)
∂x ∂y ∂z
The equations of motion for fluid provide the relationship between the
applied forces on a unit volume of fluid and the change of momentum.
The basic equation for derivation of the equation of motion is Newton’s
second law.
July 31, 2017 16:58 ws-book9x6 10666-main page 12

12 Introduction

For incompressible fluid the equation of motion are called the Euler
equations4
  
∂u ∂u ∂u ∂u ∂p
%w +u +v +w =−


∂t ∂x ∂y ∂z ∂x






by [Link] on 11/12/25. Re-use and distribution is strictly not permitted, except for Open Access articles.

  
∂v ∂v ∂v ∂v ∂p 
%w +u +v +w =− , (1.8)
∂t ∂x ∂y ∂z ∂y 



  
∂w ∂w ∂w ∂w ∂(p + %w gz) 

%w +u +v +w =−



∂t ∂x ∂y ∂z ∂z
or in vector notation 
Ocean Surface Waves Downloaded from [Link]


∂u 1 2 p
+∇ |u| + gz + + (curlu) × u = 0, (1.9)
∂t 2 %w
in which
rotu = curl u = ∇ × u
     
∂w ∂v ∂u ∂w ∂v ∂u
= − i+ − j+ − k, (1.10)
∂y ∂z ∂z ∂x ∂x ∂y
where u, v and w are the velocity components on x, y and z axes, respec-
tively. The vector rot u is known as the vorticity. When the vorticity is
different from zero, the motion is defined as rotational.
Some flows in the ocean are rotational due to presence of viscosity fric-
tion and the Earth’s rotation; they also reveal a lot of diffusion. Extension
of the Euler equations to viscous fluids is given in the form of the Navier–
Stokes equations
  
∂u ∂u ∂u ∂u ∂p 2
%w +u +v +w =− + µ∇ u 

∂t ∂x ∂y ∂z ∂x






  
∂v ∂v ∂v ∂v ∂p 2

%w +u +v +w =− + µ∇ v . (1.11)
∂t ∂x ∂y ∂z ∂y 



  
∂w ∂w ∂w ∂w ∂(p + %w gz)


2 
%w +u +v +w =− + µ∇ w  
∂t ∂x ∂y ∂z ∂z
In many oceanic motions the influence of the viscous terms is quite negli-
gible. In such case, the Lagrange theorem indicates that if, at some initial
instant, the vorticity vanishes everywhere in the field of flow, the motion is
irrotational. Therefore from Eq. (1.10) is rot u = 0 or
∂w ∂v ∂u ∂w ∂v ∂u
= , = , = . (1.12)
∂y ∂z ∂z ∂x ∂x ∂y
4 Leonhard Euler (1707–1783) was a Swiss mathematician and physicist who made very

important discoveries in many branches of mathematics, mechanics and fluid dynamics.

July 31, 2017 16:58 ws-book9x6 10666-main page 13

1.3. Fundamentals of Description of Regular Waves 13

The above conditions can be satisfied when the velocity u is represented as

the gradient of a scalar function, the velocity potential Φ
u = ∇Φ. (1.13)
Then, in virtue of the continuity equation (1.6), the potential Φ obeys the
Laplace equation
by [Link] on 11/12/25. Re-use and distribution is strictly not permitted, except for Open Access articles.

∂2Φ ∂2Φ ∂2Φ

∇2 Φ = + + = 0. (1.14)
∂x2 ∂y 2 ∂z 2
It should be noted that velocity potential Φ(x, z) and stream function
ψ(x, z) constitute in a vertical plane (x, z) the analytical function known
as complex velocity potential w(x, z) = Φ(x, z) + iψ(x, z), satisfying the
Ocean Surface Waves Downloaded from [Link]

Cauchy–Riemann relations
∂Φ(x, z) ∂ψ(x, z) ∂Φ(x, z) ∂ψ(x, z)
= , =− . (1.15)
∂x ∂z ∂z ∂x

1.3 Fundamentals of Description of Regular Waves

1.3.1 Linear description of waves

The deterministic description is a natural starting point for the development
of random wave models. Thus, let us consider a wave train observed at some
point P (x, y). The surface displacement of a wave travelling at an angle θ
to the x axis may be presented as
ζ(x, y, t) = a cos [k(x cos θ + y sin θ) − ωt + ϕ] , (1.16)
in which a is the wave amplitude, ϕ is a phase shift and k is a wavenumber
related to wave frequency ω by the classical dispersion relation
ω 2 = gk tanh(kh). (1.17)
Solution of Eq. (1.17) in terms of the wavenumber k is given by Eq. (3.111).
The simplest and the most natural representation of the confused sea
surface would be the linear superposition of many harmonics travelling in
various directions. A simple illustration of such a superposition is given in
Fig. 1.2, where 13 elementary components sum to form a final wave surface
displacement. Thus, using Eq. (1.16), the surface displacement may be
written as
n=N
X
ζ(x, y, t) = an cos [kn (x cos θn + y sin θn ) − ωn t + ϕn ] . (1.18)
n=1
The direction θn and phase ϕn cover the range −π, π, and wave amplitudes
and frequencies are from the range 0 ≤ an ≤ ∞ and 0 ≤ ωn ≤ ∞.
July 31, 2017 16:58 ws-book9x6 10666-main page 14

14 Introduction
by [Link] on 11/12/25. Re-use and distribution is strictly not permitted, except for Open Access articles.
Ocean Surface Waves Downloaded from [Link]

Fig. 1.2 Superposition of 13 spectral components and resulting spectrum.

The integral version of Eq. (1.18) would be

Z ∞Z π
ζ(x, y, t) = 2 a(ω, θ) cos [k(x cos θ + y sin θ) − ωt + ϕ(ω, θ)] dωdθ,
0 −π
(1.19)
July 31, 2017 16:58 ws-book9x6 10666-main page 15

1.3. Fundamentals of Description of Regular Waves 15

under the assumption that (Yaglom, 1962; Borgman, 1972)

N
1X
a(ω, θ) = an δ(ω − ωn )δ(θ − θn ), (1.20)
2 n=1

where δ(x) denotes Dirac’ delta function5 (Lighthill, 1975). Using Euler’s
by [Link] on 11/12/25. Re-use and distribution is strictly not permitted, except for Open Access articles.

identity
1
cos α = [exp(iα) + exp(−iα)] , (1.21)
2
we rewrite Eq. (1.19) in the form
Z ∞Z π
Ocean Surface Waves Downloaded from [Link]

ζ(x, y, t) = < a(ω, θ) exp [iϕ(ω, θ)]

−∞ −π

× exp [ik(x cos θ + y sin θ) − ωt] dωdθ. (1.22)

Let us define the following cumulative function
Z ω Z θ
A(ω, θ) = a(ω 0 , θ0 ) exp [iϕ(ω 0 , θ0 )] dω 0 dθ0 . (1.23)
−∞ −π

Then, Eq. (1.22) becomes

Z ∞ Z π
ζ(x, y, t) = exp [i (k · x − ωt)] dA(ω, θ), (1.24)
−∞ −π

where dA(ω, θ) is a differential element in the two-dimensional space (ω, θ),

k is a two-dimensional vector (k cos θ, k sin θ) and x is the vector (x, y).
For later convenience it will be useful to express Eq. (1.24) as a Fourier
transform in three dimensions, (x, y, t). Thus, we define
u = k cos θ, v = k sin θ, w = −ω (1.25)
and
Z u Z v Z w
Â(u, v, w) = â(u0 , v 0 , w0 )dw0 dv 0 du0 , (1.26)
−∞ −∞ −∞

in which
h i
1/2
â(u, v, w) = A(ω, θ) δ |ω| − (gk tanh(kh)) . (1.27)
5 The Dirac’s delta has two useful properties

∞ when x = x0
δ(x − x0 ) =
0 when x 6= x0
and Z ∞
f (y)δ(x − y)dy = f (x).
−∞
July 31, 2017 16:58 ws-book9x6 10666-main page 16

16 Introduction

Using Eq. (1.25), the Eq. (1.24) may be rewritten as

Z ∞Z ∞Z ∞
ζ(x, y, t) = exp [i(ux + vy + wt)] dÂ(u, v, w). (1.28)
−∞ −∞ −∞
Due to the restriction Eq. (1.27) imposed on frequency ω and wavenumber
k, a triple integral in Eq. (1.28) is in fact a double integral only.
by [Link] on 11/12/25. Re-use and distribution is strictly not permitted, except for Open Access articles.

The advantage of the integral representation (1.28) is that it makes

application of the technique of Fourier integrals possible. This means that
for given measurements of ζ(x, y, t) over the (x, y, t) space, the function
â(u, v, w) can be computed in principle by an inverse Fourier transform. It
can be seen from Eq. (1.27) that â(u, v, w) has nonzero values only when
Ocean Surface Waves Downloaded from [Link]

(u, v, w) satisfy the dispersion relation (1.17).

1.3.2 Nonlinear description of waves

[Link] Second-order Stokes waves
Linear description of ocean waves given by Eqs. (1.18) or (1.28) properly
reflects only ocean waves with a small steepness. For higher steepness of
regular waves, the second-order solution was first developed by (Stokes,
1847, 1980). However, for waves being superposition of components of var-
ious steepnesses and phases, an extension of the Stokes solution, which can
also capture the effects of water depth, and directional spreading, is re-
quired. Such second-order solution for the infinite water depth was given
by Longuet-Higgins (1963), and later by Sharma and Dean (1979) for ar-
bitrary water depth. Their description was subsequently reproduced by
Forristall (2000) as follows
ζ(x, t) = ζ (1) (x, t) + ζ (2) (x, t), (1.29)
in which
n=N
X
ζ (1) (x, t) = an cos(kn · x − ωn t + ϕn ) (1.30)
n=1
and
N N
1 XX h i
ζ (2) (x, t) = ai aj K (−) cos(ψi − ψj ) + K (+) cos(ψi + ψj ) ,
4 i=1 j=1
(1.31)
where
(−)
Dij − (ki · kj + Ri Rj )
K (−) = p + (Ri + Rj ), (1.32)
Ri Rj
July 31, 2017 16:58 ws-book9x6 10666-main page 17

1.3. Fundamentals of Description of Regular Waves 17

(+)
Dij − (ki · kj − Ri Rj )
K (+) = p + (Ri + Rj ), (1.33)
Ri Rj

√ p
p √ 
(−) Ri −
Rj Rj (k2i − Ri2 ) − Ri (k2j − Rj2 )
by [Link] on 11/12/25. Re-use and distribution is strictly not permitted, except for Open Access articles.

Dij = √ p
2 (−) (−)

Ri − Rj − kij tanh(kij h)
√ p
2
2 Ri − Rj (ki · kj + Ri Rj )
+ √ p
2 (−) (−)
, (1.34)
Ri − Rj − kij tanh(kij h)
Ocean Surface Waves Downloaded from [Link]

√
√ 
Ri (k2j − Rj2 ) + Rj (k2i − Ri2 )
p p
(+) Ri + Rj
Dij = √ p
2 (+) (+)
Ri + Rj − kij tanh(kij h)
√ p
2
2 Ri + Rj (ki · kj − Ri Rj )
+ √ p
2 (+) (+)
, (1.35)
Ri + Rj − kij tanh(kij h)

(−) (+)
kij = |ki − kj |, kij = |ki + kj |, (1.36)

Ri = |ki | tanh(|ki |h) = ωi2 /g, ψi = ki · x − ωi t + ϕi . (1.37)

For infinite water depth, Eq. (1.31) reduces to the result of Longuet-Higgins
(1963), except that the latter equation is missing a factor of 1/2.
The double sum in Eq. (1.31) contains all the second-order bound waves
that are generated by interactions between the free waves.6 The positive
interaction terms (ψi + ψj ) produce the sharpening of the crests and flat-
tening of the troughs. The negative interaction terms (ψi − ψj ) yield the
set-down of the water level under wave groups. It is clearly seen in Fig. 1.3
that the bound waves do not satisfy the linear dispersion relation, as
s
h
(ωi ± ωj ) 6= (ki ± kj ) tanh[(ki ± kj )h]. (1.38)
g

At present, the second-order wave description is widely used in the off-

shore industry (Birknes and Bitner-Gregersen, 2003; Bitner-Gregersen and
Hagen, 2004).
6 It should be noted that the wave energy spectrum does not change in time as the

wave field evolves. Therefore in the second-order solution, the nonlinear dynamics of
free waves and the modulational instability are not included.
July 31, 2017 16:58 ws-book9x6 10666-main page 18

18 Introduction
by [Link] on 11/12/25. Re-use and distribution is strictly not permitted, except for Open Access articles.
Ocean Surface Waves Downloaded from [Link]

Fig. 1.3 Dispersion relation for free and bound waves.

[Link] Higher-order Stokes waves

For single regular wave profile, the Stokes type solution to a higher order can
be found. If we assume the regular, plane wave propagating over constant
water depth in x-direction, we obtain the following fifth-order wave surface
profile (Fenton, 1985)

1
ζ(x, t) = ε cos ψ + ε2 B22 cos 2ψ + ε3 B31 (cos ψ − cos 3ψ)
k
+ε4 (B42 cos 2ψ + B44 cos 4ψ) + ε5 [−(B53 + B55 ) cos ψ
+ B53 cos 3ψ + B55 cos 5ψ]} + 0(ε6 ), (1.39)

where ε = k (H/2) and ψ = kx − ωt. The set of coefficients Bij is given in

Fenton (1985).
The fifth-order Stokes waves are often used for assessment of the max-
imum surface displacement and maximum water velocity and acceleration
beneath the high wave crest. The velocity and acceleration are required for
drag force and inertia force calculations. Therefore, for completeness we
provide the formulae for the horizontal and vertical velocity components,
July 31, 2017 16:58 ws-book9x6 10666-main page 19

1.3. Fundamentals of Description of Regular Waves 19

as well as for the Eulerial and Lagrangian accelerations (Fenton, 1985):

– horizontal velocity component
 g 1/2 X5 X i
u(x, z, t) = C0 Aij εi jk cosh[jk(z + h)] cos(jψ) (1.40)
k3 i=1 j=1
by [Link] on 11/12/25. Re-use and distribution is strictly not permitted, except for Open Access articles.

– vertical velocity component

 g 1/2 X5 X i
w(x, z, t) = C0 Aij εi jk sinh[jk(z + h)] sin(jψ), (1.41)
k3 i=1 j=1

where
Ocean Surface Waves Downloaded from [Link]

ψ = kx − ωt. (1.42)
The wavenumber k is given by the formula
2π
− √ + C0 (kh) + ε2 C2 (kh) + ε4 C4 (kh) = 0. (1.43)
T gk
The coefficients Aij and Ci can be found in Fenton (1985).
(E)
Using the above formulae we obtain the Eulerian az and Lagrangian
(L)
az accelerations as follows
∂ 2 ζ(x, t) ∂w(x, ζ, t)
a(E)
z = 2
= (1.44)
∂t ∂t
and
Dw ∂w ∂w ∂w ∂w ∂u ∂u
a(L)
z (x, ζ, t) = = +u +w = +u −w at z = ζ.
Dt ∂t ∂x ∂z ∂t ∂z ∂x
(1.45)
For illustration, in Fig. 1.4, the comparison of the surface displacements
for linear, second-order and fifth-order approximation is given with the
following wave parameters used in calculations: water depth h = 20 m,
wave length L = 250 m, wave height H = 6 m, H/h = 0.3 and L/h = 12.5.

[Link] Nonlinear Schrödinger (NLS) equations and modulational

instability
In the Stokes approximations, even as high as the fifth-order, the evolu-
tion of waves in time or in space is not included and the energy spectrum
remains constant. However, experiments in both natural and laboratory
conditions demonstrate clearly that wave train in general is unstable and
undergoes modulations due to nonlinear character of wave motion. When
July 31, 2017 16:58 ws-book9x6 10666-main page 20

20 Introduction
by [Link] on 11/12/25. Re-use and distribution is strictly not permitted, except for Open Access articles.
Ocean Surface Waves Downloaded from [Link]

Fig. 1.4 Comparison of linear and higher order approximations of surface displacement.

the ocean surface waves are regarded as narrow spectrum waves, the non-
linear evolution of the wave train in time is studied by using the cubic
Nonlinear Schrödinger (NLS) equation (Peregrine, 1983a; Kharif and Peli-
novsky, 2003; Mei et al., 2005; Onorato et al., 2006). The NLS equation
describing the complex envelope of sea elevations for arbitrary water depth
(in a frame of reference moving with the group velocity) is
∂A ωp ∂ 2 A ωp kp2
+ iσ 2 + iβ |A|2 A = 0, (1.46)
∂t 8kp ∂x2 2
where kp is the carrier wavenumber, ωp is the respective angular frequency
and σ and β are functions of the nondimensional wavenumber kp h
cosh(2kp h)
σ = 2 − ν12 + 8(kp h)2 , (1.47)
sinh2 (2kp h)

8 + cosh(4kp h) − 2 tanh2 (kp h)

β=
8 sinh4 (kp h)
[2 cosh2 (kp h) + 0.5ν1 ]2
− (1.48)
sinh2 (2kp h)(kp h)/ tanh(kp h) − ν12 /4
July 31, 2017 16:58 ws-book9x6 10666-main page 21

1.3. Fundamentals of Description of Regular Waves 21

and
kp h
ν1 = 1 + 2 . (1.49)
sinh(2kp h)
When kp h → ∞, coefficients σ and β tend to 1.
The complex amplitude A is a slowly varying function of x and t with
by [Link] on 11/12/25. Re-use and distribution is strictly not permitted, except for Open Access articles.

a third-order accuracy in the wave slope ka, and the surface elevation ζ(x, t)
is given by
1
ζ(x, t) = (A(x, t) exp[i(kp x − ωp t)] + c.c.) , (1.50)
2
where c.c. denotes the complex conjugate.
Ocean Surface Waves Downloaded from [Link]

The modulational instability leads to decomposition of an initially ho-

mogeneous wave train into a system of envelope quasi-solitons, and the
spatial distribution of wave energy displays significant intermittency (see
Fig. 1.5).
The wave steepness and wave dispersion play a critical role for an evo-
lution and stability of the water waves (Benjamin and Feir, 1967). The
popular measure of the relative importance of nonlinearity and wave disper-
sion is so called Benjamin–Feir index (BFI), introduced initially by Janssen
(2003). To develop the Benjamin–Feir index we define following Onorato
et al. (2001, 2006) the nondimensional quantities: Ã = A/a0 , x̃ = x∆K,
t̃ = t(∆K/kp )2 σωp /8, where ∆K is a typical spectral band-width and a0
is a typical wave amplitude.

Fig. 1.5 Formation of a highly energetic wave group in a slowly modulated wave train.
The time (0 on the left and 320 on the right) is normalised by the fundamental wave
period and the coordinate is in radius of the fundamental mode (adapted from Kharif
and Pelinovsky, 2003).
July 31, 2017 16:58 ws-book9x6 10666-main page 22

22 Introduction

After substituting Ã, x̃ and t̃ into Eq. (1.46) we obtain

2
∂ 2 Ã

∂ Ã 2kp a0 β
+i + i |Ã|2 Ã = 0. (1.51)
∂ t̃ ∂ x̃2 ∆K/kp σ
The Benjamin–Feir index BFI is defined as the square root of the coefficient
by [Link] on 11/12/25. Re-use and distribution is strictly not permitted, except for Open Access articles.

that multiples the nonlinear term in above equation, i.e.

r
2kp a0 |β|
BFI = . (1.52)
∆K/kp σ
p
The term |β|/σ represents the influence of the water depth. When wa-
ter depth increases, this term goes to one and for shallow water depth it
Ocean Surface Waves Downloaded from [Link]

becomes zero. For values of kp a0 < 1.36, the coefficient in front of the non-
linear term in the NLS equation changes sign and the equation becomes
stable with respect to side perturbations.
As the BFI increases the nonlinearity increases, and the maximum in-
stability occurs when ∆K/kp = 2kp a0 . Therefore for wave field sufficiently
steep, narrow banded, and unidirectional, the random waves become un-
stable when BFI = O(∞). If the time series (or frequency spectra) are
available, the following definition of BFI is used
r
(kp a0 ) |β|
BFI = , (1.53)
∆ω/ωp σ
which can be rewritten in terms of the significant wave height Hs and peak
wavenumber kp (Bitner-Gregersen and Toffoli, 2014)
kp Hs
BFI = . (1.54)
2∆ω/ωp
The computation methods of the BFI for a given time series or wave spec-
trum were considered by Serio et al. (2005). For more information on
the Benjamin–Feir instability, the reader should consult Benjamin (1967),
Krasitskii (1994) and Janssen (2003).
Dysthe (1979) extended the cubic NLS equation to fourth order in waves
steepness and bandwidth which can be solved more easily by a pseudospec-
tral numerical method (Lo and Mei, 1985). Beside the Dysthe extension,
there are numerous other variants of NLS-type equations for special condi-
tions (Trulsen and Dysthe, 1997; Onorato et al., 2000, 2001, 2006; Kharif
and Pelinovsky, 2003; Kurkin and Pelinovsky, 2003; Toffoli et al., 2007,
2008a, 2010b, 2011, 2016).
In the laboratory practice, wave data are usually collected from sensors
located at fixed positions in the wave basin. Therefore for comparison
July 31, 2017 16:58 ws-book9x6 10666-main page 23

1.3. Fundamentals of Description of Regular Waves 23

of such data with numerical model predictions, the governing equations

have to be transferred to the spatial form (Lo and Mei, 1985). The exact
expressions for the temporal and spatial modified NLS equations (MNLS)
for complex amplitude A(x, y, t) of the directional wave field take the form
(Socquet-Juglard et al., 2005; Toffoli et al., 2010b; Zhang et al., 2016)
by [Link] on 11/12/25. Re-use and distribution is strictly not permitted, except for Open Access articles.

∂A 1 ∂A i ∂2A i ∂2A i
+ + − + |A|2 A
∂t 2 ∂x 8 ∂x2 4 ∂y 2 2

1 ∂3A 3 ∂3A 3 2 ∂A 1 2 ∂A∗

− + + |A| + A (1.55)
16 ∂x3 8 ∂x∂y 2 2 ∂x 4 ∂x
Ocean Surface Waves Downloaded from [Link]

∂ φ̄
+ iA = 0, z=0
∂x
and
∂A ∂A ∂2A i ∂2A 2 ∂3A
+2 +i 2 − + i|A| A −
∂x ∂t ∂t 2 ∂y 2 ∂t∂y 2
(1.56)
∗
∂A ∂A ∂ φ̄
− 8|A|2 − 2A2 − 4iA = 0, z = 0,
∂t ∂t ∂t
∗
where A is the complex conjugate of wave amplitude A.
The main shortcomings of NLS equation are the requirements of nar-
row wave spectrum and slow modulation of waves, which limits its practical
applicability. A more general alternatives to NLS equations are so-called
Zakharov equations (Zakharov, 1968), which assume weak nonlinearity but
with no assumptions on the spectral bandwidth. In practical computa-
tions, the Zakharov equation is discretised in terms of a fixed number of
free wave modes which are subject to its own evolution equations. An
important computational method for solving the Zakharov equation was
formulated by West et al. (1987) and Dommermuth and Yue (1987) as the
Higher Order Spectral Method (HOSM). This relatively efficient numeri-
cal method is a pseudo-spectral type method where a series expansion in
the wave slope of the vertical velocity about the free surface is used. The
method does not contain any additional assumptions, except for the surface
being single-valued (i.e. no wave breaking). High efficiency and accuracy
makes HOSM a good tool for studying the nonlinear wave dynamics. Mei
et al. (2005) described in detail mathematical formulation and numerical
implementation of the HOSM for simulation of nonlinear wave–wave, wave–
current and wave–bottom interactions. The method has been subsequently
extended to include the presence of atmospheric forcing, finite depth varia-
tions, variable current and moving submerged bodies. In recent years, the
July 31, 2017 16:58 ws-book9x6 10666-main page 24

24 Introduction

HOSM has been applied by DNV GL (Norway) for studies on the rogue
waves (Bitner-Gregersen and Gramstad, 2015). For more detail review of
HOSM, the reader should consult Tanaka (2001).
Further extension of fully nonlinear wave modelling was reported in a se-
ries of papers by Bateman et al. (2001, 2003, 2012). The developed model
by [Link] on 11/12/25. Re-use and distribution is strictly not permitted, except for Open Access articles.

for realistic ocean spectra incorporates spread of wave energy across the fre-
quency and directional domains, well as the evolution of the largest waves.
This procedure provides full-field (both spatial and temporal) description
of the evolving water surface ζ(x, y, t), together with the associated wa-
ter particle kinematics u(x, y, z, t). Given a valid set of initial conditions,
Ocean Surface Waves Downloaded from [Link]

the applied calculation procedure is essentially exact and has been used to
provide sample calculations involving the JONSWAP spectra in relatively
deep water. However, model is not restricted to this flow regime and can
be successfully applied for any water depths. Although the modified non-
linear Schrödinger equations can be employed to describe the modulations
of a regular and irregular wave trains over some wave cycles, they are not
appropriate to the description of extreme waves arising in broad-band sea
states.
The engineering practise frequently requires the information concerning
very nonlinear ocean waves, including the limiting breaking waves. For uni-
directional waves propagating in a two-dimensional fluid, the time-marching
procedure allowing the evolution of the wave-field has been developed using
the inverse method in which potential function ϕ and stream function ψ
are adopted as the dependent variables rather than coordinates (x, z) of
the common physical plane method (Longuet-Higgins and Cokelet, 1976;
Cokelet, 1977). The solution is valid to the limiting breakers. However,
this type of solution is seldom applied in engineering design because it can
not easily be expanded to include the effects of directional spreading and
it requires significant computing resources.
The wave kinematics and loads due to extreme waves are very impor-
tant for offshore engineering operations, offshore structures and ships per-
formance. An accurate assessment of the maximum water particle velocities
beneath the high wave crest is required for drag force calculations, while
the maximum accelerations beneath the steepest section of the wave profile
are critical for inertial force calculations. The wave models used to deter-
mine the water particle kinematics associated with a measured time history
of the surface elevation ζ(t) are traditionally based either on a nonlinear
regular wave theory (for example a higher-order Stokes’ theory) or on an
unsteady linear wave theory (a Fourier transform approach).
July 31, 2017 16:58 ws-book9x6 10666-main page 25

1.3. Fundamentals of Description of Regular Waves 25

This is not the case for rough waves and near breaking waves when the
kinematics is not well predicted by any existing theory. A fully non-linear
unsteady wave model has been developed by Clamond and Grue (2001a, b)
which offers a substantial improvement in the modelling of steep wave kine-
matics. The model is based on potential theory and assumes the integration
by [Link] on 11/12/25. Re-use and distribution is strictly not permitted, except for Open Access articles.

of the kinematic and dynamic boundary conditions at the free surface of

very steep (rough) waves. The theoretical results compare favourably with
various available sets of experimental data. Measurements and computa-
tions show some deep water waves with fluid velocity up to 75% of the phase
speed (Grue et al., 2003). The PIV and LDA data demonstrate excellent
Ocean Surface Waves Downloaded from [Link]

agreement with the fully non-linear theory. The experimental data corre-
spond to the JONSWAP spectrum with ωp2 Hs /2g = 0.15 and the range
of the wave slope is 0.21 < εs < 0.34. A surprising result is that the
exponential profile ekz compares well with all measurements in deep water.
Accelerations at the sea surface are required to estimate the inertia
forces on offshore installations. Very few measurements exist of the real ac-
celeration below steep waves (see for example Bonmarin and Kjeldsen, 2000;
Grue and Jensen, 2006). In particular, Grue and Jensen reported a maxi-
mum negative vertical acceleration value at wave crest of about 1.1g, while
the upward vertical acceleration grows to about 1.5g in the front face of the
wave at the base below the overturning jet of the plunging breakers. The
fully non-linear theoretical model for accelerations in unsteady waves shows
very good agreement with experiments. It was shown that the convective
term is of the same order of magnitude and of opposite sign to the local
acceleration, and it cannot be neglected when estimating the acceleration
and forces (Jensen et al. 2007).

[Link] Boussinesq and Korteweg–de Vries equations for waves in

shallow water
In the coastal zone, the typical wave length, L, is much larger than water
depth h. Moreover, the wave height, H, is an appreciable fraction of h, so
the Ursell number becomes too large to allow the Stokes approach to be
used. This is also true for tsunami waves approaching the coast or for the
internal waves (Massel, 2015).
Let us define two small parameters δ and ν:
– relative wave length δ = kh
– relative wave amplitude ν = A/h.
July 31, 2017 16:58 ws-book9x6 10666-main page 26

26 Introduction

Therefore the Ursell number U becomes

 2
H L
U= (1.57)
h h
or
by [Link] on 11/12/25. Re-use and distribution is strictly not permitted, except for Open Access articles.

A
ν
U = 2(2π) h 2 ≈ 75 2 .
2
(1.58)
(kh) δ
The parameters δ and ν allow to distinguish between three different impor-
tant models of the long waves:
Ocean Surface Waves Downloaded from [Link]

a) nonlinear shallow water equation

ν  δ 2 , ν = O(1), U  75
b) Boussinesq and Korteweg–de Vries equation
O(ν) ≈ O(δ 2 ) < 1, U ∼ O(75)
c) linearised long wave equation
ν  δ 2 , U  75.

In the following we will further consider the case of the Boussinesq and
Korteweg–de Vries equations.7 The essence of the method related to these
equations lies in eliminating the z-dependence and adopting the character-
istic horizontal velocity either as the velocity at z = 0 or at the bottom
z = −h or the velocity averaged over the depth. When we use the hori-
zontal velocity ū averaged over the water depth, we obtain the Boussinesq
equations as follows (Whitham, 1974)

∂ζ ∂
+ [(h + ζ)ū] = 0 

∂t ∂x


. (1.59)
∂ ū ∂ ū ∂ζ 1 3 ∂ 3 ū 

+ ū +g = h


∂t ∂x ∂x 3 ∂x2 ∂t
For waves propagating in the given direction, Eq. (1.59) is reduced to the
Korteweg–de Vries (KdV) equation (Korteweg and de Vries, 1895). The
KdV equation in terms of the surface elevation takes the form
h2 p ∂3ζ
 
∂ζ p 3 ζ ∂ζ
+ gh 1 + + gh = 0, (1.60)
∂t 2 h ∂x 6 ∂x3
7 Joseph Valentin Boussinesq (1842–1929) was a French mathematician and physicist

who made significant contributions to the theory of hydrodynamics, vibration and heat.
He developed the mathematical theory to support Russell’s observations of the solitary
waves, and introduced the KdV equation in 1877 which was rediscovered by Diederik
Korteweg and Gustaw de Vries in 1895.
July 31, 2017 16:58 ws-book9x6 10666-main page 27

1.3. Fundamentals of Description of Regular Waves 27

while in terms of the velocity ū we have

∂ ū h2 p ∂ 3 ū
 
∂ ū p 3
+ gh + ū + gh = 0. (1.61)
∂t 2 ∂x 6 ∂x3
The KdV equation describes the evolution of weakly dispersive and weakly
by [Link] on 11/12/25. Re-use and distribution is strictly not permitted, except for Open Access articles.

nonlinear waves. It has the periodic and stable solutions, namely the soli-
tary waves and cnoidal waves. In particular, the surface elevation of the
solitary wave for small value of ν parameter is
r !
−2 3H
ζ(x, t) = H cosh (x − Cp t) , (1.62)
4h3
Ocean Surface Waves Downloaded from [Link]

in which
 1/2   
p H p 1 H
Cp = gh 1 + ≈ gh 1 + . (1.63)
h 2 h
The cnoidal wave form of the surface displacement resulting from the KdV
equation provides an opportunity to study dynamics of the simple linear
sinusoidal waves as well as solitary waves within the same framework of
the Jacobian elliptic functions (Abramowitz and Stegun, 1975). Thus the
cnoidal wave shape becomes (Massel, 1989)
    
1 E 2 2K
ζ(x, t) = H(x) 1− − m + cn (x − Cp t) , (1.64)
m K L
in which H(x) is the wave height, L is the wave length and cn(x) is the
Jacobian elliptic function, K and E are the complete elliptic integrals of
the first and second kind, respectively. They are functions of the elliptic
parameter m from the range 0 ≤ m ≤ 1 (Abramowitz and Stegun, 1975).
The phase speed Cp is given by
  
H 1 E
Cp2 = gh 1 + 2−3 −m . (1.65)
h m K
To define an elliptic parameter m, the following relationship is used
   2
H L 16
= mK2 . (1.66)
h h 3
It should be noted that if m = 0, the Eq. (1.64) becomes
H 2π
ζ(x, t) = cos (x − Cp t), (1.67)
2 L
which is a linear sinusoidal wave. For elliptic parameter m = 1, a limit of the
cnoidal wave is the solitary wave of infinite wave length — see Eq. (1.62).
For more systematic review of the solitary and cnoidal waves the reader
August 14, 2017 10:22 ws-book9x6 10666-main page 28

28 Introduction
by [Link] on 11/12/25. Re-use and distribution is strictly not permitted, except for Open Access articles.
Ocean Surface Waves Downloaded from [Link]

Fig. 1.6 Comparison of solitary and cnoidal wave displacements.

should refer to the paper by Fenton (1979). For illustration, the comparison
of the solitary and cnoidal wave profiles is shown in Fig. 1.6, in which the
cnoidal profile corresponds to the elliptic parameter m = 0.99.
There are several proposals for the range of application of various wave
theories. The applicability of the second order Stokes’ perturbation method
and cnoidal wave theory for determination of the limiting wave height in
shallow water is considered in some detail in Section 7.2.8. For the practical
purposes, the ranges of validity of wave theories has been given by Druet
(1978) and Chakrabarti (1987) in terms of wave shallowness (h/T 2 ) and
wave steepness (H/T 2 ).

1.4 Methods of Description of Ocean Random Waves

1.4.1 Preliminaries
In looking at ocean surface waves one notices both their randomness and
their quasi-regularity. The wave profile changes constantly with time and
in a random fashion. Consequently, the properties of waves are not readily
August 14, 2017 10:22 ws-book9x6 10666-main page 29

1.4. Methods of Description of Ocean Random Waves 29

by [Link] on 11/12/25. Re-use and distribution is strictly not permitted, except for Open Access articles.
Ocean Surface Waves Downloaded from [Link]

Fig. 1.7 Methods of processing of random time series ζ(t).

defined on a wave-by-wave basis. As the fundamental property of surface

waves, induced by wind, is their irregularity, the prediction of wave parame-
ters can be achieved only through stochastic analysis of the sea data, which
span three basic domains, namely frequency domain, probability domain,
and time domain, as was demonstrated in Fig. 1.7 for random time series
ζ(t).
Frequency analysis deals mainly with an evaluation of the distribution
of wave energy (potential or kinetic) among various frequencies and direc-
tions. There are two main methods for the development of the frequency
spectrum. The traditional method is based on the Fourier transform of
the auto-correlation function which is a measure of the connection between
the two values ζ(t) and ζ(t + τ ) of the random variable ζ. The theoretical
background of this transformation is provided by the Wiener–Khinchine
Theorem. The transformation of the auto-correlation function gives the
spectral density function of the particular variable, and the transformation
of cross-correlation function results in cross-spectral density function, which
is a complex function. The real part of the cross-spectral function is called
the co-spectrum and the imaginary part is called the quadrature spectrum.
A more general spectral representation of surface waves is obtained when
the frequency, as well the directional energy distribution, is taken into ac-
count. The resulting spectrum is called the frequency-directional spectrum.
July 31, 2017 16:58 ws-book9x6 10666-main page 30

30 Introduction

The second spectral method is based on the straightforward transferring of

the time series into its Fourier components. This technique, known as Fast
Fourier Transform (FFT), was first introduced in a new form by Cooley
and Tukey (1965). It reduces the number of calculations from a number
proportional to n2 (n is a number of samples) to a number approximately
by [Link] on 11/12/25. Re-use and distribution is strictly not permitted, except for Open Access articles.

proportional to n log n and has revolutionised the spectral analysis of time

series. One major weakness of the Fourier-based spectral methods is the
assumption of linear superposition of wave components, which, in fact,
are phase-coupled via the nonlinear dynamics present in ocean waves. It
should be noted that the Fourier transform provides an information on the
Ocean Surface Waves Downloaded from [Link]

frequency content, however, the information on the frequency localisation in

time is essentially lost in the process. In addition, Fourier spectral analysis
should, strictly speaking, be used for periodic and stationary processes
only. When the time localisation of the spectral components is needed, the
transform of time series which provides the time-frequency representation
of the signal is the Wavelet transform (see Section [Link]).
In the probability domain, the particular wave parameters, e.g. surface
displacements at a given time, wave amplitudes, wave heights, wave periods,
etc. are considered as elementary random events. The probability approach
is easy to understand when we are dealing with digitised data. The digitised
data of a particular parameter form a set of random realisations of a random
variable, when the time sequence of the parameter is exterminated. The
final results in this approach are expressed in terms of probability density
functions, distribution functions and moments (see Chapter 4).
The simplest statistics are obtain when one assumes that the observed
wave field is a result of a linear superposition of a large number of dynami-
cally independent waves. This is the basis of the Gaussian model, in which
two first moments are sufficient for the complete statistical description of
the wave field. However, in the real ocean, due to nonlinear interactions
between spectral harmonics and energy dissipation processes, a substantial
departure from the Gaussian model is observed. Ocean waves should be
treated in many cases as example of the non-Gaussian stochastic processes.
For processing of the nonstationary time series within time domain,
the Hilbert transform becomes very useful tool. The Hilbert transform of
a real-valued time domain signal ζ(t) is another real-valued time domain
signal ξ(t) such that
η(t) = ζ(t) + iξ(t) = A(t) exp(iθ(t)), (1.68)
where ξ(t) = H[ζ(t)] is the Hilbert transform of function ζ(t), and function
η(t) forms an analytical signal, also in time domain. This is a substantial
July 31, 2017 16:58 ws-book9x6 10666-main page 31

1.4. Methods of Description of Ocean Random Waves 31

Table 1.2 A comparative summary of Fourier, Wavelet and Hilbert trans-

forms methodology (adapted from Huang and Shen, 2005).
Fourier Wavelet Hilbert

Basis a priori a priori adaptive

Frequency convolution convolution differentiation
by [Link] on 11/12/25. Re-use and distribution is strictly not permitted, except for Open Access articles.

global regional local,

uncertainty uncertainty certainty
Presentation energy- energy-time- energy-time-
frequency frequency frequency
Nonlinear no no yes
Nonstationary no yes yes
Ocean Surface Waves Downloaded from [Link]

Feature no discrete no; yes

Extraction continuous yes
Theoretical base theory complete theory complete empirical

difference when comparing with the Fourier transform which transforms

signal ζ(t) into a complex-valued frequency domain signal. Representation
(1.68) provides the best local fit of an amplitude and phase for function
ζ(t) varying in time, as well as give some information on the instantaneous
frequency nature. A comparative summary of Fourier transform, Wavelet
transform and Hilbert transform methods are given in Table 1.2.

1.4.2 Basic definitions and analysis of random time series

Let us begin with an ensemble of k wave records {ζk (t)} taken under identi-
cal macroscopic conditions, i.e. position on the ocean surface, water depth,
mean wind velocity, air-sea temperature, etc. Even under identical condi-
tions, we cannot expect that these wave records will be identical or even
closely similar in detail. The family {ζk (t)} represents k realisations of the
stochastic process ζ(t). For a given k, ζ(t) is a function of time t, while
when t = t1 , ζk (t1 ) is a random variable.
Stochastic processes may belong to one of three categories: a) stationary
and ergodic, b) stationary, and c) non-stationary. Random process (or
random function) is stationary in the wide sense if

E [ζ(t)] = ζ = const, (1.69)

and

K(t1 , t2 ) = K(t1 − t2 ) = E [ζ (t1 ) ζ (t2 )] = K(τ ), τ = t1 − t2 , (1.70)

July 31, 2017 16:58 ws-book9x6 10666-main page 32

32 Introduction

in which E[ ] denotes the mean of ζ and K( ) is an auto-correlation func-

tion. Strictly speaking, a random process is stationary if all statistical
moments are invariant. These two definitions of stationarity coincide when
ζ is a Gaussian, in which case all the statistics of ζ are completely de-
termined by the first and second moments. This strict definition is often
by [Link] on 11/12/25. Re-use and distribution is strictly not permitted, except for Open Access articles.

relaxed and stationarity in wide sense is assumed.

In general, using the ensemble of wave records {ζk (t)}, we can develop
any function of ζ, say F , to give F {ζk (t)}. To be more specific, we select
time t = t1 in the family {ζk (t)}. When F is the ζ value itself, then
averaging F {ζk (t1 )} over k results in the ensemble mean of the process at
Ocean Surface Waves Downloaded from [Link]

t = t1 , i.e.
PN
k=1 ζk (t1 )
E[F {ζk (t1 )}]k = E[ζk (t1 )]k = lim . (1.71)
N →∞ N
The condition N → ∞ is only conceptual since N is usually finite.
2
When F {ζk (t1 )} ≡ [ζk (t1 )] , then averaging F {ζk (t1 )} over k leads to
the variance at t = t1
PN 2
k=1 [ζk (t1 )]
n o
2
E[F {ζk (t1 )}]k = E [ζk (t1 )] = lim . (1.72)
k N →∞ N
Let us now define F as

1
 if a < ζk (t1 ) ≤ b
F {ζk (t1 )} = (1.73)

0 otherwise.

Averaging over the ensemble E[F {ζk (t1 )}]k can be interpreted as the en-
semble probability that ζk (t1 ) falls on the interval from a to b at t = t1 .
Repetition of the above averaging for the different times helps us to
obtain the different numerical values for the statistics. However, the re-
peated observations technique which provides us with an ensemble of k wave
records can be attained in laboratory wave tanks, but is inapplicable to ob-
servations of wave phenomena in field experiments. To overcome these
difficulties, the ergodic theorem is usually invoked. This allows the ensem-
ble averages to be replaced with time averages. The ergodic theorem states
that (Ochi, 1998): if ζ(t) is an ergodic stationary random function, then the
statistics obtained by ensemble averages at a given time t = t∗ are identical
to the corresponding statistics computed by the time averaging for any given
realisation k = k∗ . Hence, the ergodic stationary process should satisfy the
July 31, 2017 16:58 ws-book9x6 10666-main page 33

1.4. Methods of Description of Ocean Random Waves 33

following equality for record of time T

PN
k=1 F {ζk (t = t∗ )}
E[F {ζk (t = t∗ )}]k = lim = E[F {ζk=k∗ (t)}]t
N →∞ N
Z T
1
F {ζk=k∗ (t)} dt.
by [Link] on 11/12/25. Re-use and distribution is strictly not permitted, except for Open Access articles.

= lim (1.74)
T →∞ 2T −T

If above condition is limited to the mean and covariance, then the stochastic
process is called a weakly ergodic random process. Therefore we have
N Z T
1 X
xk (t) = lim x(t)dt = x (1.75)
Ocean Surface Waves Downloaded from [Link]

N T →∞ 0
k=1

and
N Z T
1 X
[xk (t) − x][xk (t + τ ) − x] = lim [x(t) − x][x(t + τ ) − x]dt. (1.76)
N T →∞ 0
k=1

We can say that, while stationary processes are a limited subset of random
processes, ergodic processes in their turn are an even more limited subset
of stationary processes.
The significance of the ergodic theorem is that it enables us to develop
the statistics of the process ζ(t) using one, sufficiently long realisation.
However, it is never possible to demonstrate ergodicity for ocean waves
since experiments cannot be exactly repeated in the ocean as they can in
the laboratory.
It can be shown that for a stationary wave process ζ(t) to be ergodic
it is sufficient that its auto-correlation function K(τ ) satisfies the following
condition
K(τ ) = 0 at τ → ∞. (1.77)
Let us now demonstrate the applicability of the ergodic theorem and con-
dition (1.77) for a simple process. We assume that we have an ensemble of
records of a process {ζk (t)} = zk . This means, that for a particular k, the
process ζk (t) is constant and equal zk . It is clear that process is station-
ary. For any chosen time t, any statistic, for example the mean, calculated
across the ensemble gives the same number. However, when any single
record ζk=k∗ (t) is chosen at random and its time average calculated (Kins-
man, 1965)
Z T
1
E[ζk=k∗ (t)]t = lim ζk=k∗ (t)dt, (1.78)
T →∞ 2T −T
July 31, 2017 16:58 ws-book9x6 10666-main page 34

34 Introduction

then, it is obvious that

E[ζk (t = t∗ )]k 6= E[ζk=k∗ (t)]t . (1.79)
Hence, this extremely simple process is stationary, but not ergodic. The
condition (1.77) is clearly not satisfied as
by [Link] on 11/12/25. Re-use and distribution is strictly not permitted, except for Open Access articles.

Z T
1
K(τ ) = E [ζk=k∗ (t) · ζk=k∗ (t + τ )] = lim zk2∗ dt = zk2∗ . (1.80)
T →∞ 2T −T

In the following we assume that the ergodic property holds for the stochastic
process discussed in this text. Therefore, a single record ζ(t) will be used
instead of an ensemble of records {ζk (t)}.
Ocean Surface Waves Downloaded from [Link]

The representation (1.28) is quite general but still a deterministic one.

We will use it now for evaluation of the random wave model. We assume
that ζ(t) is a stationary ergodic process satisfying the conditions (1.69)
and (1.70). Without loosing generality we can say that ζ = 0. Under
the above assumptions there exists stochastic process Ã(u, v, w), such that
the process ζ(x, y, t) has the representation in form of the Fourier–Stjeltjes
integral8 (Yaglom, 1962)
Z ∞Z ∞Z ∞
ζ(x, y, t) = exp [i(ux + vy + wt)] dÃ(u, v, w). (1.81)
−∞ −∞ −∞

A three-dimensional stochastic process Ã(u, v, w) possesses uncorrelated

increments which means that differences [ Ã(u00 , v 00 , w00 ) − Ã(u0 , v 0 , w0 ) ] for
(u00 , v 00 , w00 ) and (u0 , v 0 w0 ) within a given cube in (u, v, w) space are uncor-
related with differences within any other non-overlapping cube.
The integrals (1.28) and (1.81) look similar, but they are significantly
different. The former representation is a deterministic integral, while the
latter is a representation of random function and should be defined by
probabilistic limits.
The general representation given by Eq. (1.81) depends only on the
assumption that process ζ(x, y, t) is ergodic and conditions (1.69) and (1.70)
are valid. It also is important to note that this representation exists without
assuming that the sea surface is composed of a superposition of linear waves.
When the surface displacements ζ(x, t) is assumed to be a stationary
random function of time and position, Eq. (1.81) can be represented as
Z Z
ζ(x, t) = exp[i(kx − ωt)]dÃ(k, ω), (1.82)
k ω
8 Fourier–Stjeltjes
integrals are a generalisation of both Fourier series and Fourier in-
tegrals. Applying the Fourier–Stjeltjes integral we can represent functions which are
not periodic and which are too large at infinity to be represented by ordinary Fourier
integral.
July 31, 2017 16:58 ws-book9x6 10666-main page 35

1.4. Methods of Description of Ocean Random Waves 35

in which
0 for k, ω 6= k0 , ω 0

∗ 0 0
E[dÃ(k, ω)dà (k , ω )] = (1.83)
Ψ̂(k, ω)dkdω for k = k0 , ω = ω 0 .
It should be noted that for stationary and homogeneous wave field,
wavenumber-frequency spectrum Ψ̂(k, ω) is real and positive for all values
by [Link] on 11/12/25. Re-use and distribution is strictly not permitted, except for Open Access articles.

of k and ω.
If the wave components obey the dispersion relation appropriate to in-
finitesimal gravity waves as

1/2 h
1/2 i
w 2 = g u2 + v 2 tanh u2 + v 2 h , (1.84)
Ocean Surface Waves Downloaded from [Link]

we should impose a restriction on wavenumber k and wave frequency ω. In

particular, the parameters (u, v, w) can be replaced by the two parameters
(ω, θ), or alternatively by k = (u, v), and dÃ(u, v, w) is zero whenever the
above dispersion relation is not satisfied. Therefore Eq. (1.81) becomes the
Fourier–Stjeltjes integral representation in space (ω, θ) or (kx , ky )
Z ∞Z π
ζ(x, y, t) = exp [i (k · x − ωt)] dà (ω, θ) (1.85)
−∞ −π
or
Z ∞ Z ∞
ζ(x, y, t) = exp [i (k · x − ωt)] dB (k) . (1.86)
−∞ −∞
The models implied by Eqs. (1.85) and (1.86) permit both random phases
and random amplitudes.
By definition, the cross-correlation function K(X, Y, τ ) is
K(X, Y, τ ) = E [ζ(x, y, t) · ζ ∗ (x + X, y + Y, t + τ )]
Z ∞Z ∞Z π Z π
= E {exp [i(ux + vy) − iωt]
−∞ −∞ −π −π

× exp [−i [u0 (x + X) + v 0 (y + Y )] + iω 0 (t + τ )]

× dA(ω, θ) dA∗ (ω 0 , θ0 )} , (1.87)

in which (∗) denotes the conjugate of the complex variable.
The dA(ω, θ) = dÃ(ω, θ) and dA∗ (ω 0 , θ0 ) = dÃ∗ (ω 0 , θ0 ) and these are
uncorrelated unless u = u0 , v = v 0 and ω = ω 0 . Using the fact that variance
of ζ(x, y, t) is finite, it can be shown that (Yaglom, 1962)
E [dA(ω, θ) dA∗ (ω 0 , θ0 )] = Ŝ(ω, θ) δ(ω − ω 0 )δ(θ − θ0 )dωdω 0 dθdθ0 , (1.88)
in which Ŝ(ω, θ) is a two-dimensional spectral density function of surface
waves.
July 31, 2017 16:58 ws-book9x6 10666-main page 36

36 Introduction

From Eq. (1.88) and the known relation for Dirac delta (Lighthill, 1975)
Z ∞
f (y) δ(x − y)dy = f (x), (1.89)
−∞

Eq. (1.87) yields

by [Link] on 11/12/25. Re-use and distribution is strictly not permitted, except for Open Access articles.

Z ∞ Z π
K(X, Y, τ ) = exp [−i (uX + vY ) + iωτ ] Ŝ(ω, θ)dωdθ
−∞ −π
Z ∞ Z π
= exp [−i k · X + iωτ ] Ŝ(ω, θ)dωdθ, (1.90)
−∞ −π

where X = (X, Y ).
Ocean Surface Waves Downloaded from [Link]

In a similar way, from Eq. (1.87) we have

Z ∞Z ∞
K(X, Y, τ ) = exp [−i k · X + iωτ ] Ψ̂ (k) dkx dky , (1.91)
−∞ −∞

in which Ψ̂ (k) is a wavenumber spectrum.

Let us now consider wave observation at a single point P (x, y) where
X = Y = 0. Expression (1.90) becomes the auto-correlation function
Z ∞Z π
K(τ ) = Ŝ(ω, θ)eiωτ dωdθ (1.92)
−∞ −π
or
Z ∞
K(τ ) = S̃(ω)eiωτ dω, (1.93)
−∞

in which
Z π
S̃(ω) = Ŝ(ω, θ)dθ. (1.94)
−π

The S̃(ω) is a frequency spectral density function or frequency spectrum

of the ocean surface. It represents the distribution of wave energy in the
frequency domain. Hence, the variance of surface waves is
Z ∞Z π Z ∞
E ζ 2 = σζ2 = K(0) =
 
Ŝ(ω, θ)dωdθ = S̃(ω)dω. (1.95)
−∞ −π −∞

By taking the inverse Fourier transform, we have

Z ∞
1
S̃(ω) = K(τ )e−iωτ dτ. (1.96)
2π −∞
The function S̃(ω) is called the symmetric function, as S̃(ω) = S̃(−ω). This
spectrum is usually used in theoretical considerations.
July 31, 2017 16:58 ws-book9x6 10666-main page 37

1.4. Methods of Description of Ocean Random Waves 37

In experimental practise, the so called nonsymmetric spectrum S(ω) is used.

It can be seen from Eqs. (1.93) and (1.96) that for ergodic stationary real
random process ζ(t) with zero mean, the spectral density function S̃(ω)
and the auto-correlation function K(τ ) are a Fourier transform pair. This
is called the Wiener–Khinchine Theorem.9 A few specific examples of auto-
by [Link] on 11/12/25. Re-use and distribution is strictly not permitted, except for Open Access articles.

correlations K(τ ) and spectral densities S̃(ω) are given in Tables 1.3 and
1.4.
For most processes, K(τ ) → ∞. The decay time which characterise this
behaviour is called correlation scale. For example, from Table 1.3 follows
that a Markov process (see Chapter 4) has the auto-correlation function as
Ocean Surface Waves Downloaded from [Link]

 
|τ |
K(τ ) = K0 exp − . (1.97)
τ0
For τ = τ0 we have K(τ0 )/τ0 = 1/e ≈ 0.368. Taking the limit τ0 → 0 in
Eq. (1.97) gives the auto-correlation of a white noise random process (see
Table 1.3).
The auto-correlation function K(τ ) for surface elevation ζ(t) is a real
and even function, i.e. K(τ ) = K(−τ ). Therefore, the spectral density
S̃(ω) is also a real and even function, and Eqs. (1.93) and (1.96) can be
written as follows Z
∞ Z ∞
K(τ ) = S̃(ω)eiωτ dτ = 2 S̃(ω) cos(ωτ )dω (1.98)
−∞ 0
and Z ∞
1 ∞
Z
1
S̃(ω) = K(τ )e−iωτ dτ = K(τ ) cos(ωτ )dτ. (1.99)
2π −∞ π 0
We define nonsymmetric spectrum S(ω) as

 2S̃(ω) for ω ≥ 0
S(ω) = (1.100)
0 for ω < 0


and
1
S̃(ω) = [S(ω) + S(−ω)]. (1.101)
2
Therefore, Eqs. (1.98) and (1.99) take the forms frequently used in appli-
cations Z ∞
2 ∞
Z
K(τ ) = S(ω) cos(ωτ )dω S(ω) = K(τ ) cos(ωτ )dτ. (1.102)
0 π 0
9 Norbert Wiener proved this theorem for the case of a deterministic function in 1930

and Aleksandr Khinchine formulated an analogous result for stationary stochastic process
in 1934. The Wiener–Khinchine Theorem states that for a wide-sense process with an
auto-correlation function K(τ ) finite at every lag τ , there exists a monotonic function
S̃(ω) in the frequency domain.
July 31, 2017 16:58 ws-book9x6 10666-main page 38

38 Introduction

Table 1.3 Auto-correlation functions.

Z ∞
Process K(τ ) = S̃(ω)eiωτ dω
by [Link] on 11/12/25. Re-use and distribution is strictly not permitted, except for Open Access articles.

−∞

Analytical expression Graphical representation

Ocean Surface Waves Downloaded from [Link]

White noise K0 δ(τ )

 
τ
sin τ0
Finite-band K0 τ
τ0
white noise

K0 exp −ατ 2


Gaussian

 
|τ |
Markov K0 exp − τ
0

Cosinusoidal K0 cos (ω0 τ )

signal with
random phase
July 31, 2017 16:58 ws-book9x6 10666-main page 39

1.4. Methods of Description of Ocean Random Waves 39

Table 1.4 Spectral density functions.

1
Z ∞
S̃(ω) = K(τ )e−iωτ dτ
2π −∞
by [Link] on 11/12/25. Re-use and distribution is strictly not permitted, except for Open Access articles.

Analytical expression Graphical representation

Ocean Surface Waves Downloaded from [Link]

K0
π

K0 τ0 for |ω| ≤ τ0−1

0 for |ω| > τ0−1

ω2
 
K0
√ exp −
πα 4α

2K0 τ0
π 1 + (ωτ0 )2

K0
2
[δ(ω − ω0 ) + δ(ω + ω0 )]
July 31, 2017 16:58 ws-book9x6 10666-main page 40

40 Introduction

Let us consider for a simple harmonic wave ζ(x, t) = a cos(kx − ωt) with
ω = ω0 (see Tables 1.3 and 1.4). Thus we obtain
K0
K(τ ) = K0 cos(ω0 τ ), S̃(ω) = [δ(ω − ω0 ) + δ(ω + ω0 )]. (1.103)
2
For ω0 > 0 we have
by [Link] on 11/12/25. Re-use and distribution is strictly not permitted, except for Open Access articles.

a2
S(ω) = K0 δ(ω − ω0 ) with K0 = = ζ 2. (1.104)
2
Let us consider now the wavenumber spectrum Ψ̂ (k) for X = Y = 0. From
Eq. (1.91) we have
Ocean Surface Waves Downloaded from [Link]

Z ∞Z ∞ Z ∞Z π
K(τ ) = eiωτ Ψ̂ (k) dkx dky = eiωτ Ψ̂ (k, θ) kdkdθ,
−∞ −∞ 0 −π
(1.105)
in which k = |k|.
In deep water, the dispersion relation (1.17) becomes ω 2 = gk and
ω3
kdk = 2 2 dω. After substituting this into Eq. (1.105) we obtain
g
Z ∞Z π  3 
2ω
K(τ ) = eiωτ Ψ̂ (k, θ) |k=ω 2 /g dωdθ. (1.106)
−∞ −π g2
Now Eq. (1.106) can be compared with Eq. (1.92), i.e.
2ω 3
Ŝ(ω, θ) = Ψ̂ (k, θ) (1.107)
g2
and
π
2ω 3
Z
S(ω) = Ψ̂ (k, θ) dθ. (1.108)
g2 −π

Thus, the frequency spectrum S(ω) represents an integral over the

wavenumber spectrum Ψ̂(k, θ) for all directions θ, when k is constant and
equals ω 2 /g.
The above definitions can be immediately extended to a set of two ran-
dom processes, simultaneously recorded at the same point or at two different
points. As an example, we will first consider the surface elevation ζ(t) and
the horizontal component of orbital velocity u(t) recorded at point P (0, 0).
Using the linear wave theory we express velocity u(t) as corresponding to
surface elevation ζ(t), given by Eq. (1.85), i.e.

Z ∞ Z π
gk cos θ cosh k(z + h)
u(0, 0, z, t) = exp (−iωt) dÃ(ω, θ), (1.109)
−∞ −π ω cosh kh
July 31, 2017 16:58 ws-book9x6 10666-main page 41

1.4. Methods of Description of Ocean Random Waves 41

in which z is a level at which velocity u is recorded.

The cross-correlation function Kζu (τ ) takes the form

∞ ∞ ∞ ∞
gk 0 cos θ0
Z Z Z Z
Kζu (τ, z) = E [ζ(t) u(t + τ )] =
by [Link] on 11/12/25. Re-use and distribution is strictly not permitted, except for Open Access articles.

−∞ −∞ −π −π ω0
cosh k(z + h)
× E {exp [i(ω 0 − ω) + iω 0 τ ]}
cosh kh
×E {dA(ω, θ)dA∗ (ω 0 , θ0 )} (1.110)
Ocean Surface Waves Downloaded from [Link]

or
Z ∞ Z π 
gk cos θ cosh k(z + h)
Kζu (τ, z) = Ŝ(ω, θ)dθ eiωτ dω. (1.111)
−∞ −π ω cosh kh

Defining Sζu (ω, z) as

Z π
gk cos θ cosh k(z + h)
Sζu (ω, z) = Ŝ(ω, θ)dθ, (1.112)
−π ω cosh kh
the Wiener–Khinchine Theorem yields
Z ∞
Kζu (τ, z) = Sζu (ω, z)eiωτ dω (1.113)
−∞

and
Z ∞
1
Sζu (ω, z) = Kζu (τ, z)e−iωτ dτ. (1.114)
2π −∞

In general Sζu (ω, z) is a complex function, i.e.

Sζu (ω, z) = Cζu (ω, z) + i Qζu (ω, z), (1.115)
in which the real part Cζu (ω, z) is referred to as the co-spectrum, while the
imaginary part Qζu (ω, z) is referred to as the quadrature spectrum. The
amplitude spectrum of Sζu (ω, z) becomes
q
2 2
|Sζu (ω, z)| = [Cζu (ω, z)] + [Qζu (ω, z)] , (1.116)
and the phase spectrum is given by
 1/2
−1 Qζu (ω, z)
ϕ(ω, z) = tan . (1.117)
Cζu (ω, z)
The various properties of the co-and quadrature spectra can be found else-
where (for example, Bendat and Piersol, 1986; Ochi, 1998).
July 31, 2017 16:58 ws-book9x6 10666-main page 42

42 Introduction

It should be noted that the two-dimensional spectrum of surface waves

Ŝ(ω, θ) is commonly represented as a product of a frequency spectrum S(ω)
and directional spreading function D(θ, ω), i.e.
Ŝ(ω, θ) = S(ω) D(θ, ω). (1.118)
Function D(θ, ω) will be discussed in detail in Section 3.4.
by [Link] on 11/12/25. Re-use and distribution is strictly not permitted, except for Open Access articles.

After substituting Eq. (1.118) into Eq. (1.112) and using Eq. (1.115),
we obtain Z π 
gk cosh k(z + h)
Cζu (ω, z) = S(ω) D(θ; ω) cos θ dθ 

ω cosh kh

−π
. (1.119)

Ocean Surface Waves Downloaded from [Link]


Qζu (ω, z) = 0


The vanishing of the quadrature spectrum Qζu (ω, z) is a consequence of the

assumption that sea surface oscillation ζ(t) and wave-induced horizontal
velocity u(t) are given in the same point and that there is no time delay
between them. This conclusion is valid for the symmetrical directional
spreading function, i.e. when D(θ) = D(−θ).
We now consider the case when surface elevation is simultaneously
recorded in two points, i.e. point P (x, y, ) and P1 (x + X, y + Y ). Eq. (1.90)
for cross-correlation function yields
Z ∞Z π
K(X, Y, τ ) = exp {−ikd cos (θ − θ0 ) + iωτ } Ŝ(ω, θ)dωdθ,
−∞ −π

√ (1.120)
in which d = X 2 + Y 2 and arctan θ0 = Y /X. Now the Wiener–Khinchine
Theorem gives
Z π
S(ω; d, θ0 ) = exp [−ikd cos (θ − θ0 )] Ŝ(ω, θ)dθ. (1.121)
−π
Using the identity (Abramowitz and Stegun, 1975)
∞
X
exp [ikd cos (θ − θ0 )] = m im Jm (kd) cos m(θ − θ0 ), (1.122)
m=0
in which 0 = 1 and m = 2 (m ≥ 1), Jm (x) is a first kind of Bessel
function of m order and separating the real and imaginary parts we obtain
co-spectrum C(ω; d, θ0 ) and quadrature spectrum Q(ω; d, θ0 ) in the form
∞
X
C(ω; d, θ0 ) = J0 (kd) + 2 (−1)m J2m (kd) cos [2m(θ − θ0 )] (1.123)
m=1
and
∞
X
−Q(ω; d, θ0 ) = 2 (−1)m J2m+1 (kd) cos [(2m + 1)(θ − θ0 )] . (1.124)
m=0
July 31, 2017 16:58 ws-book9x6 10666-main page 43

1.4. Methods of Description of Ocean Random Waves 43

1.4.3 Wave energy balance in spectral form

Wave generation models provide the linkage between the wavenumber-
frequency spectrum of air pressure fluctuations Ψ̂a (k, ω) and the wavenum-
ber spectrum of surface oscillations Ψ̂ (k, x, t) in a simpler case when waves
propagate in deep water or in water of finite but constant depth. Generally,
by [Link] on 11/12/25. Re-use and distribution is strictly not permitted, except for Open Access articles.

when water depth changes in space, and island or artificial obstacles are
present, the wind generated waves will undergo several other processes such
as refraction, diffraction and dissipation due to bottom friction, whitecap-
ping and breaking in the surf zone.
Taking into account the random character of wave motion, it seems that
Ocean Surface Waves Downloaded from [Link]

the most appropriate method of evaluation of waves, in time and space

are the spectral methods. Resulting functions of these methods, i.e. the
spectral densities Ŝ(ω, θ) or Ψ̂(k), provide the distribution of wave energy
in the frequency or wavenumber space.
In order to develop an energy balance equation we assume for a mo-
ment that the dispersion relation does not depend on time, but is rather
a slowly changing function of space coordinates. It can also depend on local
properties such as water depth, current velocity or ambient density current,
i.e.
ω = Ω [k, f (x, h, u, %, . . .)] . (1.125)
If the medium itself is moving with velocity u, the frequency of waves
passing a field point is
1/2
ω = Ω (k, x) = σ + k · u = [gk tanh(kh)] + k · u. (1.126)
Usually the quantity ω is called the observed or apparent frequency, while
σ is the intrinsic frequency whose functional dependence on k is known as
a classical dispersion relation
σ 2 = gk tanh(kh). (1.127)
It is well known that the number of wave crests passing a given point per
unit time is determined by the kinematics conservation equations (Massel,
1989)
∂k ∂ki ∂kj
+ ∇h ω = 0, − = 0, i, j = 1, 2. (1.128)
∂t ∂xj ∂xi
Substituting Eq. (1.125) into Eq. (1.128) gives
∂ki ∂kj ∂Ω
+ Cgj (x, k, t) =− , (1.129)
∂t ∂xi ∂xi
July 31, 2017 16:58 ws-book9x6 10666-main page 44

44 Introduction

∂Ω
in which Cgj = is a group velocity. Eq. (1.129) can be rewritten as
∂kj
dki ∂Ω dxi ∂Ω
=− and = . (1.130)
dt ∂xi dt ∂ki
Equations (1.130) are identical to the Hamilton equations in mechanics.
by [Link] on 11/12/25. Re-use and distribution is strictly not permitted, except for Open Access articles.

The vectors x should be interpreted as coordinates and vectors k represent

the moments, and the frequency Ω (k, x, t) is taken to be the Hamiltonian
(Whitham, 1974).
For free surface waves, the mean densities of potential and kinetic energy
are equal. Therefore, 2Ψ̂(k, x, t) is a mean spectral density of total wave
Ocean Surface Waves Downloaded from [Link]

energy in the unit volume dk · dx in phase space (k, x) at a given time t.

According to Liouville’s theorem, this energy should be conservative during
volume evolution
dΨ̂ ∂ Ψ̂ ∂xi ∂ Ψ̂ ∂ki ∂ Ψ̂
= + + =0 (1.131)
dt ∂t ∂t ∂xi ∂xi ∂ki
or
∂ Ψ̂ ∂Ω ∂ Ψ̂ ∂Ω ∂ Ψ̂
+ − = 0, (1.132)
dt ∂ki ∂xi ∂xi ∂ki
when Eq. (1.130) is used.
If the wave field is subjected to processes of generation, dissipation,
nonlinear interaction between spectral components and other possible in-
teractions with atmospheric boundary layer and various ocean movements
(currents, internal waves, etc.), Eq. (1.132) should be supplemented by
a source-sink term at the right-hand side, i.e.
∂ Ψ̂ ∂Ω ∂ Ψ̂ ∂Ω ∂ Ψ̂
+ − = Q (k, x, t) , (1.133)
dt ∂ki ∂xi ∂xi ∂ki
P
in which Q = i Qi .
The first term at the left-hand side of Eq. (1.133) expresses the local
evolution of the spectrum in time, while the second term represents the
evolution of the spectrum for the horizontally non-homogeneous wave field.
This term shows that energy is transported at the group velocity. The third
term reflects the effects of refraction and shoaling due to a non-horizontal
bottom or due to current. The spectral energy balance equation is also
known as the radiative transfer equation or the transport equation. In
particular, an analogy between the sea and a kinetic gas, which explains
the name kinetic equation, is used by Russian oceanographers.
In engineering and oceanographic practise the two-dimensional
frequency-directional spectrum Ŝ (ω, θ; x, t) is frequently used, rather than
July 31, 2017 16:58 ws-book9x6 10666-main page 45

1.4. Methods of Description of Ocean Random Waves 45

the wavenumber spectrum Ψ̂ (k; x, t) (see Chapter 3). To evaluate spectrum

Ŝ (ω, θ; x, t), we adopt the following variables transformation
k1 = χ (ω, x) cos θ, k2 = χ (ω, x) sin θ, (1.134)
where function χ (ω, x) = k is a solution of the dispersion relation.
by [Link] on 11/12/25. Re-use and distribution is strictly not permitted, except for Open Access articles.

Therefore variables transformation (1.134) now becomes

Ψ̂ (k; x, t) = J Ŝ (ω, θ; x, t) = f (ω, θ; x, t) , (1.135)
in which J is a Jacobian of the transformation, i.e.
 −1
∂(ω, θ) ∂(k1 , k2 )
Ocean Surface Waves Downloaded from [Link]

J= = . (1.136)
∂(k1 , k2 ) ∂(ω, θ)
Thus
−1
∂χ



 ∂ω cos θ, −χ sin θ 

 Cg
J= = . (1.137)
 ∂χ  χ
sin θ, −χ cos θ

 

∂ω
After substituting Eq. (1.135) into Eq. (1.133) we obtain (Krasitskii, 1974)
 
∂f ∂f ∂f ∂ ln χ ∂ ln χ ∂f
+ Cg cos θ + Cg sin θ + Cg cos θ − sin θ
∂t ∂x ∂y ∂y ∂x ∂θ

Cg
= Q (ω, θ, x, t) . (1.138)
χ
As Cp = ω/χ = ω/k, we can rewrite Eq. (1.138) as follows
∂ ∂ ∂
(Cp Cg Ŝ) + Cg cos θ (Cp Cg Ŝ) + Cg sin θ (Cp Cg Ŝ)
∂t ∂x ∂y
 
Cg ∂Cp ∂Cp ∂
+ sin θ − cos θ (Cp Cg Ŝ)
Cp ∂x ∂y ∂θ

= Cp Cg Q (ω, θ, x, t) , (1.139)
P
in which Q = i Qi . Some solutions of Eq. (1.138) (or (1.139)) will be
discussed in Chapters 6 and 7.
Equations (1.138) and (1.139) can be applied for a dispersion relation
not depending on time, when an ambient current is uniform in space. If
waves propagate on a non-stationary and non-uniform current u (x, t), the
intrinsic frequency σ may vary in space and time. In such a case, the
July 31, 2017 16:58 ws-book9x6 10666-main page 46

46 Introduction

spectral energy density will not be conserved. However Bretherton and

Garrett (1969) showed that the quantity called wave action density defined
as N = Ψ̂ (k; x, t)/σ is conserved in a moving medium. Therefore, instead
of Eq. (1.133), we apply a more general principle of conservation of wave
action density N in the form
by [Link] on 11/12/25. Re-use and distribution is strictly not permitted, except for Open Access articles.

∂N ∂Ω ∂N ∂Ω ∂N Q
+ − = . (1.140)
∂t ∂ki ∂xi ∂xi ∂ki σ
In order to illustrate the relationship between wave action density N and
energy density Ψ̂ let us consider the simplest example of a steady unidi-
rectional wave train in deep water on a variable current u(x) (source term
Ocean Surface Waves Downloaded from [Link]

Q ≡ 0). Then, Eq. (1.140) reduces to

  
∂ u + Cg
Ψ̂ = 0. (1.141)
∂x σ
For gravity waves in deep water σ = g/Cp and Cg = 1/2Cp . Thus
 
1
Ψ̂ u + Cp Cp = const. (1.142)
2
Although the flux of wave action is constant, the flux of wave energy is
not. For example, as Cp decreases in an adverse current, the energy flux
of the wave motion increases because of the work done by ‘radiation stress’
(Longuet-Higgins and Stewart, 1964).
The basic difficulty in solving Eq. (1.133) or (1.140) is an evaluation of
the function Q. In general, Q is a function of wavenumber k and a func-
tional of Ψ̂, i.e. Q = Q[k, Ψ̂]. Function Q can also be dependent on other
parameters which are functions of x and t (for example, wind velocity).
The theory of nonlinear weak interactions for wind-induced waves was
first formulated by Hasselmann (1962, 1968). Assuming that
1 ∂ Ψ̂ 1 ∂ Ψ̂
 1 and  1, (1.143)
k ∂ki ω ∂t
with the accuracy of the second approximation in interaction equations,
Hasselmann evaluated the general form of source/sink function Q as follows
(Hasselmann, 1968; Willebrand, 1975)
i=9
X
Q (k) = Qi , (1.144)
i=1

in which
Q1 = α, (1.145)
July 31, 2017 16:58 ws-book9x6 10666-main page 47

1.4. Methods of Description of Ocean Random Waves 47

Q2 = β Ψ̂ (k) , (1.146)
Z
Q3 = Ψ̂ (k) γ̂ (k, k0 ) Ψ̂ (k0 ) dk0 , (1.147)
by [Link] on 11/12/25. Re-use and distribution is strictly not permitted, except for Open Access articles.

Z
Q4 = −δ Ψ̂ (k) +  (k, k0 ) Ψ̂ (k0 ) dk0 , (1.148)

Z h
Q5 = T1 Ψ̂ (k0 ) Ψ̂ (k00 ) Ψ̂ (k − k0 − k00 )
i
− T2 Ψ̂ (k) Ψ̂ (k0 ) Ψ̂ (k00 ) dk0 dk00 ,
Ocean Surface Waves Downloaded from [Link]

(1.149)

Q6 = −νij ki kj Ψ̂ (k) , (1.150)

Q7 = −γ̂ Ψ̂ (k) . (1.151)

The coefficients α, β, δ, , γ̂ depend on the properties of the interacting fields.
Similar expressions for source/sink function Q in Eq. (1.139) are obtained
after appropriate variables transformation.
The terms Q1 and Q2 determine the wind input which may be added
or substracted from the momentum of existing waves. In particular, term
Q1 describes energy flux due to turbulent fluctuations of surface pressure
(Phillips’ mechanism). Term Q2 corresponds to Miles’ coupling mechanism
between the wave field and the mean boundary layer flow. The term Q3 is
a nonlinear correction of Miles’ theory and Q4 is responsible for interaction
of the wave motion with atmospheric turbulence. The term Q5 expresses
nonlinear interaction between spectral components when energy is moved
from one part of the spectrum to another. In general, nonlinear interactions
in the wave field are weak, as the evolution time of the wave field is much
greater than the characteristic wave period and the evolution in space is
much greater than the characteristic wavelength. Therefore, the process of
resonant wave-wave interaction may be conservative with no net loss or gain
of energy over the spectrum. However, in some regions and at some times,
the interactions may be quite strong. For example, strong interactions are
induced by wave breaking.
July 31, 2017 16:58 ws-book9x6 10666-main page 48

48 Introduction

The term Q6 represents the energy changes due to interaction of waves with
the mean current, while Q7 expresses the energy changes due to interaction
of waves with turbulence in the water layer. The two other components of
source function Q (Q8 and Q9 ) describe the strong nonlinearities. The term
Q8 is related to energy dissipation due to bottom friction and Q9 repre-
by [Link] on 11/12/25. Re-use and distribution is strictly not permitted, except for Open Access articles.

sents the dissipation of wave energy due to wave breaking. Wave breaking is
a strongly nonlinear interaction in which wave energy is moving from break-
ing waves to turbulence and to waves longer than waves corresponding to
the spectrum peak.
The most important components of function Q are the terms which
Ocean Surface Waves Downloaded from [Link]

describe the waves generation by wind: Qin = Q1 + Q2 , the nonlinear

interaction between wave components: Qnl = Q5 and energy dissipation
due to bottom friction and due to wave breaking: Qdiss = Q8 + Q9 . Further
discussion of the particular source Qi components, their importance for
wave prediction models and their last improvements is given in Chapters 6
and 7.

1.4.4 Representation of wave properties in time-frequency

space
[Link] Wavelet transform approach
Most of the oceanographic signals are the time-domain signals. However, in
many cases the most distinguished information is hidden in the frequency
spectrum which provides the energy associated with a given frequency. The
frequency spectrum of the signal can be obtained by the Fourier transform.
The Fourier transform yields an information on how much but not when
(in time) the particular frequency components exist. Such information is
sufficient in a case of the stationary signals as the frequency content of such
signals does not change in time and all frequency components exist all the
time. When the time localisation of the spectral components is required, the
transform of time series which provides the time-frequency representation
of the signal should be developed. A transform of such type is the Wavelet
transform which gives full time-frequency representation of the time series
(Torrence and Compo, 1998). In general, the Wavelet transform (WT) of
the signal, x(t), is defined as a following inner product

Z ∞
W T (τ, b) = hgτ b | xi = x(t)gτ∗b (t; τ, b)dt. (1.152)
−∞
July 31, 2017 16:58 ws-book9x6 10666-main page 49

1.4. Methods of Description of Ocean Random Waves 49

The family of continuously translated and dilated wavelets is generated

from mother wavelet g(t)
 
1 t−τ
gτ b (t; τ, b) = √ g , (1.153)
b b
by [Link] on 11/12/25. Re-use and distribution is strictly not permitted, except for Open Access articles.

where τ is the translation parameter, corresponding to the position of the

wavelet as it is shifted through the signal, b is the scale dilation parameter
determining the width of the wavelet and the (*) indicates the complex con-
jugate. The scale b > 1 dilates (or stretches out) the signals, whereas scale
b < 1 compresses the signal. The wavelet coefficients, W T (τ, b), represent
Ocean Surface Waves Downloaded from [Link]

the correlation (in term of the time-scale functions) between the wavelet
and a localised section of the signal. If the signal has a major component
of the frequency corresponding to the given scale, then the wavelet at this
scale is close to the signal at the particular location and the correspond-
ing wavelet transform coefficient, determined at this point, has a relatively
large value. Therefore, the Wavelet transform is a sort of a microscope with
magnification 1/b and location given by parameter τ , while the optics of
the microscope is characterised by the function gτ b (t; τ, b). For the wavelet
which has the mother wavelet status, the function g(t) must be regular
function such that it decays to zero in the limit |t| → ∞ and it has zero
mean.
The Wavelet transform should reflect the type of features which are
present in the time series. For time series with sharp steps, a boxcar-like
wavelet should be chosen, while for smoothly varying time series a smooth
function is more appropriate. However, if the wavelet power spectra are
not of the primary interest, the choice of wavelet function is not critical.
One of the most extensively used mother wavelet is the Morlet’s wavelet
 
1
g(t) = exp − t2 exp(ict). (1.154)
2
Eq. (1.154) represents a plane wave of frequency c, modulated by a Gaussian
envelope of the unit width. In the oceanographic applications, the other
mother wavelets, such as the orthogonal wavelets (Yamada and Ohkitani,
1990; Mori and Yasuda, 1994), Paul’s wavelet or DOG wavelet (deriva-
tive of a Gaussian) (Torrence and Compo, 1998) are also used. Using the
representation (1.153), the Morlet wavelet takes the form (Massel, 2001b)
"  2 #  
1 1 t−τ t−τ
gτ b (t) = √ exp − exp ic . (1.155)
b 2 b b
July 31, 2017 16:58 ws-book9x6 10666-main page 50

50 Introduction
by [Link] on 11/12/25. Re-use and distribution is strictly not permitted, except for Open Access articles.
Ocean Surface Waves Downloaded from [Link]

Fig. 1.8 The Morlet’s wavelet Eq. (1.155) of a constant value of c and different scales
b (a) b = 0.5, c = 2π; (b) b = 5.0, c = 2π (adapted from Massel, 2001b).

The frequency nature of the parameter c is clearly seen when we take c = 2π.
The Eq. (1.155) becomes
"  2 #  
1 1 t−τ 2π
gτ b (t) = √ exp − exp i (t − τ ) . (1.156)
b 2 b b
 
2π
The term exp i (t − τ ) represents now the plane sinusoidal wave with
b
a frequency 2π/b or with the scale dilation b which can be treated as a pe-
riod.
The real and imaginary parts of the wavelet (1.156) for various scale
parameters b are illustrated in Fig. 1.8. As the b increases, a width of the
Gaussian spreads in time from its centre value. If the scale parameter b is
constant and frequency parameter c increases, the number of oscillations
over the span of the function also increases. This effect is clearly seen by
comparison of Fig. 1.8 and Fig. 1.9.
The wavelet analysis is done in a such a way that the initial signal is
multiplied with a window (the wavelet). The window width is changed as
the transform is computed for every single spectral component. This is
the most significant characteristic of the Wavelet transform. Changing the
window width influences the resolution of the transform. It is schematically
illustrated schematically in Fig. 1.10. Every box corresponds to the value of
the Wavelet transform in the translation-scale plane. Although the widths
and heights of the boxes change, their areas are constant (for example, see
two hatched boxes). Each box represents the same portion of the time-
frequency plane, but gives different ratios of time and frequency. At low
frequency (high value of scale b), the height of the box is small, but its width
July 31, 2017 16:58 ws-book9x6 10666-main page 51

1.4. Methods of Description of Ocean Random Waves 51

by [Link] on 11/12/25. Re-use and distribution is strictly not permitted, except for Open Access articles.
Ocean Surface Waves Downloaded from [Link]

Fig. 1.9 The Morlet’s wavelet Eq. (1.155) of a constant scale b and different c parameters
(a) b = 5.0, c = π/2; (b) b = 5.0, c = π (adapted from Massel, 2001b).

is large, i.e. there is better frequency resolution and poor time resolution.
On the other hand, at higher frequencies (low scale b), the height of the
boxes increases and the width of the boxes [Link] corresponds to
a very narrow window when the frequency resolution is going poorer and
the time resolution is going better.

Fig. 1.10 Resolution scheme in the time-frequency plane (adapted from Massel, 2001b).
July 31, 2017 16:58 ws-book9x6 10666-main page 52

52 Introduction

Now we will examine some energy properties of the wavelet. Firstly, the
Wavelet transform conserves the energy, i.e.
Z ∞ Z ∞Z ∞
2 −1 2
|x(t)| dt = C | W T (τ, b) | b−2 dτ db, (1.157)
−∞ 0 0
by [Link] on 11/12/25. Re-use and distribution is strictly not permitted, except for Open Access articles.

in which
Z ∞  
−1 2
C = ω −1 |G(ω)| dω, (1.158)
−∞

and the G(ω) is the Fourier transform of function g(t).

Using the Wavelet transform and the coefficient C we can define various
Ocean Surface Waves Downloaded from [Link]

wavelet energy spectra and spectral densities. In particular, so called time-

scale energy density becomes
2
| W T (τ, b) |
E1 (τ, b) = . (1.159)
b
By integrating Eq. (1.159) versus scale b we obtain the local energy density
(Farge, 1992)
Z ∞
db
E2 (τ ) = C −1 E1 (τ, b) . (1.160)
0 b
On the other hand, the integration of Eq. (1.159) versus time τ gives the
global wavelet energy spectrum E3 (b) as
Z ∞
E3 (b) = E1 (τ, b)dτ. (1.161)
0

Torrence and Compo (1998) noted that the smoothed Fourier spectrum
approaches the global wavelet spectrum more closely when the amount of
necessary smoothing decreasing with increasing scale.
Finally, the total energy of the time series x(t) becomes
Z ∞
db
E = C −1 E3 (b) . (1.162)
0 b
After substituting Eqs. (1.159) and (1.161) into Eq. (1.162) we get
Z ∞Z ∞ Z ∞Z ∞
dτ db 2 dτ db
E = C −1 E1 (τ, b) = C −1 | W T (τ, b) | ,
0 −∞ b 0 −∞ b2
(1.163)
what constitutes the conservation of energy (1.157).
In Figs. 1.11 and 1.12, the contours of |W T | and global wavelet energy
spectrum E3 of the “chirp” signal are given. The frequency components,
July 31, 2017 16:58 ws-book9x6 10666-main page 53

1.4. Methods of Description of Ocean Random Waves 53

contributed to the signal, and time limits are as follows (Massel, 2001b)
  
2πt

 1.0 · sin for 5 s < t < 20 s,
2





  
 2πt
ζ(t) = 1.0 · sin for 20 s < t < 40 s, (1.164)
by [Link] on 11/12/25. Re-use and distribution is strictly not permitted, except for Open Access articles.


 5

  
2πt


 1.0 · sin for 40 s < t < 50 s.


10
Maximum values of the |W T | correspond to the scale parameters equiva-
lent to the periods of signal components. However, the maximum values
Ocean Surface Waves Downloaded from [Link]

of |W T | are shifting to the higher scales as the periods of the initial signal
components increase. Fig. 1.11 demonstrates that Wavelet transform is ca-
pable to provide the time and frequency information simultaneously, giving
a full time-frequency representation of the signal. In contrast to the Fourier
transform, the Wavelet transform allows exceptional localisation both, in
the time and in the scale (frequency) domains.

Fig. 1.11 Contours of the |W T | for the “chirp” type signal (adapted from Massel,
2001b).
July 31, 2017 16:58 ws-book9x6 10666-main page 54

54 Introduction
by [Link] on 11/12/25. Re-use and distribution is strictly not permitted, except for Open Access articles.
Ocean Surface Waves Downloaded from [Link]

Fig. 1.12 Global wavelet energy spectrum E3 for “chirp” signal (adapted from Massel,
2001b).

The global wavelet spectrum E3 (b) in Fig. 1.12 is plotted versus the quantity
f = 1/b which can be treated as a measure of frequency. The dominant
contribution into the global wavelet spectrum is provided by the signal
components with frequencies f = 0.1 Hz, 0.2 Hz and 0.5 Hz. The spectrum
E3 mimics the commonly used frequency spectrum. However it should be
pointed out that there is not a straightforward correspondence between E3
and frequency spectrum based on the Fourier transform.
Wavelets have been successfully implemented in signal and image pro-
cessing, ordinary and partial differential equations theory, numerical anal-
ysis and communication theory. On the other hand, the application of
the Wavelet transform to the fields of ocean engineering and oceanography
is not frequent. Meyers et al. (1993) demonstrated the usefulness of the
wavelet transform in studying dispersion of Yanai waves in a reduced grav-
ity equatorial model. In a paper by Shen and Mei (1994), the continuous
Wavelet transform was developed to analyse the energy balance in the equi-
librium spectral sub-range of wind-generated gravity waves. Scott (2005)
estimated the statistics of steep waves using the wavelet methodology to
obtain wave gauge data. Mori and Yasuda (1994), and Liu and Babanin
(2004) applied the Wavelet transform to detect wave growth and breaking
July 31, 2017 16:58 ws-book9x6 10666-main page 55

1.4. Methods of Description of Ocean Random Waves 55

in the time series. A sudden surface jump associated with breaking wave
is regarded as a shock wave and the wavelet spectrum is defined to de-
tect the occurrence of the surface jumps. The approach is shown to be
able to produce the same breaking statistics as field measurements of wave
breaking conditions based on the detection of whitecaps at fixed points of
by [Link] on 11/12/25. Re-use and distribution is strictly not permitted, except for Open Access articles.

observation.
Lin and Liu (2004) presented an application of the Wavelet transform for
the analysis of freak waves events (in-depth discussion of freak waves is given
in Chapter 8). The local wavelet energy density of the wave train is used
successfully to determine the groupiness degree of waves and its variation
Ocean Surface Waves Downloaded from [Link]

cross-shore (Liu, 2000; Dong et al., 2008), as well as for the separation
of wave fields into incident and reflected waves (Ma et al., 2010). Such
separation can be performed in the real time domain, also for non-stationary
wave trains.
Elsayed (2006a) and Dong et al. (2008) demonstrated that the Wavelet
transform has become a powerful tool for study of the nonlinear wave–wave
interactions of ocean waves. Let us first define the wavelet power spectrum
integrated over a finite time interval T
Z
Exx = W Tx∗ (τ, b)W Tx (τ, b)dτ, (1.165)
T

where W Tx (τ, b) is the wavelet transform of a given time series x(t).

In a similar way, the wavelet cross-spectrum takes the form
Z
Exy = W Tx∗ (τ, b)W Ty (τ, b)dτ, (1.166)
T

when W Ty (τ, b) is the wave transform of the second time series.

To the contrast to classical bispectrum (see Section 3.2.4), they defined
the wavelet based cross-bispectrum as a triple product of Wavelet trans-
forms as follows
Z
Byxx (b1 , b2 ) = W Ty∗ (τ, b) · W Tx (τ, b1 ) · W Tx (τ, b2 )dτ, (1.167)
T

where the scale dilation parameters b1 , b2 and b must satisfied the frequency
sum rule
1 1 1
= + . (1.168)
b b1 b2
The wavelet based cross-bispectrum is the measure of the phase coupling
in the interval T that occurs between wavelet components of time scales b1
and b2 of the time series x(t) and wavelet component b of time series y(t),
July 31, 2017 16:58 ws-book9x6 10666-main page 56

56 Introduction

satisfying Eq. (1.168). In practical applications, the following normalised

squared wavelet cross-bicoherence is used
2
2 | Byxx (b1 , b2 ) |
[byxx (b1 , b2 )] = R 2 R 2 .
[ T | W Tx (τ, b1 )W Tx (τ, b2 ) | dτ ][ T | W Ty (τ, b) | dτ ]
(1.169)
by [Link] on 11/12/25. Re-use and distribution is strictly not permitted, except for Open Access articles.

In particular, Elsayed (2006b) showed that phase coupling occurs between

wind speed and wave height over a certain range of frequencies during wave
growth in a Mistral event. However this range is different from one segment
to another due to the nonstationarity of the time series.
Ocean Surface Waves Downloaded from [Link]

[Link] The Hilbert transform representation of wave signal

The wavelet technique, based on Fourier spectral methodology with ad-
justable frequency-dependent window functions — the mother wavelets —
provides temporal-frequency resolution for non-stationary signals. How-
ever, the Fourier methodology interprets wave nonlinearity in terms of har-
monic generation with the spectral energy leaking into the higher-frequency
components. On the other hand, the Hilbert transform interprets wave non-
linearity as frequency modulation, and the spectral energy remains near the
base frequencies.10
To determine the function ξ(t), let us assume for a moment that a real-
valued function ζ(t) extends over the range −∞ < t < ∞. Then its Hilbert
transform becomes (Bendat and Piersol, 1986)
Z∞
ζ(t)
ξ(t) = H[ζ(t)] = P dt1 , (1.170)
π(t − t1
−∞

in which P denotes the Cauchy principle value. Eq. (1.170) indicates that
ξ(t) is the convolution integral of ζ(t) and 1/πt.
If the signal ζ(t) is represented in the form of a series, we can write
∞
X
ζ(t) = [an cos(nωt) + bn sin(nωt)] (1.171)
n=0
and
∞
X
ξ(t) = [an sin(nωt) − bn cos(nωt)]. (1.172)
n=0
10 The Hilbert transform is named after David Hilbert (1862–1943), a German mathe-

matician who developed many important tools used in modern mathematical physics.
For example, he introduced the extension of the real signal into the complex plane such
that it satisfies the Cauchy–Riemann equations (see Eq. (1.15)).
July 31, 2017 16:58 ws-book9x6 10666-main page 57

1.4. Methods of Description of Ocean Random Waves 57

The digital computation of the Hilbert transform was explained by Bendat

and Piersol (1986) and Massel (2007), where also was demonstrated that
the Hilbert transform of the velocity potential function appears to be the
stream function of the water motion (see also Eq. (1.15)). Here for our
applications we consider a record of N values (n = 0, 1, 2, ..., N − 1) of
by [Link] on 11/12/25. Re-use and distribution is strictly not permitted, except for Open Access articles.

function ζ(t) sampled at ∆t. Thus the Fourier transform Z(k∆f ) of ζ(t)
becomes
N −1  
X 2πkn
Z(k∆f ) = ∆t ζ(n∆t) exp −i , (1.173)
n=0
N
Ocean Surface Waves Downloaded from [Link]

where ∆f = 1/(N ∆t). The values of Z(k∆f ) are needed only from k = 0
up to k = (N/2), where the Nyguist frequency occurs. Now using the values
of Z(k∆f ) we can obtain the ξ(n∆t) values as follows
 
N/2  
X 2πkn 
ξ(n∆t) = 2∆f =  Z(k∆f ) exp i , (1.174)
N
k=0

in which = denotes the imaginary part of the given value.

Combining the real-valued function ζ(t) with its Hilbert transform ξ(t)
we obtain a complex analytical signal in the time domain
η(t) = ζ(t) + iξ(t) = A(t) exp(iΘ(t)), (1.175)
in which A(t) is interpreted as an instantaneous amplitude
p
A(t) = ζ 2 (t) + ξ 2 (t) (1.176)
and θ(t) is an instantaneous phase
 
ξ(t)
θ(t) = arctan . (1.177)
ζ(t)
From Eq. (1.175) the real part of the analytical signal is the given signal
ζ(t), i.e.
ζ(t) = <(η(t)) = A(t) cos(θ(t)). (1.178)
Therefore instantaneous frequency ω can be determined as
dθ(t)
ω= . (1.179)
dt
In contrast to the instantaneous envelope of a signal, the notion of an in-
stantaneous frequency is somewhat controversial. When applying Fourier
analysis, we need at least one full oscillation to define the local frequency
July 31, 2017 16:58 ws-book9x6 10666-main page 58

58 Introduction

value. Therefore, a process shorter than a full wave cannot be consid-

ered in a signal analysis. However, a straightforward application of the
Hilbert transform sometimes leads to the appearance of positive and neg-
ative frequency values. To overcome this difficulty, Huang et al. (1998)
introduced the empirical mode decomposition method. They found that at
by [Link] on 11/12/25. Re-use and distribution is strictly not permitted, except for Open Access articles.

any given time, there is only one frequency value which can only represent
one component. To satisfy this requirement, some limitations on the data
are necessary.
Probably the most obvious restriction is the assumption that the process
under the consideration is of narrow band, when the bandwidth parame-
Ocean Surface Waves Downloaded from [Link]

ter ν (see Eq. (3.9)) is much smaller than 1, and the expected numbers
of extremes and zero crossings have to be equal. However, the bandwidth
definition given by Eq. (3.9) has a global sense and it lacks precision. This
means that filtering the data within the bandwidth requirement does not
necessarily guarantee a meaningful instantaneous frequency. In order to
avoid such difficulties Huang et al. (1998) suggested applying the Gabor
(1946) statement that for any function to have a meaningful instantaneous
frequency, the real part of its Fourier transform must have only a positive
frequency. Subsequently Huang et al. (1998) modified this requirement by
imposing the restriction that the function should be symmetric with respect
to the local zero mean, and have the same numbers of zero crossings and
extremes. To define the local zero mean, the local mean of the envelopes
determined by the local symmetry is used. The first restriction is equiva-
lent to the traditional narrow band requirements for a stationary Gaussian
process. However, the second requirement in fact modifies the classical
global requirement to a local one in a such a way that the instantaneous
frequency will not have the unwanted fluctuations induced by asymmetric
wave forms.
A modification of the Hilbert transform, is known as the Hilbert–Huang
transform (HHT) — see Huang and Shen (2005). Pascoal et al. (2005) and
Veltcheva (2002) discussed the performance of the HHT and determined an
adequate stopping criterion in order to obtain physically meaningful results,
and to explain the cross-variation of the group characteristics.
The applications of the Hilbert transform methodology for the wave
signals in coastal zone is given in Chapters 9 and 10.