1.Introduction

1.1 Generals

TOMAWAC is a scientific software which models the changes, both in the time and in the spatial
domain, of the power spectrum of wind-driven waves and wave agitation for applications in the oceanic
domain, in the intracontinental seas as well as in the coastal zone. The model uses the finite elements
formalism for discretizing the sea domain; it is based on the computational subroutines of the TELEMAC
system as developed by the EDF R&D’s Laboratoire National d'Hydraulique et Environnement (LNHE).
The acronym TOMAWAC being adopted for naming the software was derived from the following English
denomination:

TELEMAC-based Operational Model Addressing Wave Action Computation

TOMAWAC is one of the models making up the TELEMAC system [Hervouet, 2007], which addresses the
various issues that are related to both free surface (either river- or sea-typed) and underground flows, as well
as the associated physical processes: bed-load transport, water quality, etc.

1.2 Implementing TOMAWAC

TOMAWAC models the sea states by solving the balance equation of the action density directional spectrum.
To serve that purpose, the model should reproduce the evolution of the action density directional spectrum at
each node of a spatial computational grid.
In TOMAWAC the wave directional spectrum is split into a finite number of propagation frequencies fi and
directions . The balance equation of wave action density is solved for each component (). The model
is said to be a third generation model (e.g. like the WAM model [WAMDI, 1988] [Komen et al., 1994]), since it
does not require any parameterization on the spectral or directional distribution of power (or action density).
Each component of the action density spectrum changes in time under the effects of the software-modelled processes.

1.3 TOMAWAC general purposes

TOMAWAC can be used for three types of applications:

1. Wave climate forecasting a few days ahead, from wind field forecasts. This real time type of application is rather directed to weather-forecasting institutes such as Météo-France, whose one mission consists of predicting continuously the weather developments and, as the case may be, publishing storm warnings.

2. Hindcasting of exceptional events having severely damaged maritime structures and for which field records are either incomplete or unavailable.

3. Study of wave climatology and maritime or coastal site features, through the application of various, medium or extreme, weather conditions in order to obtain the conditions necessary to carry out projects and studies (harbour constructions, morphodynamic coastal evolutions, …).

During the development of the TOMAWAC model, the LNHE laboratory has been interested mainly on
the last two types of applications. It considered also the possibility to carry out research activities focused on
the following topics:

wave-currents and wave-storm surge interactions, especially in those places where tide plays a significant role,

coastal morphodynamics,

probability of floods in coastal zone,

coastal structure stability and coast protection,

assimilation of wind or wave satellite data during computation…

2.Representing waves in TOMAWAC

2.1 General definition of waves

As stated in paragraph 1.1, the purpose of the TOMAWAC software consists of modelling the generation and
the spatio-temporal evolution of waves at the surface of the seas or of the oceans. Then, the main physical
process of interest is the wave or the sea states, these two terms being used interchangeably in this
document.

The word waves, generally means all the wind driven free surface waves propagating at the surface of the
ocean and the period of which (denoted as T) typically ranges from 2.5 to 25 s, or even, equivalently, whose
frequency f=1/T ranges from 0.04 to 0.4 Hz.

The sea state may take various forms, depending on whether the sea is still and quiet or, on the contrary, in
a stormy phase, whether the waves are being formed (the so-called wind sea) or, on the contrary, are
coming from the ocean after travelling several hundreds or thousand kilometres (the so-called “swell”).

2.2 Plane monochromatic waves

The most commonly used way to introduce wave modelling consists of considering simple sinusoidal waves
(they are often called regular waves). It is a monochromatic (one period or frequency) and plane (one
propagation direction) wave. The free surface elevation, which is denoted as , depends on the position (x,
y) of the point being considered in space, as well as on time t. It is written as:

is the wave amplitude (in meters) and corresponds to the distance from the wave crest
and the mean level at rest. The wave height, being measured from the crest to the trough
of the wave, is used as well: H=2a.

is the wave frequency (in rad/s). The period (in seconds) or the frequency f (in
hertz) = 1/T = is used as well.

is the wave number (in rad/m). The wavelength (in meters): is used as well. The
wave number k is yielded by the free surface wave linear dispersion relation, according to
frequency w and depth d:

is the wave propagation direction (in radians). Conventionally, this direction is measured
herein clockwise with respect to Y axis.

is the wave phase (in radians).

The energy per unit area of these progressive waves (which consists of kinetic energy and potential energy
in halves) amounts to:

wherein:

g is the gravity acceleration (g ≈ 9.81 ) is the water density (in ) ( for seawater).

2.3 Random multidirectional waves

A first representation of waves at the surface of the ocean is possible through the sinusoidal expression
being used in the preceding paragraph. When watching an actual sea state, however, not all the waves have
the same features, whether it is in terms of height, period or propagation direction. As a matter of fact, the
free surface wave energy is distributed over a range of frequencies (waves are then said to be irregular or
random) and over a range of propagation directions (waves are then called multidirectional). Mathematically,
that irregularity is expressed by writing that a real sea state results from the superposition of an infinite (or
large) number of elementary sinusoidal components (i.e. monochromatic and uni-directional components).
Thus, a random multidirectional wave field can be modelled through a superposition method, considering M
plane monochromatic components:

A major point in the above expression concerns the phase distribution of elementary wave components.
The approach used in the TOMAWAC model assumes that these phases are randomly distributed over the
[0; ] range with a uniform probability density. The various wave components are then independent, i.e. a
linear or phase averaged representation is used.
With the linear representation featuring TOMAWAC and using the random phase hypothesis, the energy per
unit area of random multidirectional waves can then be expressed as:

It is noteworthy, however, that the distortions of shallow water wave profiles cannot be modelled with such a
representation. This is because, as the water depth decreases, the non-linear processes linked to wave
propagation and wave interactions with the sea bottom get some importance. The waves become steeper
and dissymmetrical: they depart from a sinusoidal profile. A fine modelling of these non-linear effects
involves non-linear wave theories (3rd- or 5th-order Stokes waves, cnoidal waves, …) and/or so-called
phase resolving» propagation models modelling the evolution of each wave from a train, with a spatial
discretization of 20-50 points per wavelength (Boussinesq, Serre equations, …).

TOMAWAC is a phase averaged model: it is therefore a priori hardly suitable for modelling these non-linear
effects when the wave profile can no longer be considered as the superposition of a number of independent
sinusoidal components. In Section 4, however, it will be explained how the non-linear effects can be
processed and represented through source terms.

2.4 Sea state directional power spectrum

Real waves were introduced in the previous chapter as a discrete sum of elementary components. Actually,
the power spectrum over both frequencies and propagation directions is a continuous function. The relevant
variable for describing that sea state power spectrum is the directional spectrum of wave energy which is
also known as wave directional spectrum of energy and will henceforth be denoted as .

Correspondence with the discrete case of the previous section is set considering the following equivalence:

In case of a wave propagation in a zero-current medium, a balance equation of the wave energy directional
spectrum can be written taking into account some source and sink terms for energy generation or energy
dissipation.

2.5 Directional spectrum of sea state variance

The preferred variable for sea state representation and modelling is rather the variance density directional
spectrum.
This function, noted as F(f,q) and expressed in m2.Hz1.rad-1 is simply derived from the directional spectrum
of wave energy by the relation:

Then, in particular, we have:

The relation linking the variance density directional spectrum and the free surface elevation is then written in
the following pseudo-integral form:

It should be reminded that the phases are randomly distributed in that expression over the range [0 ; 2p] with
a uniform probability density. As regards the amplitude of each elementary component, it is related to the
variance density directional spectrum by:

Among these moments, the 0-order moment is equal to the variance of the free surface elevation:

In particular, that moment m0 affects the determination of the significant spectral wave height Hmo (equal to
the significant height H1/3 assuming that the wave heights are distributed according to a Rayleigh’s law) by
the relation:

The average frequencies and and are also used and computed as follows:

Further derived parameters can be computed from the variance density directional spectrum (see e.g. in
[AIRH, 1986]).

2.6 Sea state directional spectrum of wave action

In the general case of wave propagation in an unsteady medium (sea currents and/or levels varying in time
and space), the directional spectrum of the variance density is no longer kept and a new quantity should be
introduced, namely the directional spectrum of wave action.

That action density spectrum is related to the directional spectrum of variance density by the relation:

wherein denotes the relative or intrinsic angular frequency, i.e. the angular frequency being observed in a
coordinate system moving at the velocity of current. Such a frequency is different from the absolute angular
frequency w observed in a fixed system of coordinates. The two frequencies are linked by the Doppler effect
relation in the presence of a current

2.7 Selecting the directional spectrum discretization variables

The directional spectra of wave energy, variance or action shall generally be considered as functions
depending on five variables:

time t,

the pair of coordinates proving the spatial position of the point being considered. In TOMAWAC, these coordinates can be expressed either in a Cartesian coordinate system (x, y) or in a spherical coordinate system (latitude, longitude) according to the dimension of the computational domain.

the pair of variables applied for directional spectrum discretization, for which several solutions are theoretically possible:

The directional spectra output by TOMAWAC, however, are always expressed in . The
equations solved by TOMAWAC are thoroughly reviewed in section 4.

3. Hypotheses and application domain of TOMAWAC

3.1. Application domain of the model TOMAWAC

TOMAWAC is designed to be applied from the ocean domain up to the coastal zone. The limits of the
application range can be determined by the value of the relative depth d/L, wherein d denotes the water
height (in metres) and L denotes the wave length (in metres) corresponding to the peak spectral frequency
for irregular waves.

The application domain of TOMAWAC includes:

the oceanic domain, characterized by large water depths, i.e. by relative water depths of over 0.5 The dominant physical processes are: wind driven waves, whitecapping dissipation and non-linear quadruplet interactions.

the continental seas and the medium depths, characterized by a relative water depth ranging from 0.05 to 0.5. In addition to the above processes, the bottom friction, the shoaling (wave growth due to a bottom rise) and the effects of refraction due to the bathymetry and/or to the currents are to be taken into account.

The coastal domain, including shoals or near-shore areas (relative water depth lower than 0.05). For these shallow water areas, such physical processes as bottom friction, bathymetric breaking, non-linear triad interactions between waves should be included. Furthermore, it could be useful to take into account the effects related to unsteady sea level and currents due to the tide and/or to the weather-dependent surges.

Through a so-called finite element spatial discretization, one computational grid may include mesh cells
among which the ratio of the largest sizes to the smallest ones may reach or even exceed 100. That is why
TOMAWAC can be applied to a sea domain that is featured by highly variable relative water depths ; in
particular, the coastal areas can be finely represented.

The application domain of TOMAWAC does not include the harbour areas and, more generally, all those
cases in which the effects of reflection on structures and/or diffraction may not be ignored.

3.2. Wave interactions with other physical factors

Several factors are involved in the wave physics and interact to various extents with the waves changing
their characteristics. The following main factors should be mentioned:

The fine modelling of the interactions between these various physical factors and the waves is generally
rather complex and several research projects are currently focused on it. Within the application domain as
defined in the previous paragraph, TOMAWAC models the following interactions:

wave-bathymetry interaction: the submarine relief data input into TOMAWAC are constant in time, but the sea level can change in time. In addition to the effects of the sea level variations in time, TOMAWAC allows to take into account refraction, shoaling, bottom friction and bathymetric breaking.

wave-atmosphere interaction: this interaction is the driving phenomenon in the wave generation, takes part in energy dissipation processes (whitecapping, wave propagation against the wind…) and is involved in the energy transfer. To represent the unsteady behaviour of this interaction, TOMAWAC requires 10 m wind fields (specification of the couple of horizontal velocity components) with a time step matched to the weather conditions being modelled. These wind fields can be provided either by a meteorological model or from satellite measurements.

wave-current interaction: the sea currents (as generated either by the tide or by oceanic circulations) may significantly affect the waves according to their intensity. They modify the refractive wave propagation direction, they reduce or increase the wave height according to their propagation direction in relation to the waves and may influence the wave periods if exhibiting a marked unsteady behaviour. In TOMAWAC, the current field is provided by the couple of horizontal components of its average (or depth-integrated) velocity at the nodes of the computational grid. TOMAWAC allows to model the frequency changes caused either by the Doppler effect or by the unsteady currents, as well as by a non-homogeneous current field.

The physical processes modelled in TOMAWAC

Those interactions being taken into account by TOMAWAC have been reviewed and a number of physical
events or processes have been mentioned in the previous paragraph. These processes modify the total
wave energy as well as the directional spectrum distribution of that energy (i.e. the shape of the directional
spectrum of energy). So far, the numerical modelling of these various processes, although some of them are
now very well known, is not yet mature and keep on providing many investigation subjects. Considering the
brief review of physical interactions given in the previous paragraph, the following physical processes are
taken into account and digitally modelled in TOMAWAC:

—> Energy source/dissipation processes:

wind driven interactions with atmosphere. Those interactions imply the modelling of the wind energy input into the waves. It is the prevailing source term for the wave energy directional spectrum. The way that spectrum evolves primarily depends on wind velocity, direction, time of action and fetch (distance over which the wind is active). It must be pointed out that the energy which is dissipated when the wind blocks the waves is not taken into account in TOMAWAC.

whitecapping dissipation or wave breaking, due to an excessive wave steepness during wave generation and propagation.

These various processes are numerically modelled as presented in Part 4.

It should be remembered that, due to the hypothesis adopted in paragraph 3.1 about the TOMAWAC
application domain, the following physical processes are not addressed by the model (non-exhaustive list)

diffraction by a coastal structure (breakwater, pier, etc…) or a shoal, resulting in an energy transfer towards the shadow areas beyond the obstacles blocking the wave propagation.

reflection (partial or total) from a structure or a pronounced depth irregularity.

4. Mathematical modelling procedures used by TOMAWAC

4.1. Scope of sea state modelling

The directional spectrum of wave action density, as defined in paragraph 2.6, is considered as a function of
five variables:

using, as discretization variables:

the position vector for spatial location in a Cartesian coordinate system

the wave number vector for directional spectrum discretization, q denoting the wave propagation direction (direction in which the waves travel).

the time t.

Under the hypotheses made on the wave representation (see in paragraph 2.6) as well as on the
model application domain and the modelled physical processes (see in paragraph 3.3), an equation of
evolution of the directional spectrum of wave action can be written in the following form (see in
[Willebrand, 1975] [Phillips, 1977] [Bretherton, 1969] for a detailed demonstration of the way that equation is
arranged):

The equation expresses that, in the general case of waves propagating in a non-homogeneous, unsteady
environment (currents and/or sea levels varying in time and space), the wave action is preserved to within
the source and sink terms (designated by the term Q).

The following notation is also used in :

In that form (conservative writing in the form of a flux), equation (4.1) can be transposed to other coordinate
systems and, for instance, , or else can be used for the discretization of directional
spectrum [Komen et al., 1994] [Tolman, 1991].

Working in (x, y, kx, ky), however, makes it possible to remain in the canonical coordinate system and to
write, for the propagation equations (also named Hamilton’s equations):

wherein results from the Doppler relation applied to the wave dispersion relation for the general case with
current:

wherein:

is the absolute angular frequency observed in a fixed coordinate system.

is named absolute frequency.

denotes the current velocity (depth-integrated).

denotes the intrinsic or relative angular frequency, which is observed in a coordinate system
moving at the velocity . It is given by the dispersion relation in the zero-current case:

is named intrinsic or relative wave frequency.

d denotes the water height.

Through the Hamilton’s equations (4.2.a and 4.2.b), it can be demonstrated that we have:

or when defining :

The evolution equation can then alternatively be written in the following form (the so-called transport
form):

is the relative (or intrinsic) group velocity of waves, i.e. as is observed in a coordinate system moving at
the velocity of the current:

The relative (or intrinsic) phase velocity C of waves is also introduced :

The sea state spectral modelling will then consist of solving the evolution equations (4.1) or (4.6.a),
using the kinematic equations (4.7.a – 4.7.d).
The transport equation formulation (4.6.a) or (4.6.b) has been adopted in TOMAWAC, since it is closely
related to other equations applied in hydraulics, which have already been treated at the LNHE and for which
methods and a know-how have been developed long ago.

As regards the discretization variables being used in TOMAWAC, we have already mentioned in paragraph
2.7 that :

spatial discretization can be based either on a Cartesian coordinate system in (x, y) or on a spherical coordinate system at the Earth’s surface in (, ) = (longitude, latitude).

the x-axis (in the Cartesian coordinate system) or the -axis of longitudes (in the spherical

coordinate system) is assumed to be horizontal, directed to the right, whereas the y-axis (in
the Cartesian coordinate system) or the -axis of latitudes (in the spherical coordinate
system) is assumed to be vertical, upwardly directed. Then, in spherical coordinates, the
vertical axis points at the north, whereas the horizontal axis points to the East.

In either case, the wave propagation directions are defined with respect to the vertical axis

in the clockwise direction.

These conventions are illustrated below in Figure 4.1. Those equations that correspond to the two
spatial discretizations options are developed in the next paragraphs.

4.2. Equation solved

4.2.1. Equations solved in a Cartesian spatial coordinate system

By switching the variable from (x, y, kx, ky) to (x, y, , ), it can be shown that the following relation exists for
the directional spectrum of wave action as expressed in both coordinate systems :

The evolution equation (4.6.b) is then written as :

with the following transfer rates, as computed from the linear wave theory:

The operators and refer to the computation of a function gradient in directions that are respectively
normal and tangential to the characteristic curve with the direction :

Besides, using the dispersion relation (4.4), is can be demonstrated that :

The spatial transfer rates and (equations 4.12.a and 4.12.b) model the spatial wave propagation and
the shoaling. The directional transfer rate (equation 4.12.c) models the refraction-induced change of wave
propagation direction. Refraction is generated by the spatial variations of those properties of the environment
in which the waves propagate and can result either from a bathymetric variation (first term in 4.12.c) or from
current gradients (second term in 4.12.c). The relative frequency transfer rate (equation 4.12.d) models
the relative frequency changes resulting from sea level variations both in space and time and/or from current
variations in space.

It is noteworthy that this last term is zero in the case of zero-current and of no variation of sea level in time : the advection equation is then reduced to a three-dimensional equation.

Lastly, as regards the source terms, it should be mentioned that changing the coordinate system and using
the factor allows to switch from the term Q to a term that is directly expressed in terms of the
directional variance spectrum with a variance . The content of that term is explained in paragraph
4.2.3.

4.2.2. Equations solved in a spherical spatial coordinate system

By switching the variables from (x, y, kx, ky) to , it can be shown that the following relation exists
for the directional spectrum of wave action as expressed in both coordinate systems:

R denotes the Earth’s radius (R ≈ 6400 km) and, once more, and are respectively the longitude and the
latitude of the point being considered.

The evolution equation (4.6.b) is then written as :

with the following transfer rates:

As in the previous case, the operators nd refer to the computation of a function gradient in
directions that are respectively normal and tangential to the characteristic curve with the direction :

As previously, the spatial transfer rates and (equations 4.18.a and 4.18.b) model the wave propagation
in space and the shoaling. In that coordinate system, the directional transfer rate has an additional term.

(the first term in equation (4.18.c)) compared to the case in Cartesian coordinates. That term results from the
propagation in spherical coordinates, in such a way that waves are located with respect to the North change
during the propagation over a large circle arc at the Earth’s surface [WAMDI, 1988] [Komen et al., 1994].
Both second and third terms (equation 4.18.c) model the refraction caused respectively by bathymetry
and currents. The relative frequency transfer rate (equation 4.18.d) models the changes of relative
frequency resulting from variations of the sea level or of the current in both space and time. It is noteworthy
that this last term is zero in the case of zero current and of no variation of the sea level in time: the advection
equation is then reduced to a three-dimensional equation.

4.2.3. TOMAWAC source and sink terms

4.2.3.1. Generals

The source and sink terms that compose and in the right-hand members of evolution equations (4.11)
and (4.17) of directional spectrum of wave action gather the contributions from the physical processes listed
in paragraph 3.3. for the application domain of TOMAWAC:

These source and sink terms are numerically modelled and parameterized as detailed in the next
paragraphs. For most of these processes, several models or formulations are proposed and available in
TOMAWAC.

4.2.3.2. Wind input (term Qin)

Three wind generation models are available in TOMAWAC. The model to be activated is selected through
the keyword WIND GENERATION in the steering file, which can take four values, namely:

Beside those exponential growth-type wind generation models, a linear growth model is also available in
TOMAWAC, which has been proposed by Cavaleri & Malanotte-Rizzoli [Cavaleri and Malanotte-
Rizzoli, 1981] (see paragraph 4.2.3.2.4). The model can be activated through the keyword LINEAR WAVE
GROWTH, and can be used together with one of the three above mentioned models. Its main feature is that
it permits to start a wave simulation from a nil wave spectrum (whereas the three above mentioned models
need some initial energy level for the wave spectrum to grow under wind action).

4.2.3.2.1. Option 1 for wind input: Janssen’s model

With that option, the model implemented for the wind input term is based upon the Janssen’s works
[Janssen, 1989] [Janssen, 1991]; Janssen proposed a quasi-linear theory for modelling the ocean/atmosphere interactions. The linear growth term is ignored and only an exponential energy growth is
taken into account, following Miles’ results [Miles, 1957].

A quasi-linear source term is obtained as a function of the directional variance spectrum:

with the following notations:

is the ratio of air and water specific gravities ().

is the wave phase velocity

is the local wind direction (direction in which it blows)

is the friction velocity, being linked to the surface stress by the following relation:

is a constant allowing to offset the growth curve

The operator 'max' in the source term expression limits the wave generation for the propagation directions
included within angular sector with respect to the local wind direction . For the wave directions
making an angle in excess of 90° with respect to th e wind direction , the wind input term is void.
In the Janssen’s model [Janssen, 1991], the Miles’ parameter is a function of two non-dimensional variables:

the wave age:

the wind profile parameter:

It is written as

where

is the Von Karman’s constant

denotes a coefficient set to 1.2 by Janssen [Janssen, 1991].

denotes the roughness length

denotes the non-dimensional critical height:

The Janssen’s model [Janssen, 1989] [Janssen, 1991] is characterized by the method it uses for computing
and . The surface stress is a function depending, on the one hand, on the wind velocity U10 at 10 m
and, on the other hand, on the sea state roughness through the wave stress . It is obtained by solving the following system of equations:

The solution of the system of equations through a Newton-Raphson’s iterative method yields the surface
stress , the friction velocity and the roughness length .

The initial value of friction velocity being used in the iterative algorithm is obtained considering a
constant drag coefficient:

where: by default.

The wave stress itself is computed from the variance spectrum F (via the source term Qin) using the
following relation:

That integral is numerically computed over the discretized portion of the spectrum and a parametrization,
based upon a decrement of variance in , is used for the high frequencies portion of the spectrum.

In fact, that source term has eight parameters, namely:

coefficient (corresponding to the keyword WIND GENERATION COEFFICIENT in the steering file). Its default value is taken as 1.2, in accordance with the Janssen’s proposal [Janssen, 1991] and the value adopted in the model WAM-Cycle 4.

air specific gravity (corresponding to the keyword AIR DENSITY in the steering file. Its default value is taken as 1.225 kg/m3.

water specific gravity (corresponding to the keyword WATER DENSITY in the steering file). Its default value is taken as 1,000 kg/m3.

constant (corresponding to the keyword CHARNOCK CONSTANT in the steering file). Its default value is taken as 0.01, in accordance with the Janssen’s proposal [Janssen, 1991] and the standard value adopted in the model WAM-Cycle 4.

constant (corresponding to the keyword VON KARMAN CONSTANT in the steering file). Its default value is taken as 0.41, i.e. the typical value.

initial drag coefficient (corresponding to the keyword WIND DRAG COEFFICIENT in the steering file). This drag coefficient is provided for initializing the iterative computation of friction velocity u*. Its default value is taken as .

offset constant (corresponding to the keyword SHIFT GROWING CURVE DUE TO WIND in the steering file). Its default value is taken as 0.011, in accordance with the value adopted in the model WAM-Cycle 4.

elevation at which the wind is recorded (corresponding to the keyword WIND MEASUREMENTS LEVEL in the steering file). Its default value is taken as 10 m: its corresponds to the typical value and to the value being adopted in the above explanations.

4.2.3.2.2. Option 2 for wind input: Snyder et al. model

In that option, the model implemented for the wind input term is based upon the works conducted by Snyder
et al. [Snyder et al., 1981], as amended by Komen et al. [Komen et al., 1984] to take into account the friction
velocity u* instead of the wind velocity at 5 m. It corresponds to the formulation being used in the cycle 3
release of WAM model. The formulation is simpler than the Janssen’s theory which Option 1 is based upon
(see in preceding paragraph):

As in Option 1, the linear growth term is ignored and only an exponential energy growth is taken into
account, following the Miles’ results [Miles, 1957]:

The shear velocity value u* used is obtained considering a drag coefficient linearly depending on the wind
velocity:

where:

That source term only uses two parameters, namely:

air density (corresponding to the keyword AIR DENSITY» in the steering file. Its default value is taken as 1.225 kg/m3.

water density (corresponding to the keyword WATER DENSITY in the steering file). Its default value is taken as 1,000 kg/m3.

4.2.3.2.3. Option 3 for wind input: Yan's model

The Yan's model [Yan, 1987] consists of a combination of u*/C and (u*/C)2 terms. It is valid over a wide range
of frequencies and wind speeds:

To select this model, the keyword WIND GENERATION must be set to 3 in the steering file.

This source term makes use of four parameters. The default values of those parameters correspond to the
coefficients proposed by Westhuysen [Westhuysen et al., 2007].

The coefficient D, corresponding to the keyword YAN GENERATION COEFFICIENT D, has a default value of ;

The coefficient E, corresponding to the keyword YAN GENERATION COEFFICIENT E, has a default value of ;

The coefficient F, corresponding to the keyword YAN GENERATION COEFFICIENT F, has a default value of ;

The coefficient H, corresponding to the keyword YAN GENERATION COEFFICIENT H, has a default value of .

4.2.3.2.4. Linear wave growth: Cavaleri and Malanotte-Rizzoli model

The linear growth mechanism described by Phillips [Phillips, 1957], [Phillips, 1958] is useful to initialise wave growth. If this term is neglected, it is necessary to set a non-zero sea-state as initial condition in order
to enable the wave energy spectrum to grow.

The term that has been implemented in TOMAWAC is the linear wave growth term of Cavaleri & Malanotte-
Rizzoli [Cavaleri & Malanotte-Rizzoli, 1981], as formulated by Tolman [Tolman, 1992]:

where is the friction wind velocity, the wind direction and is a peak frequency called Pierson-
Moskowitz frequency [Pierson & Moskowitz, 1964], defined as:

To select this model, the keyword LINEAR WAVE GROWTH must be set to 1 in the steering file. This model
does not require any input parameter.

4.2.3.3. Whitecapping-induced dissipations (term Qds)

Two models are available in TOMAWAC. The whitecapping or the free surface slope-induced breaking is
activated through the keyword WHITE CAPPING DISSIPATION in the steering file; the keyword can take
three values, namely:

denotes the total variance, denotes the average intrinsic frequency and denotes the average wave
number; they are respectively computed as followings:

The formulas for computing the average frequency and the average wave number are derived from those in
use in WAM-cycle 4 [Komen et al., 1994]. These averages are not directly weighted by the variance
spectrum, since it was found, when WAM-cycle 3 [WAMDI, 1988] was being developed, that the above
expressions yielded more stable results than the conventional weighted averages. Lastly, it should be
pointed out that in TOMAWAC, the above average quantities are computed not only on the discretized
portion of the variance spectrum, but also analytically on the high frequency portion (up to +) considering a
decreasing variance in .

That source term has two parameters:

constant (corresponding to the keyword WHITE CAPPING DISSIPATION COEFFICIENT in the steering file). Its default value is taken as 4.5, in accordance with the proposal made by Komen et al. [Komen et al., 1984] and the standard value adopted in the model WAM-Cycle 4.

weighting parameter (corresponding to the keyword WHITE CAPPING WEIGHTING COEFFICIENT in the steering file). Its default value is taken as the 0.5 average value.

4.2.3.3.2. Option 2 for whitecapping: Westhuysen dissipation model

The Westhuysen dissipation model [Westhuysen et al., 2007] is based on a saturation-based model
formulation, which defines the term as depending on the saturation threshold .

The expression proposed by Westhuysen is:

where:

and

The variable w is set equal to 25.

This model is implemented in TOMAWAC in its most recent version, as formulated by Westhuysen
[Westhuysen, 2008], which combines the terms of Komen [Komen et al., 1984] (Qds
K) with that of Westhuysen [Westhuysen et al., 2007] as follows:

This model is selected by setting the keyword WHITE CAPPING DISSIPATION to 2 in the steering file.

This source term makes use of 4 parameters. Their default values correspond to the coefficients proposed
by Westhuysen [Westhuysen, 2008]:

The coefficient , corresponding to the keyword WESTHUYSEN DISSIPATION COEFFICIENT, has a default value of ,

The coefficient , corresponding to the keyword SATURATION THRESHOLD FOR THE DISSIPATION, has a default value of ,

The coefficient , corresponding to the keyword WESTHUYSEN WHITE CAPPING DISSIPATION, has a default value of 3.29,

The coefficient , corresponding to the keyword WESTHUYSEN WEIGHTING COEFFICIENT, has a default value of 0.0.

4.2.3.4. Bottom friction-induced dissipations (term Qbf)

A single model is available in TOMAWAC. The bottom friction-induced dissipation is activated through the
keyword BOTTOM FRICTION DISSIPATION in the steering file; the keyword can take two values, namely:

1. no bottom friction-induced dissipation (default value)

2. expression obtained during the JONSWAP campaign (Hasselmann et al. [Hasselmann et al., 1973]) and taken up by Bouws and Komen [Bouws, 1983].

The model used for the bottom friction-induced energy losses is an empirical expression globally
representing the various contributions from the wave-bottom interaction (percolation, friction…):

That (linear) expression is programmed in TOMAWAC using the following alternative form, which involves
the dispersion relation:

That source term has a single parameter:

constant (corresponding to the keyword BOTTOM FRICTION COEFFICIENT in the
steering file). Its default value is taken as , in accordance with what had been
obtained during the JONSWAP campaign [Hasselmann et al., 1973] and with the standard value
being used in the model WAM-Cycle 4.

4.2.3.5. Bathymetric breaking-induced dissipations (term Qbr)

In TOMAWAC, four parametric formulas are proposed for reproducing the effects of the bathymetric
breaking-induced energy dissipation on the spectrum. The bathymetric breaking-induced dissipation is
activated through the keyword DEPTH-INDUCED BREAKING DISSIPATION in the steering file; the keyword
can take five values:

0 No breaking-induced dissipation (default value)

1 Battjes and Janssen’s model [Battjes, 1978]

2 Thornton and Guza’s model [Thornton, 1983]

3 Roelvink’s model [Roelvink, 1993]

4 Izumiya and Horikawa’s model [Izumiya, 1984]

The first three models are parametric spectral models developed for studying the random waves,
whereas the fourth one is a turbulence model initially developed for studying the regular waves.
The general principle of the parametric spectral models consists in developing an expression for the total
dissipation of wave energy by combining a rate of breaker-induced dissipation with a breaking probability.

Whatever model is adopted, the directional spectrum version of the bathymetric breaking-induced dissipation term is based on the assumption that breaking does not affect the energy frequency and direction
distributions.

4.2.3.5.1.Battjes and Janssen’s model (1978)

The Battjes and Janssen’s breaking model [Battjes, 1978] is based on the analogy with a hydraulic jump.
Besides, it assumes that all the breaking waves have a height , which is of the same order of magnitude
as the water depth. The total energy dissipation term is expressed as follows

where denotes the maximum wave height being compatible with the water depth, is the fraction of
breaking wave, is a characteristic wave frequency and is a numerical constant of order 1.
can be computed either through the relation:

or through a relation derived from the Miche’s criterion

where is linked to by the linear wave dispersion relation.
is estimated, according to the Battjes and Janssen’s theory, as a solution of the implicit equation:

In TOMAWAC, that equation can be solved either in a dichotomous way or through explicit approximations
as proposed by Dingemans [Dingemans, 1983]. The latter are expressed as follows when putting:

version 1:

version 2:

version 2:

The directional spectrum version of the sink term is based on the assumption that breaking does not modify the frequency and directional distribution of energy. The source term Qbr is then written as:

Three constants can be modified using keywords:

constant corresponds to the keyword DEPTH-INDUCED BREAKING 1 (BJ) COEFFICIENT ALPHA in the steering file. Its default value is taken as 1, in accordance with the value as recommended by Battjes and Janssen [Battjes, 1978].

constant corresponds to the keyword DEPTH-INDUCED BREAKING 1 (BJ) COEFFICIENT GAMMA1 in the steering file. Its default value is taken as 0.88, in accordance with the value as recommended by Battjes and Janssen [Battjes, 1978].

constant corresponds to the keyword DEPTH-INDUCED BREAKING 1 (BJ) COEFFICIENT GAMMA2 in the steering file. Its default value is taken as 0.8, in accordance to the value as recommended by Battjes and Janssen [Battjes, 1978].

The following keywords are for selecting the model options:

The characteristic wave frequency is selected through the keyword DEPTH-INDUCED BREAKING 1 (BJ) CHARACTERISTIC FREQUENCY. Six values are possible:

1. average frequency: (refer to equation (4.29.b))

2. average frequency: , computed from the spectrum moments and (default value)

3. average frequency: , computed from the spectrum moments and

4. discrete peak frequency:

5. peak frequency computed through the Read’s method to order 5:

6. peak frequency computed through the Read’s method to order 8:

The computation mode for breaking probability (exact computation or utilization of an explicit approximation as proposed by Dingemans [Dingemans, 1983]) is selected through the keyword DEPTH-INDUCED BREAKING 1 (BJ) QB COMPUTATION METHOD. By default, version 2 of the explicit formulations as proposed by Dingemans [Dingemans, 1983] is used (see above). For applications, it is recommended not to modify the value of that keyword.

The computation mode for the maximum height compatible with the local water depth, , is selected through the keyword DEPTH-INDUCED BREAKING 1 (BJ) HM COMPUTATION METHOD. Two values are possible:

1. Relation: (default value)

2. Miches’ relation (see in (4.39) above)

4.2.3.5.2. Thornton and Guza’s model (1983)

The Thornton and Guza’s breaking model [Thornton, 1983] is based on the analogy with a hydraulic jump
and on two types of breaking wave height distribution. The energy sink term is written according to the
breaking wave height distribution being retained:

function 1:

function 2:

is the characteristic wave frequency (average frequency, , or peak frequency) and B is a
parameter ranging from 0.8 to 1.5 (its default value in TOMAWAC is B =1.0 ). The maximum wave height
compatible with the water depth, , is governed by the parameter .

The breaking model as proposed by Thornton and Guza can then be parameterized by the user via the
following 4 keywords:

5. peak frequency computed through the Read’s method to order 5: (default value)

6. peak frequency computed through the Read’s method to order 8:

DEPTH-INDUCED BREAKING 2 (TG) COEFFICIENT B, corresponding to the B variable. Its default value in the model is taken as 1.

« DEPTH-INDUCED BREAKING 2 (TG) COEFFICIENT GAMMA, corresponding to the g variable. Its default value in the model is taken as 0.42.

4.2.3.5.3.Roelvink’s model (1993)

The Roelvink’s breaking model [Roelvink, 1993] is based on the analogy with a hydraulic jump and on two
types of wave height distribution in the breaking zone (Weibull or Rayleigh). The energy sink term is written
according to the wave height distribution in the breaking zone:

Weibull’s distribution:

The coefficient is usually set to 0.65.

Rayleigh’s distribution:

denotes the characteristic wave frequency (average frequency, , or peak frequency), is a
numerical constant of order 1, is the proportional control factor between the allowable wave height and the
water depth (by default, ) and N is an exponent in the wake breaking weighting function (typically
N=10).

Thus, the Roelvink’s breaking model can be parameterized by the user via the following 5 keywords:

« DEPTH-INDUCED BREAKING 3 (RO) COEFFICIENT ALPHA, corresponding to the variable. Its default value in the model is taken as 1.0.

« DEPTH-INDUCED BREAKING 3 (RO) COEFFICIENT GAMMA, corresponding to the variable. Its default value in the model is taken as 0.54.

« DEPTH-INDUCED BREAKING 3 (RO) COEFFICIENT GAMMA2, corresponding to the variable. Its default value in the model is taken as 0.65.

« DEPTH-INDUCED BREAKING 3 (RO) WAVE HEIGHT DISTRIBUTION provided for retaining either a Weibull distribution (4.43) if the (default) value of the parameter is 1 or a Rayleigh distribution (4.45) if the parameter value is 2.

2. average frequency: , as computed from the spectrum moments m0 and m1

3. average frequency: , as computed from the spectrum moments m0 and m2

4. discrete peak frequency:

5. peak frequency as computed through the Read’s method to order 5: (default value)

6. peak frequency as computed through the Read’s method to order 8:

4.2.3.5.4.Izumiya and Horikawa’s turbulence model (1984)

Izumiya and Horikawa [Izumiya, 1984] sought an estimate of the dissipation by breaking-induced turbulence
in the case of regular waves. From the Reynolds’ equations and only considering a one-dimensional
condition, they obtained an expression of the breaking-induced dissipation of wave energy in the following
form:

E denotes the total wave energy, Cg and c are respectively group and phase velocities associated to the
characteristic wave frequency (average frequency , or peak frequency), is a parameter
governing the magnitude of the energy dissipation to be determined. For any profile, a shoal may induce the
wave reforming. In order to take that process into account, Izumiya and Horikawa express the factor in
terms of a deviation from a steady state:

where is a dimensionless quantity in the form:

From laboratory data, Izumiya and Horikawa set to 9 10-3 and to 1.8.

Assuming that the breaking does not affect the frequency and direction distribution of energy, the dissipation
term is lastly written as:

Thus, the breaking model as proposed by Izumiya and Horikawa can be parameterized by the user through
the three following keywords:

« DEPTH-INDUCED BREAKING 4 (IH) COEFFICIENT BETA0, corresponding to the variable. The default value in the model is 1.8.

« DEPTH-INDUCED BREAKING 4 (IH) COEFFICIENT M2STAR, corresponding to the variable. The default value in the model is 0.009.

2. average frequency: , as computed from the spectrum moments m0 and m1

3. average frequency: , as computed from the spectrum moments m0 and m2

4. discrete peak frequency:

5. peak frequency as computed through the Read’s method to order 5: (default value)

6. peak frequency as computed through the Read’s method to order 8:

4.2.3.6. Non-linear quadruplet interactions (term Qnl)

Three non-linear quadruplet interactions models are available in TOMAWAC. The non-linear quadruplet
interactions are activated through the keyword NON-LINEAR TRANSFERS BETWEEN FREQUENCIES in
the steering file; the keyword can take four values, namely:

4.2.3.6.1.Option 1 for non-linear quadruplet interactions: DIA method

The method and its implementation in TOMAWAC have been the subject of a specific report [Benoit, 1997]
which the reader is invited to refer to for further information. The major teachings of the DIA method are
summarized below.

The exact expression of the deep water interactions term as set by Hasselmann [Hasselmann, 1962]
[Hasselmann, 1962], expressed herein for convenience as a function of the wave number vector, is
analogous to a Boltzmann integral:

The energy exchanges, in that integral (a priori rather uneasily computable), take place between quadruplets
meeting the resonance conditions:

as evidenced by the two Dirac functions in the integral.

G denotes the value of the coupling term for the resonant quadruplet interactions

Establishing and computing its expression is also an awkward task. Hasselmann [Hasselmann, 1962]
proposed a computation mode that was also taken up and given a more concise form by such other authors
as Webb [Webb, 1978].

The exact computation of the above Boltzmann integral is too complex and time-consuming for such a sea
state operational model as TOMAWAC (see e.g. [Hasselmann, 1985]). That is why, starting from the
experiment as developed in WAM [WAMDI, 1988] [Komen et al., 1994], TOMAWAC uses the DIA (« Discrete
Interaction Approximation) approximate computation method as proposed by Hasselmann et al.
[Hasselmann et al., 1985]. Whereas the exact computation requires the summation of the contributions from
a great number of quadruplets, the approximate computation implements only a small number of quadruplet
configurations which are all similar.

That standard interaction configuration is defined as follows:

two of the wave numbers are alike: , which also involves that the two related frequencies are identical:

the other two frequencies s3 and s4 are defined by:

Through the value = 0.25, a good correlation with the exact computation of the integral [Hasselmann et al., 1985] could be achieved. That value is used in the model WAM [WAMDI, 1988] [Komen et al., 1994] and is taken up in TOMAWAC.

since the wave vectors and should observe the resonance condition, it can be shown they are featured by angles = 11.5° and = -33.6° with respect to the common direction of (refer to [Hasselmann et al., 1985]).

Furthermore, the mirror image is taken into account by considering the vectors as symmetrical with respect to the direction of .

The standard interaction configuration (in full line) and its mirror image (in dotted line) are shown
schematically in Figure 4.2.

With this standard configuration, the non-linear interaction term for all four resonant wave numbers is
written as [Hasselmann et al., 1985]:

With such a computation method, the vector \vec{k} scans all the discretization nodes of the directional spectrum
mesh. The number of configurations being considered is then twice as large as the number of points in that
mesh. In relation to the full computation, the 5-dimensional space (three integration dimensions and two
dimensions for \vec{k_4} of all the possible resonant quadruples is reduced to a 2-dimensional space.
In a finite water depth, from exact computations of the Boltzmann integral, Herterich and Hasselmann
[Herterich, 1980] suggested to make a deep water computation based on the previous method, then to
multiply it by a coefficient R representing the effect of the finite water height:

Coefficient R is a function of the normalized water height and is expressed as follows:

The average wave number was defined in the previous paragraph (see in (4.29.c)).
In its authors’ opinion, that relation remains valid as long as . It is used as such in TOMAWAC for the
finite water depth computations.

That source-term has a single parameter:

constant (corresponding to the keyword STANDARD CONFIGURATION PARAMETER of the steering file). Its default value is taken as 0.25, in accordance with the proposal made by Hasselmann et al. [Hasselmann et al., 1985] and with the standard value in the model WAMCycle 4.

The MDIA method (multiple DIA) is an extension of the DIA algorithm. We use here the version proposed by
Tolman [Tolman, 2004].

This method can give very reasonable results in simple situations, but in case of unsteady or rapidly
changing sea-state conditions (e.g. in the situation of abrupt changes of wind direction) it results in significant
qualitative and quantitative differences when compared with exact methods [Benoit, 2005].

The MDIA method consists of using a quadruplet depending on 2 parameters, l and m, defined as follows:

The and values proposed by Tolman that allow to best estimate the source term in the case of 4
interacting quadruplets are shown in Table 4.1.

To select this model in TOMAWAC the keyword NON-LINEAR TRANSFERS BETWEEN FREQUENCIES
must be set to 2 in the steering file.

This model does not require any other parameter: the values of the l and m parameters are set as constants
in the code. However they can be modified when considering a larger number of interacting quadruplets.

The Gaussian Quadrature Method (GQM) is based on the use of Gaussian quadratures for the different
numerical integrations arising in evaluating Equation 4.48. This technique, proposed by Lavrenov
[Lavrenov, 2001], has been developed and optimised to adequate results regarding both precision and CPU
time [Benoit, 2005], [Gagnaire-Renou, 2009], [Gagnaire-Renou et al., 2010]

Several steps are needed to transform Equation 4.48 into an expression that can be integrated via Gaussian
quadrature method. They can be summarized as follows (for a detailed description reference can be made to
[Gagnaire-Renou, 2009]):

1. Elimination of the Dirac function on the wave numbers , by imposing (see the resonance condition, Equation 4.49). Equation 4.50 is therefore reduced to an integral with 4 dimensions, including a single Dirac function on the frequency .

2. Variable change, to work with instead of , and reformulation of the equation in terms of variance density (F) instead of wave action density (N).

3. Integration over and elimination of the Dirac function on the frequency. A 3-dimension integral is obtained, without any Dirac function.

4. Final expression of the non-linear transfer term: the variables , and , are expressed as functions of , and , and of et , which term depends on. The variables , et , defined by

and depending only on ( et ) and ( et ), are used as well.

where is defined as .

Equation 4.55 is then integrated using different quadrature methods:

Gauss-Legendre or Gauss-Chebyshev quadratures are used for the integration over s2, depending on the ea values determining number and type of singularities.

Gauss-Chebyshev quadratures are used for the integration over q1.

The integration over s1 is realized using the trapezoidal rule.

Three different GQM method resolutions have been tested:

A “fine” resolution, considered as the exact calculation of the non-linear transfer term, as no improvement in the results is noticed when further increasing the method resolution.

An “intermediate” resolution.

A “coarse” resolution, whose parameters are given as default values in TOMAWAC, which represents the best compromise between accuracy of the solution and CPU time.

The configurations that do not effect significantly the global computation of Qnl4 are neglected. This configuration selection allows to reduce the CPU time. The threshold values set as default in TOMAWAC reduce the number of configuration:

by 21% in the “fine” resolution case

by 34% in the “intermediate” resolution case

by 64% in the “coarse” resolution case

The GQM method is selected by setting to 3 the keyword NON-LINEAR TRANSFERS BETWEEN
FREQUENCIES in the steering file.

This method makes use of 6 parameters. The default values of those parameters correspond to the “coarse”
resolution case:

The three keywords SETTING FOR INTEGRATION ON OMEGA1, SETTING FOR INTEGRATION ON OMEGA2 and SETTING FOR INTEGRATION ON THETA1 determine the number of integration points over the three variables s1, q1 and s2 and their default values are respectively 3, 3 and 6 (“coarse” resolution). The values 1, 4, 8 and 2, 8, 12 correspond respectively to the “intermediate” and “fine” resolution cases.

The three keywords THRESHOLD0 FOR CONFIGURATIONS ELIMINATION, THRESHOLD1 FOR CONFIGURATIONS ELIMINATION and THRESHOLD2 FOR CONFIGURATIONS ELIMINATION affect the percentage of discarded configurations. Their default values are respectively 0, 1010 and 0.15. For the “intermediate” and “fine” resolution cases, the first two values are the same, and the threshold2 values is equal respectively to 0.01 and 0.001.

4.2.3.7. Non-linear transfers between triads (Qtr term)

4.2.3.7.1.LTA (Lumped Triad Approximation) model

A parametric model allowing to take into account the non-linear triad interactions in the averaged-phase
models has been proposed by Eldeberky and Battjes [Eldeberky, 1995]. The LTA model is a parametric
approach that is based on the Madsen and Sorensen’s deterministic spectral model [Madsen, 1993].
Simplifying hypotheses are introduced for reducing the computation cost. Thus, a parametric formulation is given for the biphase as a function of the Ursell’s number and the model is restricted to the self-interactions.
The source term is written as:

is the model adjustment coefficient; c and Cg denote the phase and group velocities, respectively.

R is the self-interaction coefficient

The biphase is given by the relation

where Ur denotes the Ursell's number

with being the significant spectral height and being the average wave time.
denotes the negative and positive contributions of the self-interactions. Since denotes the
positive contributions to the first upper harmonic, it should be positive. The negative values of are
replaced by the zero value. In the numerical integration of the energy equation, the source term for the triad
interactions is generally only computed for frequencies that are lower than Rfmfm (Ris [Ris, 1997]) where .

4.2.3.7.2. SPB model

The SPB model was developed by Becq [Becq, 1998] from the extended Boussinesq equations as
proposed by Madsen and Sorensen [Madsen, 1992]. The model is for simulating the effects induced by the
collinear and non-collinear interactions of spectral components. The source term is written as:

F denotes the variance spectrum in terms of frequencies and directions, . and respectively correspond to the sum and difference interactions. K is the model adjustment parameter.

Since the model was designed for taking into account the energy transfers for all the possible triad
configurations within the spectrum, the computation times are very long. In order to shorten these
computation times, the interactions can be restricted over a range of spectral components that are included
within a given angular sector. Thus, directional limits can be user-prescribed.

5. Discretizations used in TOMAWAC

The main aspects concerning the numerical discretization in TOMAWAC are presented and discussed herein
for the two spatial variables (paragraph 5.1), for the two spectro-angular variables (paragraph 5.2) and for
the time domain (paragraph 5.3).

5.1. Spatial discretization

The spatial coordinate system, whether it is Cartesian or spherical, is a planar two-dimensional domain that
is meshed by means of triangular finite elements. Only the maritime portion of the computational domain is
meshed, so that all the computational points of the spatial grid are provided with a water depth that is strictly
above zero. Through this discretization technique, the mesh size may naturally be variable over the spatial
domain, particularly enabling to get a fine grid in the areas of specific interest, featured either by complex
geometries (straits, intracontinental seas, bays…) or by high bathymetric gradients. Furthermore, that spatial
grid may include one or more islands.

The number of discretization points is only limited by the RAM capacities of the computing machine. The
equation solved by TOMAWAC does not prescribe a priori any conditions about the number of grid points per
wave length. The density of spatial discretization points is left at the user’s will. It should match, however,
both spatial and temporal scales of variation of the physical characteristics of the domain being studied, in
particular bathymetry and wind field.

In the general case, this spatial grid is realised on a workstation using one of the mesh generators
associated to the TELEMAC system (refer to the 7.2.2 for further details about the preparation of the grid).
Two examples of spatial grids developed for TOMAWAC for simulated storms in the North Atlantic Ocean,
the Channel and the North Sea are illustrated in Figure 5.1.

5.2. Spectro-angular discretization

5.2.1. Frequency discretization

In TOMAWAC, the frequency domain is discretized considering a series of NF frequencies in a geometric
progression:

with n ranging from 1 to NF

The minimum frequency is then f1 and the maximum frequency is .

In order to define the frequency discretization, the user should specify as an input into the steering file:

the frequency number: NF (corresponding to the keyword NUMBER OF FREQUENCIES in the steering file)

the minimum frequency: (in Hertz) (corresponding to the keyword MINIMAL FREQUENCY in the steering file)

the frequential ratio: q (corresponding to the keyword FREQUENTIAL RATIO in the steering file)

5.2.2. Directional discretization

The interval of propagation direction [0, 360°] is discretized into ND evenly distributed directions, so that
these directions are:

with m ranging from 1 to ND

In order to define the directional discretization, the user should specify as an input into the steering file:

the direction number: ND (corresponding to the keyword NUMBER OF DIRECTIONS in the steering file).

The direction convention selected for the input/output directional variables: either nautical or counterclockwise (corresponding to the keyword TRIGONOMETRICAL CONVENTION in the steering file, the default value of which is NO). The nautical convention sets the wave propagation directions (towards which the waves are propagating) in relation to the true North or the vertical axis and opposite to the counterclockwise direction. The counterclockwise convention sets the wave propagation directions in relation to the horizontal axis.

Note that the convention selected for computing the directions within the FORTRAN model always defines the propagation directions in the clockwise direction from the true North, even though the keyword TRIGONOMETRICAL CONVENTION = YES !

5.2.3. Spectro-angular grid

A two-dimensional grid for spectro-angular discretization is achieved by combining the above defined
frequency and directional discretizations. That grid has NF.ND points.

A polar representation is used in TOMAWAC, where the wave frequencies are measured radially and where
the propagation direction corresponds to the value of the angle in relation to the axis selected by the user as
(vertical or horizontal) origin. An example of a spectro-angular grid having 25 frequencies and 12 directions
is illustrated in Figure 5.2.

5.3. Temporal discretization

In TOMAWAC, each computation begins at the internal date 0, to which an actual date being defined by the
keyword DATE OF COMPUTATION BEGINNING in the steering file can be associated. That date is
specified as per the yymmddhhmm format which corresponds to the moment dd/mm/yy at hh:mm (for
example, 9505120345 corresponds to May 12, 1995 at 3.45).

The evolution equation of the directional spectrum of wave action density is integrated with a constant time
step which is expressed in seconds through the keyword TIME STEP in the steering file. Sub-iterations of
that time step can also be made for computing the source terms (refer to paragraph 6.3). That number of
sub-time steps per time step is defined in the steering file through the keyword NUMBER OF ITERATIONS
FOR THE SOURCE TERMS.

6. Numerical methods used in TOMAWAC

6.1. General solution algorithm

As stated in Section 4, the equation to be solved by the TOMAWAC software is a transport (convection)
equation with source terms that can be written in the following general form:

Both functions F and Q are functions of five variables and depend, e.g. in Cartesian coordinates, on x, y, ,
and t. The above equation is then to be solved on a four-dimensional grid in (x, y, , ) and is a transport
vector, which is a dimension-4 vector in the general case. It is reduced, however, to a three-dimensional
vector ( is zero) when there is neither a current nor a variation of water depth in time.

Equation (6.1) is solved in TOMAWAC through a fractional step method, i.e. the convection and the source
term integration steps are solved successively and separately. Thus, the following steps are successively
solved from a current state at the date , in which the variance spectrum is known in all points:

a convection step without source terms (refer to paragraph 6.1.2):

discretized as follows:

from which a value of , then of , intermediate after the convection step, is derived

a source term integration step (refer to paragraph 6.1.3):

discretized as follows:

since coefficient B is time independent.

The variance density spectrum for a time step (time ) is then obtained. That operation is
then repeated for the next time step and as many times as necessary for covering the simulation period
being considered.

6.2. Processing the propagation step

The propagation step is solved in TOMAWAC by means of the method of characteristics which is used at the LNHE for processing various convection equations (refer for example to [Esposito, 1981]). The
application of that method to TOMAWAC has a specific feature: the method should be applied to a
dimension-4 space in the general case and to a dimension-3 space when there is no current and the depth is
constant over time; furthermore the domain in propagation directions is periodic.

It should be reminded that equation (6.3) without source terms is processed in that step, being discretized as
follows for a time step :

The convector field , whose expression was given in Section 4, is not time dependent when there is no
tide, just like factor B (refer to paragraph 4). The equation to be processed can be simplified as follows in that
case:

This is a major advantage, since the characteristics can be traced back only once, at the beginning of the
simulation. It is sufficient to store the origin of the characteristic pathlines and to retrieve them whenever the
convection step is called. For each quadruplet of the discretized spatial and spectro-angular
variables, the characteristic curve is traced back to the time step and the “arrival” point ,
which is called foot of the characteristic pathline, is stored. Actually, the numbers of the discretization
elements (triangular elements for the spatial grid and quadrangular elements for the spectro-angular grid)
including that foot of the characteristic pathline, as well as the linear interpolation coefficients allowing to
obtain the values in that point from the values at the apices of the elements (barycentric coordinates), are
kept. Thus, the convection step can be reduced in the form:

That step requires a short computation time since it only consists of an interpolation operation over
each time step, once the characteristics have been traced back at the beginning of a computation.
When there is a tide, the principle remains unchanged, but the characteristics should be traced back after
every depth and current update.

Such a method has the advantage of being unconditionally stable, enabling to revoke the condition that
requires a Courant number below 1 and which is implemented, for example, in the upstream off-centred firstorder
propagation scheme being used in the WAM-cycle 4 model [WAMDI, 1988] [Komen et al., 1994]. The
finite element grid generation technique is provided for achieving a locally finer computational grid in order to
represent irregular bathymetric features or an irregular coastline. Thanks to the applied propagation scheme,
the time step does not necessarily have to be much reduced, so that reasonable computation times can be
kept. It should actually be clear that, rather than the propagation step, the source term integration step
(particularly the computation of non-linear interactions) does consume most of the computation time. As
regards the numerical schemes in which the propagation step implies a shorter time step when making the
grid finer (e.g. as in the case of the WAM model), the overall computation time happens to become much
longer because of the source terms and the model becomes less attractive for the practical applications.
Owing to the method of characteristic, on the contrary, the TOMAWAC model allows to overcome that
restriction and is therefore attractive even for grids with a rather fine spatial resolution.

The method of characteristic, however, has some drawbacks due to the fact that, in the general case, it has
a significant level of numerical diffusion and is not conservative.

6.3. Processing the source term integration step

6.3.1. Source term integration numerical scheme

The source and sink terms in the equation of variance density spectrum evolution are integrated using semi-implicit scheme:

where the exponent * denotes the values of the variables after the propagation step (but before the source
term integration step) and the exponent n+1 denotes the values of the variables after the source term
integration step. Emphasis should be laid on the fact that the source term integration step is local, i.e. it is
carried out independently for each point in the 2D spatial grid.

That scheme is inspired by the scheme that is used in the WAM-Cycle 4 model [WAMDI, 1988] [Komen et
al., 1994] since it enables to use fairly long time steps (about 20-30 min in an oceanic environment).
It is defined: and the source and sink terms are classified as linear and non-linear terms in F:

Q = Qnl + Ql

As regards the source terms that are linear in F, note that: , hence we get:

As regards the source terms that are non-linear in F, a Taylor’s expansion is made keeping only the first-order term:

is a matrix of differential increments that is broken down into a diagonal portion [L*] and an extradiagonal portion :

substituting into the expression of , this yields:

Adding the contributions from the linear and nonlinear terms, we obtain:

The variation of the variance density spectrum due to the source terms is written as:

i.e., after substitution of the source term expressions:

The matrix between brackets in the left-hand member of the latter equation cannot be easily inverted in the
general case. The designers of the WAM model [WAMDI, 1988], however, demonstrated that the diagonal
portion [L*] usually prevails over the extra-diagonal portion [N*]. Relying on comparative tests, they conclude
that the extra-diagonal portion can be ignored in favour of the diagonal portion, even with time steps of 20
min. or so. Due to that simplification, the inversion is much easier and we finally obtain:

For the sake of convenience, that expression is rewritten as:

where :

denotes a total source term, and

denotes a source term derived with respect to F.

The contributions of the various source terms implemented in TOMAWAC and described in paragraph 4.2.3
are schematically illustrated in the Table 6.1.

The source term integration time step may be different from the propagation time step in TOMAWAC, but it
should be a sub-multiple of it. Thus, several source term integration time sub-steps per propagation time step
can be defined. That option is governed by the keyword NUMBER OF ITERATIONS FOR THE SOURCE
TERMS in the steering file. The default value of that parameter is set to 1.

It was found experimentally that the depth-induced breaking source term, which is sometimes very strong,
can still be overestimated if the time step that was selected for the source term integration is too long. In
order to avoid that, TOMAWAC gives an opportunity to make a number of time sub-steps that are specific to
that source term. These time sub-steps are in a geometric progression. In order to limit that number of time
sub-step, TOMAWAC first clips the wave height by setting a maximum Hs/d ratio (d being the depth) to 1.

Subsequently, a Euler’s explicit scheme is used at each time step:

The Hs/d ratio can be modified through the keyword MAXIMUM VALUE OF THE RATIO HM0 ON D
(however, this it not advisable). The number of time sub-steps is specified through the keyword NUMBER OF
BREAKING TIME STEPS. The geometric ratio is given by the keyword COEFFICIENT OF THE TIME SUBINCREMENTS
FOR BREAKING.

6.3.2. Monitoring the growth of the wave spectrum

In order to limit the possible risks of numerical instabilities related to the source term integration, TOMAWAC
is provided with a criterion for limiting the growth of the directional spectrum per source term integration time
step. That criterion is directly inspired by the criterion proposed by the WAM group [WAMDI, 1988].

The absolute variation of the variance density spectrum as it was computed by the semi-implicit scheme in
paragraph 6.1.3.1. should remain lower than a fraction of an equilibrium spectrum :

Such an expression of a growth limiter is not always valid. Recent tests demonstrated that it might be
disadvantageous for wave growth over short fetches [Herbach et al., 1996].

6.4. Processing the boundaries – Boundary conditions

6.4.1. Spatial grid

Two types of boundary conditions are considered in TOMAWAC for the finite element spatial grid:

The former corresponds to a free boundary condition, i.e. that absorbs the whole wave energy. It may be a sea boundary, hence it is assumed that the waves propagate beyond the domain and nothing enters it. It may be a solid boundary, hence it is assumed that the coast absorbs completely the wave energy (no reflection).

The latter corresponds to a prescribed value boundary condition. The whole wave spectrum is then prescribed at each point along that boundary and for each step. Energy enters into the computational domain.

6.4.2. Spectro-angular grid

As regards the propagation directions, the grid generation is periodical over the range [0 ; 360°]: he nce there
are no directional boundary conditions.

As regards the wave frequencies that are discretized, the minimum and maximum frequency markers are
considered as “open boundary limits”, where the energy can be transferred to lower or higher frequencies,
exiting the discretized frequency range.

7. Inputs-outputs

7.1. Preliminary remark

During a computation, the TOMAWAC software uses a number of files, some of which are optional, as inputs
and outputs.

The input files are:

The steering or CAS file (mandatory),

The mesh or geometry file (mandatory),

The boundary conditions or CONLIM file (mandatory),

The seabed, bottom or bathymetry file (optional),

The FORTRAN or PRINCI file (optional),

The currents file (optional),

The winds file (optional),

The previous computation file (optional),

The binary user file (optional),

The formatted user file (optional).

The output files are:

The 2D results or grid file (mandatory),

The punctual results or spectra file (mandatory),

The next computation file (optional),

The listing printout (either on the display screen or in the file, see in Appendix 2),

The binary user file (optional),

The formatted user file (optional).

7.2. Preliminary remark

7.2.1. The steering (or CAS) file

The steering file name is specified in the steering file through the keyword: STEERING FILE.

It is a text file created by means of a text editor. In a way, it serves as the computation control panel. It
includes a set of keywords to which values are assigned. If a keyword does not appear in this file, then
TOMAWAC will assign to it the default value as defined in the dictionary file (refer to the description in
APPENDIX 3). If such a default value is not defined in the dictionary, then the computation will come to a halt
and display an error message. For instance, the command NUMBER OF DIRECTIONS = 12 is for specifying
that the direction spectrum of wave action or its moments will be discretised over 12 propagation directions.

TOMAWAC reads the steering file at the beginning of the computation.

Both dictionary file and steering file are read by the so-called DAMOCLES utility which is included in
TOMAWAC. The syntactic rules of DAMOCLES should then be observed upon the creation of the steering
file These rules are described here below.

The write rules are as follows:

The keywords can be of the Integer, Real, Logical or Character format type.

The keyword sequence order in the steering file is of no importance.

Each line has a maximum of 72 characters. However, as many linefeeds as one wants are allowed provided that the keyword name does not run from one line to the next.

For the table-like keywords, the successive values are separated by a semi-colon. A number of values equal to the table dimension should not necessarily be given; in such a case,DAMOCLES returns the number of values being read. For example:

ABSCISSAE OF SPECTRUM PRINTOUT POINTS = 1.2;3.4

(that keyword is declared as a 19-valued table)

The symbols ”:” or ”=” are indiscriminately used to separate a keyword from its value. They can be either preceded or followed with any number of blanks. The value itself may appear on the next line. For example:

NUMBER OF DIRECTIONS = 12

or

NUMBER OF DIRECTIONS: 12

or else

NUMBER OF DIRECTIONS =

12

The characters occurring between a pair of ”/” on one line are regarded as comments. Likewise, the characters occurring between a ”/” and a the end of a line are regarded as comments. For example:

TYPE OF BOUNDARY DIRECTIONAL SPECTRUM = 1 / Jonswap spectrum

A whole line beginning with a ”/” in the first column is regarded as a comment, even though another ¨/¨ occurs on the line. For example: / The geometry file is ./maillage/geo

Integer writing: Do not exceed the maximum size being allowed by the machine (in a machine with 32 bit architecture, the values range from -2 147 483 647 to + 2 147 483 648. Do not enter a blank between the sign (optional for the + sign) and the number. A dot at the end of the number is tolerated.

Real writing: A dot or a comma is allowed as a decimal point, as well as the FORTRAN E and D formats (1.E-3 0.001 0,001 1.D-3 denote the same value).

Logical value writing: The values 1, YES, OUI, .TRUE., TRUE, VRAI on the one hand, and 0, NON, NO, .FALSE., FALSE, FAUX on the other hand are allowed.

Character string writing: Those strings including blanks or reserved symbols (”/”,”:”, ”=”, ”&”) should be put in single quotes ('). The value of a character keyword may include up to 144 characters. As in FORTRAN, the quotes occurring within a string should be doubled. A string may neither begin nor end with a blank. For example:

TITLE = 'HOULE D''OUEST'

In addition to the keywords, a number of directives or metacommands that are interpreted during the
sequential readout of the steering file may be used as well:

The &FIN command indicates the end of file (even though the file is not completed). Thus, some keywords can be disabled simply by placing them behind that command for easily making it possible to enable them again subsequently.

The &ETA command prints the list of keywords and the relevant values at the time when DAMOCLES meets that command. This display will occur at the beginning of listing printout.

The &LIS command prints the list of keywords. This display will occur at the beginning of listing printout.

The &IND command prints the detailed list of keywords. This display will occur at the beginning of listing printout.

The &STO command causes the interruption of the program, the computation does not go on.

7.2.2. The geometry file

The geometry file name is specified in the steering file through the keyword: GEOMETRY FILE.

It is a SERAFIN-formatted binary file: it can be read by FUDAA PRE-PRO or RUBENS and it can be created
by the STBTEL module from the file(s) as produced by the mesh generator. The SERAFIN format structure
is described in APPENDIX 9.

This file includes the complete information about the horizontal mesh, i.e. the number of mesh points
(variable NPOIN2), the number of elements (variable NELEM2), the X and Y vectors containing the coordinates
of all the points and, lastly, the IKLE2 vector containing the connectivity table.

Furthermore, this file may also include bathymetry information in each point of the mesh, provided that the
interpolation of the bathymetry was carried out during the execution of the STBTEL module or during the
generation of the mesh.

TOMAWAC reproduces the information regarding the geometry at the beginning of the 2D results. Any
computation results file can then be used as a geometry file when one wants to perform a further simulation
on the same mesh.

7.2.3. The boundary conditions file

The boundary conditions file name is specified in the steering file through the keyword: BOUNDARY
CONDITIONS FILE.

It is a formatted file that can be created automatically by STBTEL and can be modified by means of a text
editor. Each line in this file is assigned to one point of the boundary and listed in sequential order in terms of
the boundary node numbers. The numbering of the boundary points first delineates the domain contour in
the counterclockwise direction, then the islands in the clockwise direction.

This file is described in detail in 8.5.1.

7.2.4. The currents file

According to its type – binary or formatted- the currents file name is specified in the steering file through the
keywords: BINARY CURRENTS FILE and FORMATTED CURRENTS FILE.

It is the file from which TOMAWAC reads the current field components. The current field may be either
stationary or non-stationary. The current field will be non-stationary when the keyword CONSIDERATION OF
TIDE is set to TRUE. When the current field is stationary, the keyword CONSIDERATION OF A
STATIONARY CURRENT should be set to TRUE. By default, both keywords will be set to FALSE. When
both are set to TRUE, the keywords will be inconsistent, and the program will halt.

Several commonly used formats can be read. This selection is made through the integer keyword
CURRENTS FILE FORMAT. It can be set to a value from 1 to 4

The format is 1: it is a finite-difference-typed format (as described in Appendix 8). The file is formatted and the file name should be assigned to the keyword: FORMATTED CURRENTS FILE

The format is 2: it is a point pattern-type SINUSX format (as described in Appendix 8). This file is formatted and this file name should be assigned to the keyword: FORMATTED CURRENTS FILE. This format cannot be used for reading a non-stationary current.

The format is 3: it is a TELEMAC result file of the SERAFIN standard. It is a binary file the name and its name be assigned to the keyword: BINARY CURRENTS FILE. If the current is assumed to be stationary, then the additional keyword TIME INCREMENT NUMBER IN TELEMAC FILE should be used in order to find the time step number related to the desired record. TELEMAC data other than the current components e.g. water levels, can also be read by means of this format (refer to 8.2.5).

The format is 4: data written in a different format can be read provided that the user supplies the relevant subroutine in the relevant FORTRAN file (see 8.2.3 and 8.2.6).

7.2.5. The tidal water level file

According to its type – binary or formatted- the tidal water level file name is specified in the steering file
through the keywords: BINARY TIDAL WATER LEVEL FILE or FORMATTED TIDAL WATER LEVEL FILE.

This is the file from which TOMAWAC reads the tidal water level being referred to the INITIAL STILL WATER
LEVEL. Several commonly used formats can be read. This selection is made by means of the integer
keyword TIDAL WATER LEVEL FILE FORMAT. It can be set to a value from 1 to 3.

The format is 1: it is a finite-difference-typed format (as described in APPENDIX 8). The file is formatted and the file name should be assigned to the keyword: FORMATTED TIDAL WATER LEVEL FILE.

The format is 2: it is a TELEMAC result file of the SERAFIN standard. It is a binary file and its name should be assigned to the keyword: BINARY TIDAL WATER LEVEL FILE.

The format is 3: data written in a different format can be read provided that the user supplies the relevant subroutine in the relevant FORTRAN file (see in 8.2.6).

7.2.6. The winds file

According to its type – binary or formatted- the wind file name is specified in the steering file through the
keywords: BINARY WINDS FILE or FORMATTED WINDS FILE.

This is the file from which TOMAWAC reads the information about the wind fields. As in the case of the
current, several read formats are allowed. The integer keyword WINDS FILE FORMAT can be set to values
from 1 to 4.

The format is 1: it is a WAM-cycle 4 format type(as described in APPENDIX 8). The file is formatted and the file name should be assigned to the keyword: FORMATTED WINDS FILE

The format is 2: it is a point pattern-type SINUSX format (as described in APPENDIX 8). The file is formatted and the file name should be assigned to the keyword: FORMATTED WINDS FILE.

The format is 3: it is a TELEMAC result file of the SERAFIN standard. It is a binary file and its name should be assigned to the keyword: BINARY WINDS FILE. If the wind is assumed to be stationary, then the additional keyword TIME STEP NUMBER IN TELEMAC FILE should be used in order to find the time step number related to the desired record.

The format is 4: data written in a different format can be read provided that the user supplies the relevant subroutine in the relevant FORTRAN file (see in 8.2.4).

7.2.7. The previous computation file

This previous computation file name is specified in the steering file through the character keyword:
PREVIOUS COMPUTATION FILE.

If a NEXT COMPUTATION is done, TOMAWAC fetches this file in order to initialize the directional spectrum
of wave action at every point. This file’s format, which is specific to TOMAWAC, is described in Appendix 8. It
is a binary file.

7.2.8. The global results file

The global results file name is specified in the steering file through the keyword: GLOBAL RESULTS FILE.

This file is created when a GLOBAL OUTPUT AT THE END is requested. It saves the wave action density
directional spectrum at every point in the last time step. This file format is described in APPENDIX 8.

7.2.9. The 2D results file

The 2D results file name is specified in the steering file through the character keyword: 2D RESULTS FILE.

This is the file into which TOMAWAC writes the results of the 2-dimensional variables during the
computation. It is a binary file of the SERAFIN standard. The data contained in it are in the following order:

1- all the data about the mesh geometry;

2- the names of the variables being stored;

3- for each time step, the time and the values of the variables are given for each point of the 2D mesh.

Its content varies according to the values of the following keywords:

NUMBER OF FIRST ITERATION FOR GRAPHICS PRINTOUTS: provided for determining from which time step will the data storage desirably begin, so that the file size will not be too large.

PERIOD FOR GRAPHICS PRINTOUTS: sets the period, as a number of propagation time increments, of printouts so that the file size will not be too large.

VARIABLES FOR 2D GRAPHICS PRINTOUTS: provided for specifying the list of variables to be stored into the 2D results file. Each variable is identified by 2, 3 or 4 letters (refer to Table 7.1 that lists the available variables).

For instance, if the significant wave heights, the water depths and the average wave propagation directions
are desired,

VARIABLES FOR 2D GRAPHICS PRINTOUTS = HM0,WD,DMOY
must be entered in the steering file.

7.2.10.The punctual or spectrum results file

This file’s name is specified in the steering file through the character keyword: PUNCTUAL RESULTS FILE.

This is the file into which the directional spectra of wave action at some previously specified points are stored
by TOMAWAC during the computation. These points are selected by means of the following keywords:

ABSCISSAE OF SPECTRUM PRINTOUT POINTS and ORDINATES OF SPECTRUM PRINTOUT POINTS: they are chart keywords. The maximum number of points is 19, i.e. a maximum of 19 printout points. The spectrum will be recorded at the closest point to the specified position, no spatial interpolation is made.

This file is a SERAFIN formatted file. It first includes all the data about the spectral mesh geometry, then the names-codes of displayed points. This name-code is of the type: Fa_PT2Db, where a denotes the point’s sequence order number within the list written in the steering file and b denotes the number of the closest 2D point to the specified position. Subsequently, for each graphic printout, it contains the time and the value of the directional spectrum of wave action for each pair (direction, frequency) in the spectral mesh.

The keywords PERIOD FOR GRAPHICS PRINTOUTS and NUMBER OF FIRST ITERATION FOR
GRAPHICS PRINTOUTS are shared by the two results files; thus, the printouts are synchronous for either
file.

7.2.11.The printout listing

This file contains all the messages as generated by TOMAWAC during the computation. It is the main report
of a TOMAWAC run. Its content depends on the value of the following keyword:

PERIOD FOR LISTING PRINTOUTS: this sets the time between two time steps of message transmission. This value is given in terms of the number of iterations. For example, the following sequence:

TIME INCREMENT = 30.

PERIOD FOR LISTING PRINTOUTS = 2

will result in a print in the output listing every 60 seconds of simulation.

The listing is either displayed on the monitor or saved in a file. The file name is defined by the user at the
execution of the TOMAWAC simulation (refer to APPENDIX 1).

7.2.12.The User FORTRAN file

This User FORTRAN file name is specified in the steering file through the character keyword: FORTRAN
FILE.

The FORTRAN contains all the user-modified TOMAWAC subroutines as well as the specifically developed
routines for that computation.

This file is compiled and linked during run time in order to generate the executable being used for the
simulation.

7.2.13.The auxiliary files

Other input/output files may be used by TOMAWAC.

A binary data or results file: its name is specified through the character keyword BINARY FILE 1 (Channel unit No. 24).

A formatted data or results file: its name is specified through the character keyword FORMATTED FILE 1 (Channel unit No. 26).

These files can be used either for supplying data to the program or for allowing data to be processed that are
not available in the standard results files; obviously, the user must manage the read and write operations of
these files within the FORTRAN program.

7.2.14.The dictionary file

This dictionary file contains all the information about the keywords (French/English name, default values,
type). This file can be viewed in a text editor by the user, but it must not be modified in any way.

7.2.15.The libraries

At the beginning of a computation, the main user-written FORTRAN routine is compiled, then linked in order
to generate the executable program that is subsequently run.

BIEF library: contains the computation modules related to the finite element-typed operations (operations on both matrixes and vectors). This library is shared by all the simulation models as developed by the LNHE within the TELEMAC structure (BIEF means “BIbliothèque d'Eléments Finis”, i.e. Finite Element Library).

7.3. Binary files

Binary files are an efficient way to store data on disk. However, binary files written on different computers
may differ. TOMAWAC recognizes three types of binary files, namely:

the native binary of the computer,

IBM binary (so that a file that has been generated on an IBM computer can be read), and

IEEE binary, so that these files can be read on a workstation (provided that the suitable subroutines are set up when installing TOMAWAC on the computer).

The following keywords can be used:

GEOMETRY FILE BINARY, for the geometry file,

2D RESULTS FILE BINARY for the 2D results file.

PUNCTUAL RESULTS FILE BINARY for the punctual results file.

GLOBAL RESULTS FILE BINARY, for the global results file,

PREVIOUS COMPUTATION FILE BINARY, for the previous computation file,

CURRENTS FILE BINARY, for the currents and/or TELEMAC results file.

TIDAL WATER LEVEL FILE BINARY, for the tidal water level file,

WINDS FILE BINARY, for the winds file.

BINARY 1 FILE BINARY for binary file.

In all the cases, the default value as specified in the dictionary file is 'STD' (default value of the machine
being used). The other possible values are 'IBM' and 'I3E'.

7.4. Files standard

Almost all files that were in Serafin format in previous versions of TOMAWAC, have been given a key-word
for the file format.

If the name of the file is: “GEOMETRY FILE” (“FICHIER DE GEOMETRIE”), the new keyword will be:
“GEOMETRY FILE FORMAT” (“FORMAT DU FICHIER DE GEOMETRIE”).

This format is given in 8 characters. Three choices are possible so far:

1. ‘SERAFIN ‘(do not forget the space at the end): it is the default standard within the TELEMAC
processing chain. The format is recognized by the FUDAA PRE-PRO graphics post-processor. The
RUBENS graphics post-processor reads the SERAFIN format as well, but it won't be developed anymore
and it is bound to disappear. The SERAFIN file format is described in detail in APPENDIX 8.

2. ‘SERAFIND’: Serafin format, but with double precision. Can be used for a more accurate “computation
continued” or for more accurate validations. Neither FUDAA PRE-PRO nor Rubens can read this format.

3. ‘MED ‘: this is an EDF-CEA format used in the Salomé platform, that enables to use the post-processors
of this platform. It is based on hdf5. This new format is not activated if you use the default subroutine
med.f provided, which is mostly void. If you take instead the file med.edf and rename it med.f, med
formats will be available, but two additional libraries are necessary to use this format and have to be
specified in the systel.ini file. Full instructions will be given in further releases, this is so far for internal
use at EDF.

A new file structure has been added to library BIEF for simplifying the opening/closing and reading/writing
operations with these file formats, as well as for simplifying the coupling between programmes,. The
description of this file structure and of the operations on those files are given in APPENDIX 9.

As specified in section7.2.7, a fourth binary format exists, which is specific to TOMAWAC and is used only
for saving the results when they are used to initialize a next computation. This binary file format cannot be
read by the RUBENS post-processor, or by FUDAA PRE-PRO graphics post-processor.

7.5. Bathymetry data

The bathymetry information can be supplied to TOMAWAC at two levels:

Directly in the geometry/mesh file by a bathymetry value being assigned to each node in the mesh. In this case the bathymetry data have been processed previously, running the STBTEL module or mesh generator. For example, STBTEL reads the information from one or more bottom topography files (up to 5 files) and performs an interpolation at every point within the domain;

In the form of an irregular pattern of spot heights without any necessary relation to the mesh nodes, during the TOMAWAC computation. The interpolation is then performed directly by TOMAWAC with the same algorithm as used by STBTEL. The bathymetry file name is given by the character keyword BOTTOM TOPOGRAPHY FILE. Unlike STBTEL, TOMAWAC only handles one bottom topography file. The file can be in SINUSX format or can consist of three columns X,Y,Z.

TOMAWAC also provides an opportunity to carry out a smoothing of the bathymetry in order to get a more
consistent geometry. The smoothing algorithm can be iterated several times in order to achieve more or less
extensive smoothing. The number of iterations is set using the keyword BOTTOM SMOOTHINGS and is
carried out within the CORFON subroutine. This keyword’s default value is 0. (also refer to the programming of
the CORFON user subroutine in 8.6.1).

NOTE: the bathymetry data should preferably be supplied to TOMAWAC in the form of water depth and not
of water height. If necessary, a conversion can be performed in the CORFON subroutine.

8. Controlling the simulation

8.1. General parameterisation

The general parameterisation of the computation is controlled using the steering file.

8.1.1. Computation title

The computation case title is specified by the keyword TITLE

8.1.2. Computation time

The time data are prescribed using the two keywords TIME STEP (real) and NUMBER OF TIME
INCREMENTS (integer). The former sets the elapsed time, in seconds, between two consecutive
computation instants (but not necessarily two outputs in the results file). The number of time steps is for
setting the overall computation time (which is obviously equal to the time increment value multiplied by the
number of time steps).

An additional keyword also refers to time. This is the DATE OF COMPUTATION BEGINNING, which is used
to identify the time in relation to the date/times written in the WINDS FILE (refer to APPENDIX 8, WAM-typed
format). The convention adopted for writing the DATE OF COMPUTATION BEGINNING is yymmddhhmm.
For instance, 0311051110 corresponds to November 5th, 2003, at 11.10 AM.

Note that when a computation is resumed, the initial time of the new computation corresponds to the last
time increment in the previous computation (i.e. the computation is not resumed at time zero).

8.1.3. Spectral discretisation

The spectral discretisation is defined by the following 5 keywords:

NUMBER OF DIRECTIONS,

NUMBER OF FREQUENCIES,

MINIMUM FREQUENCY,

FREQUENTIAL RATIO,

SPECTRUM TAIL FACTOR.

It should be reminded that the directions are evenly distributed from 0 to 360 degrees. Two conventions can
be chosen by the user for expressing the wave propagation by means of the keyword TRIGONOMETRICAL
CONVENTION (logical). The trigonometrical convention locates the wave propagation from the positive X
axis and the direction of rotation is in a counterclockwise direction. The default convention is the nautical
convention (TRIGONOMETRICAL CONVENTION = NO) that locates the propagation direction in relation to
the true “North”, i.e. the Y axis. The selected direction of rotation is a clockwise direction. The direction will
then correspond to the “heading” in the sense of the navigational maps, i.e. the direction the waves are
propagating towards.

The frequencies are distributed geometrically in accordance with the following relation:

(k = 1,NF)

where is the MINIMUM FREQUENCY, r is the FREQUENTIAL RATIO and NF is the NUMBER OF FREQUENCIES

In order to take the contribution of high frequencies (higher than the maximum discretised frequency) into
account in the computations, it is assumed that the decay of the spectrum follows a law that is of the type f-n.
The keyword SPECTRUM TAIL FACTOR corresponds to the value of n.

8.1.4. Release

When generating the executable, the release number of libraries being used for editing the links is indirectly
provided by the keyword TOMAWAC RELEASE NUMBER. By default, TOMAWAC release 6.0 utilises the
6.0 releases of the TELEMAC system libraries.

8.1.5. Environment

When a vector computer is used, the CPU vector length used in the forced vectorisation technique can be
specified by means of the keyword VECTOR LENGTH. The default value is 1 and is appropriate for scalar
machines such as the present workstations. If a value of 1 is used on a vector machine, then the advantage
of the vectorisation loops (although they are few in TOMAWAC) is lost.

8.2. Computation options

8.2.1. Co-ordinate system

Cartesian co-ordinates (expressed in meters) are used by default. For domains of a large extent, working in
spherical co-ordinates may become necessary. The value of the logical keyword SPHERICAL
COORDINATES should then be set to “TRUE”. The co-ordinates are then expressed in degrees.

8.2.2. Finite depth

In nearshore areas, wave conditions will be influenced by the water depth and therefore the bottom effect
can no longer be ignored. This is the default case in TOMAWAC: the keyword INFINITE DEPTH is set to
“FALSE”. In situations where depths effects are to be explicitly ignored the keyword INFINITE DEPTH should
be set to “TRUE”.

8.2.3. Taking a stationary current into account

A stationary current can be taken into account in the TOMAWAC release 6.0. The relevant logical keyword is
CONSIDERATION OF A STATIONARY CURRENT. The current affects mainly the convection step.

The current can be specified in various ways:

When the current is either constant over the domain or can be described analytically, the ANACOS.f subroutine can be included in the FORTRAN file and modified accordingly. In this subroutine the UC and VC
are NPOIN2-sized (number of points in the horizontal mesh) vectors and correspond to the components
along the X and Y axes of the current, respectively. This is how the current is specified when the keyword
FORMATTED CURRENTS FILE or BINARY CURRENTS FILE is not specified, whereas the keyword
CONSIDERATION OF A STATIONARY CURRENT is “TRUE”.

TOMAWAC can also take into account a current provided in a binary or formatted file. The keyword BINARY
CURRENTS FILE or FORMATTED CURRENTS FILE should then be given a value (the name of the file).
Three different formats are available to read this data. The value corresponds to the keyword CURRENTS
FILE FORMAT (see 7.2.4). When the currents file is taken from TELEMAC-2D, then an additional keyword
should be specified, namely the TIME INCREMENT NUMBER IN TELEMAC FILE. This locates the desired
record.

If the predefined formats cannot be used, the COUUTI.f subroutine can be included in the FORTRAN and
modified accordingly, specifying the format 4 for the INDIC FORTRAN variable in the CAS file. The current
data are read from the file and are interpolated onto the nodes of the computation mesh.

8.2.4. Taking a wind into account

Consideration of a wind is specified by the logical keyword CONSIDERATION OF A WIND. The wind may be
either stationary or variable in time and is specified by means of the logical keyword STATIONARY WIND.
When the wind can be described analytically, the user subroutine ANAVEN.f can be used. The wind is fully
specified when the keywords FORMATTED WINDS FILE and BINARY WINDS FILE do not have any value,
whereas the keyword CONSIDERATION OF A WIND is “TRUE”.

TOMAWAC can also take into account a wind given in a binary or formatted file. In this case a value (the
name of the file) should be assigned to the keyword FORMATTED WINDS FILE or BINARY WINDS FILE.
The available formats for reading out these data correspond to the keyword WINDS FILE FORMAT (see in
section 7.2.6. and APPENDIX 8).

When these predefined formats cannot be used, the subroutine VENUTI.f can be included in the FORTRAN
file and modified accordingly, specifying the format 4 for the INDIV FORTRAN variable. On completion of
reading the winds file, the wind components are used as such if provided on the computational mesh, or
interpolated over that mesh in provided on a different grid.

NOTE: an interpolation between two different meshes of equivalent sizes is usually computationally very
expensive. Although possible, it is highly inadvisable, particularly as regards to the wind, since this is a timevarying
data item. In such cases a pre-interpolation over the computation mesh, e.g using STBTEL is
recommended, followed by the reading of the wind data in format 3. Alternatively this pre-interpolation can be
performed by means of the FASP subroutine from the utile library.

8.2.5. Recovering a TELEMAC data item

Recovering a 2D result data item from a TELEMAC-2D hydrodynamic computation might be of interest, e.g.
the value of wind-driven surge at every point. To avoid an increase in the number of files the BINARY
CURRENTS FILE is used to specify this input file. The keyword CURRENTS FILE FORMAT should then be
set to 3. This option is further specified using the logical keyword RECOVERY OF TELEMAC DATA ITEM.
The data item is located within the file through the TIME INCREMENT NUMBER IN TELEMAC FILE and the
RANK OF THE TELEMAC DATA ITEM TO BE RECOVERED that corresponds to desired variable's
sequence number in the record. is the needed data are read from the file and interpolated over the
computation mesh.

NOTE: a TELEMAC data item and the components of a current can both be read simultaneously provided
that they occur in the same file at the same record.

The recovered variable, which is interpolated over the mesh, can be utilised in the subroutine VARTEL.

8.2.6. Taking the tide into account

Tide-induced effects, i.e. unsteady/non-stationary water levels and currents can be taken in to account. The relevant logic keyword is CONSIDERATION OF TIDE.

In order to take tide into account, a current and a tide water depth that is referenced in relation to the
“INITIAL STILL WATER LEVEL” must be specified. These data can be initialized in various ways:

Should the tide be easy to describe analytically, the ANAMAR.f subroutine, can be included in the FORTRAN
file and modified accordingly. In the subroutine the terms UC and VC, ZM and DZHDT are NPOIN2-sized
(number of points in the horizontal mesh) vectors and correspond to the current components along the X and
Y axes, the tidal water level in relation to the “INITIAL STILL WATER LEVEL” and the water depth variation
in time, respectively. An analytical expression must be assigned to all of these vectors.

TOMAWAC can also take into account a current that is given in a binary or formatted file. A value (the name
of the currents file)should be assigned to the keyword BINARY CURRENTS FILE or FORMATTED
CURRENTS FILE. Two different formats are available for reading these data. This format is specified using
the keyword CURRENTS FILE FORMAT (see 7.2.4). When these predefined formats cannot be used, the
user subroutine COUUTI.f can be included in the FORTRAN file and modified accordingly. In such cases
the CURRENTS FILE FORMAT (FORTRAN variable INDIC) must be set to 4 in the CAS file. Once the data
of the currents file are read, the current components are interpolated over the computation mesh.

The tidal water level can be also provided in a binary or formatted file. A value (the name of the water level
file) should be assigned to the keyword BINARY TIDAL WATER LEVEL FILE or FORMATTED TIDAL
WATER LEVEL FILE. Two different predefined formats are available for reading this data. The format type is
specified using the keyword TIDAL WATER LEVEL FILE FORMAT (refer to 7.2.5). If the user chooses the
Serafin format (i.e. TELEMAC-2D format), then the RANK OF THE WATER LEVEL DATA IN THE
TELEMAC FILE must also be specified. When the predefined formats cannot be used, the user subroutine
MARUTI.f file can be included in the FORTRAN file and modified accordingly. In such cases the TIDAL
WATER LEVEL FILE FORMAT (the INDIM FORTRAN variable) must be set to 4 in the CAS file.

Both currents and tidal water levels will be updated upon each TIDE REFRESHING PERIOD. This keyword
corresponds to an integer multiple of the propagation TIME

It is possible to directly couple a TOMAWAC and a TELEMAC simulation (either TELEMAC-2D or
TELEMAC-3D) to represent the wave-current interactions more precisely.

In TOMAWAC, when the keywords CONSIDERATION OF A TIDE or CONSIDERATION OF A
STATIONARY CURRENT are used, the current is imposed as an input data: in this case only the effect of
the current on the waves is taken into account, but not the effect of the waves on the current. The current
imposed, therefore, is not affected by the waves and can differ from the real current interacting with the
waves.

By using a direct coupling TELEMAC-TOMAWAC it is possible to represent wave current interactions in both
directions: TELEMAC transfers to TOMAWAC the updated values of current velocities and water depths,
while TOMAWAC solves the wave action density conservation equation with reference to those current and
water depth values and returns to TELEMAC the updated values of the wave driving forces FX and FY acting
on the current.

To directly couple a TELEMAC model with a TOMAWAC model, the following conditions must be satisfied:

the TELEMAC and TOMAWAC models must have the same horizontal mesh

the TELEMAC and TOMAWAC simulations must have the same simulation time length (given by the time step multiplied by the number of time steps)

the time steps set for the two simulations must be equal or multiple of each other

current and/or water level file cannot be used as input data for the TOMAWAC simulation

the driving force along X and Y (FX, FY) must be set among the 2D output variables in the steering file

In case of direct coupling TELEMAC-TOMAWAC, TELEMAC is the main programme and calls the
TOMAWAC subroutine WAC.f, which is the main subroutine of TOMAWAC and solves the equation of
generation and propagation of the directional wave spectrum.

The TELEMAC-TOMAWAC coupling is set via four keywords in the TELEMAC steering file:

COUPLING WITH, to which the value 'TOMAWAC' must be assigned

WAVE DRIVEN CURRENT must be set to 1

TOMAWAC STEERING FILE, which specifies the name of the TOMAWAC steering file: its path must be specified with reference to the working directory of the TELEMAC steering file.

COUPLING PERIOD FOR TOMAWAC (variable PERCOU_WAC), which specifies every how many TELEMAC time steps TOMAWAC is called:

If PERCOU_WAC = 1 (default value), then the TELEMAC simulation time step is equal to or

larger than the TOMAWAC time step

If PERCOU_WAC > 1, then the TOMAWAC simulation time step is larger than the

TELEMAC time step

All the file names specified in the TOMAWAC steering file (geometry file, boundary conditions file, Fortran
file, wind file, …) must be given as paths relative to the working directory of the TELEMAC steering file.

For more information concerning the modelling options in TELEMAC, please refer to the TELEMAC
documentation.

8.2.8. Convection step

For specific validation tests, for example, it may be interesting to drop the convection step and only consider
the effect of the source terms. To do this requires assigning the “FALSE” value to the keyword
CONSIDERATION OF PROPAGATION.

8.3. Parameterising the source term integration step

8.3.1. Introduction

When it is required to take into account the source/sink terms, the logic keyword CONSIDERATION OF
SOURCE TERMS should be set to “TRUE”.

It has been shown that the source term integration may require a shorter time step than the time step that is
used for convection. The TIME STEP in the steering file corresponds to the convection time step. The source
term integration step is controlled using the integer keyword NUMBER OF ITERATIONS FOR THE SOURCE
TERMS. This keyword is set to the number of source terms integration time steps that will be conducted after
each convection step (default value = 1). The effective time-step used for source term integration is thus:
(TIME STEP)/(NUMBER OF ITERATIONS FOR THE SOURCE TERMS).

Depending on the source/sink terms, two different schemes are used for time integration:

For the source/sink terms which are dominant in large and medium water depths (namely wind input, white-capping dissipation, nonlinear quadruplet interactions and bottom friction) a scheme with variable implicitation level is used (see section 8.3.2).

For the source/sink terms which are dominant in shallow water depths (namely depth-induced breaking and, nonlinear triad interactions) an explicit scheme is used, possibly with sub-steps to cover one source term time-step (see section 8.3.3).

8.3.2. Source/sink terms in large and medium water depth

8.3.2.1. Wind input

In TOMAWAC three wind generation models and a linear wave growth model are available (see section
4.2.3.2 for details).

If the keyword WIND GENERATION is set to 0, no wind input will be taken into account. If, on the other
hand, a strictly positive value (1, 2, …) is chosen, the corresponding wind input formulation will be taken into
account.

The Janssen formulation is activated using the value 1 for the WIND GENERATION keyword. Janssen’s
formulation requires several additional data, specified by the following keywords:

AIR DENSITY,

WATER DENSITY,

WIND GENERATION COEFFICIENT,

VON KARMAN CONSTANT,

CHARNOCK CONSTANT,

SHIFT GROWING CURVE DUE TO WIND,

WIND DRAG COEFFICIENT,

WIND MEASUREMENT LEVEL.

As a general rule, the default values for these keywords shall not be modified.

The Snyder formulation is activated setting to 2 the value for the WIND GENERATION keyword and uses
two parameters, specified by the following keywords:

AIR DENSITY,

WATER DENSITY,

The Yan formulation is activated setting to 3 the value for the WIND GENERATION keyword and uses four
parameters, specified by the following keywords:

YAN GENERATION COEFFICIENT D

YAN GENERATION COEFFICIENT E

YAN GENERATION COEFFICIENT F

YAN GENERATION COEFFICIENT H

As a general rule, the default values for these keywords shall not be modified.
The linear wave growth model, based on the Cavaleri & Malanotte-Rizzoli formulation, is activated setting to
1 the keyword LINEAR WAVE GROWTH. This formulation does not require any additional data to be
specified by the user.

8.3.2.2. White capping dissipation

In TOMAWAC two whitecapping dissipation models are available (see section 4.2.3.2 for details).

If the integer keyword WHITE CAPPING DISSIPATION is set to 0, this source term will be ignored. If a
strictly positive value (1, 2, …) is selected, the corresponding formulation will be taken into account.

The Komen's formulation corresponds to the value 1 of the WHITE CAPPING DISSIPATION keyword. This
formulation requires two complementary data, specified by the following keywords:

WHITE CAPPING DISSIPATION COEFFICIENT

WHITE CAPPING WEIGHTING COEFFICIENT.

As a general rule, the default values for these keywords shall not be modified.

The Westhuysen formulation corresponds to the value 2 of the WHITE CAPPING DISSIPATION keyword.
This formulation requires four complementary data, specified by the following keywords:

WESTHUYSEN DISSIPATION COEFFICIENT'

SATURATION THRESHOLD FOR THE DISSIPATION

WESTHUYSEN WHITE CAPPING DISSIPATION

WESTHUYSEN WEIGHTING COEFFICIENT'

As a general rule, the default values for these keywords shall not be modified.

8.3.2.3. Bottom friction dissipation

If the integer keyword BOTTOM FRICTION DISSIPATION is set to 0, this source term will be ignored. If a
strictly positive value (1, 2, …) is chosen, the corresponding formulation will be taken into account. This
source term is only taken account if the keyword INFINITE DEPTH is set to “FALSE”.

The only formulation implemented is from Hasselmann (see section 4.2.3.4 for details) This formulation is
specified by setting the BOTTOM FRICTION DISSIPATION keyword to 1. This formulation requires the
specification of the keyword:

BOTTOM FRICTION COEFFICIENT.

As a general rule, the default value for this keyword shall not be modified.

8.3.2.4. Non-linear transfers between quadruplets

Three non-linear quadruplet interactions models are available in TOMAWAC (see in section 4.2.3.6). The
non-linear quadruplet interactions are activated through the keyword NON-LINEAR TRANSFERS BETWEEN
FREQUENCIES in the steering file; the keyword can take four values.

If the integer keyword NON-LINEAR TRANSFERS BETWEEN FREQUENCIES is set to 0, this source term
will not be taken into account. If a strictly positive value (1, 2, …) is chosen, the corresponding formulation will
be taken into account.

The DIA method (formulation from Hasselmann et al., 1985) corresponds to the value 1 for the NON-LINEAR
TRANSFERS BETWEEN FREQUENCIES keyword. This formulation does not require any additional data to
be specified by the user.

The MDIA method (formulation from Tolman, 2004) corresponds to the value 2 for the NON-LINEAR
TRANSFERS BETWEEN FREQUENCIES keyword. This formulation does not require any additional data to
be specified by the user.

The qausi-excat GQM method (formulation based on Lavrenov, 2001) corresponds to the value 3 for the
NON-LINEAR TRANSFERS BETWEEN FREQUENCIES keyword. This formulation requires the
specification of six keywords:

SETTING FOR INTEGRATION ON OMEGA1

SETTING FOR INTEGRATION ON THETA1

SETTING FOR INTEGRATION ON OMEGA2

THRESHOLD0 FOR CONFIGURATIONS ELIMINATION

THRESHOLD1 FOR CONFIGURATIONS ELIMINATION

THRESHOLD2 FOR CONFIGURATIONS ELIMINATION

8.3.3. Source/sink terms in shallow water depth

8.3.3.1. Time integration scheme and time step

As mentioned above, the depth-induced breaking and nonlinear triad interaction terms are time-integrated
with an explicit scheme.

As found practically, contributions from these source terms can be overestimated if the selected time step for steps which are specific to these source terms through the keyword NUMBER OF BREAKING TIME STEPS.
These time sub-steps are arranged in a geometrical progression, i.e. they are defined in the following way:

where the geometrical ratio q is specified through the keyword: COEFFICIENT OF THE TIME SUBINCREMENTS
FOR BREAKING
source term integration is too long. In order to avoid this, TOMAWAC can perform a number of time sub

In order to limit this number of time-steps, TOMAWAC first clips the wave height by setting a maximum
ratio (D is the water depth) to 1. This ratio can be modified by means of the keyword MAXIMUM
VALUE OF THE RATIO HM0 TO D. However, this is generally not advisable.

8.3.3.2. Wave breaking dissipation

If the integer keyword DEPTH-INDUCED BREAKING DISSIPATION is taken as 0, this source term will be
ignored. If a strictly positive value (1, 2, …) is chosen, the corresponding formulation will be taken into
account.

Four formulations have been implemented (see section 4.2.3.5 for details):

1: Battjes and Janssen's model (1978)

This formulation requires additional data to be provided, specified by the following keywords:

DEPTH-INDUCED BREAKING 1 (BJ) COEFFICIENT ALPHA

DEPTH-INDUCED BREAKING 1 (BJ) COEFFICIENT GAMMA1

DEPTH-INDUCED BREAKING 1 (BJ) COEFFICIENT GAMMA2

DEPTH-INDUCED BREAKING 1 (BJ)CHARACTERISTIC FREQUENCY

DEPTH-INDUCED BREAKING 1 (BJ)QB COMPUTATION METHOD

2: Thornton and Guza's model (1983)

This formulation requires additional data to be provided, specified by the following keywords:

DEPTH-INDUCED BREAKING 2 (TG) COEFFICIENT B

DEPTH-INDUCED BREAKING 2 (TG) COEFFICIENT GAMMA

DEPTH-INDUCED BREAKING 2 (TG) WEIGHTING FUNCTION

DEPTH-INDUCED BREAKING 2 (TG) CHARACTERISTIC FREQUENCY

3: Roelvink's model (1993)

This formulation requires additional data to be provided, specified by the following keywords:

DEPTH-INDUCED BREAKING 3 (RO) COEFFICIENT ALPHA

DEPTH-INDUCED BREAKING 3 (RO) COEFFICIENT GAMMA

DEPTH-INDUCED BREAKING 3 (RO) COEFFICIENT GAMMA2

DEPTH-INDUCED BREAKING 3 (RO) WAVE HEIGHT DISTRIBUTION

DEPTH-INDUCED BREAKING 3 (RO)EXPONENT WEIGHTING FUNCTION

DEPTH-INDUCED BREAKING 3 (RO)CHARACTERISTIC FREQUENCY

4: Izumiya and Horikawa's model (1984)

This formulation requires additional data to be provided, specified by the following keywords:

DEPTH-INDUCED BREAKING 4 (IH) COEFFICIENT BETA0

DEPTH-INDUCED BREAKING 4 (IH) COEFFICIENT M2STAR

DEPTH-INDUCED BREAKING 4 (IH) CHARACTERISTIC FREQUENCY

8.3.3.3. Triad interactions

If the integer keyword TRIAD INTERACTIONS is set to 0, this source term will be ignored. If a strictly positive
value (1, 2, …) is specified, the corresponding formulation will be taken into account.

Two formulations (see section 4.2.3.7) have been implemented.

1: LTA model

This formulation requires additional associated data to be specified using the following keywords:

TRIADS 1 (LTA) COEFFICIENT ALPHA

TRIADS 1 (LTA) COEFFICIENT RFMLTA

2: SPB model

This formulation requires additional associated data to be specified using the following keywords:

TRIADS 2 (SPB) COEFFICIENT K

TRIADS 2 (SPB) LOWER DIRECTIONAL BOUNDARY

TRIADS 2 (SPB) UPPER DIRECTIONAL BOUNDARY

ATTENTION: the SPB model is very time-consuming; compared to the LTA model formulation, it requires a
computational time approximately 700 times higher.

8.4. Prescribing the initial conditions

Initial conditions can be prescribed using the integer keyword TYPE OF INITIAL DIRECTIONAL
SPECTRUM;

Table 8.1 shows all the available options in TOMAWAC for computing the frequential and directional
distribution of wave action. It should be remembered that the variance density directional spectrum is
computed as the product:

where E(f) here denotes the variance density spectrum and denotes the angular distribution function.
It is reminded that the parameterised JONSWAP spectrum is defined as:

where :

and that the TMA spectrum is a depth-corrected JONSWAP-type spectrum.
A parameterised spectrum with two directional peaks can also be generated. In such case, both main
propagation directions ( and ) as well as the weighting factor ( ) between the two power peaks:

must be specified.

and , in the equation above, are automatically computed by TOMAWAC in order to normalize the
angular distribution function.
Three angular distribution functions may be chosen using the keyword INITIAL ANGULAR FUNCTION
DISTRIBUTION, which correspond to the following options:

1. model in
2. model in
3. model in

with being the main sea-state propagation direction and s being the directional spread.

The forementioned constants can be specified using the following of keywords:

The keyword SPECTRUM ENERGY THRESHOLD is used whatever option is chosen. It is only useful for
comparisons with the WAM model.

Specific initial conditions can be prescribed directly for the directional spectrum of wave action using the
condiw.f subroutine which can be called from the speini.f subroutine.

8.5. Prescribing the boundary conditions

The boundary conditions are prescribed over the relative spectrum of wave action, i.e. expressed in coordinate
system that moves with the current.

Only two kinds of boundary conditions are available in TOMAWAC.

The first one corresponds to a free boundary i.e. a boundary that fully absorbs the wave energy. It may be a
liquid boundary: it is then assumed that the waves propagate beyond the domain and nothing else enters it.
It may be a solid boundary: it is then assumed that the shore fully absorbs the wave energy.

The second one corresponds to a boundary with a prescribed value. In this case the wave action spectrum is
then strictly imposed at each point along that boundary. This boundary condition allows wave energy to enter
the computational domain.

The boundary conditions are specified using the boundary conditions (CONLIM) file, the steering (CAS) file
and the LIMWAC.f.

8.5.1. The boundary conditions file

The boundary condition (CONLIM) file is normally supplied by STBTEL (or other TELEMAC mesh
generators), but it can also be generated by means of a text editor. Each line in this file is assigned to one
point of the boundary and listed in sequential order in terms of the boundary node numbers. The numbering
of the boundary points first delineates the domain contour in a counterclockwise direction, then the islands in
the clockwise direction. The total number of edge points is noted as NPTFR.

13 values are given for each point. Only data in colums 1, 12 and 13 are used by TOMAWAC:

The 13th data column (integer variable IPTFR) corresponds to the boundary point number ranked in terms of the boundary point numbering.

The 12th data column (integer variable IPOIN) corresponds to the global number of the point in the 2D mesh.

Lastly, the 1st data column (integer variable LIFBOR) corresponds to the kind of boundary condition. Consistent with TELEMAC-2D, its value is 2 in the case of a free boundary and 5 in the case of a boundary with a prescribed value.

8.5.2. Prescribing the boundary conditions in the CAS file

Boundary conditions prescribed using the CAS file will necessarily be homogeneous all along the domain
entry boundaries.

The boundary conditions can be prescribed by means of the integer keyword TYPE OF BOUNDARY
DIRECTIONAL SPECTRUM

Table 8.1 (see section 8.4) presents all the spectrum types available in TOMAWAC. The constants given in
Table 8.1 can be prescribed using the following keywords:

: BOUNDARY SIGNIFICANT WAVE HEIGHT,

: BOUNDARY PEAK FREQUENCY

: BOUNDARY PEAK FACTOR

: BOUNDARY SPECTRUM VALUE OF SIGMA-A

: BOUNDARY SPECTRUM VALUE OF SIGMA-B

: BOUNDARY PHILLIPS CONSTANT

: BOUNDARY MEAN FETCH VALUE

: BOUNDARY MAXIMUM PEAK FREQUENCY

: BOUNDARY MAIN DIRECTION 1

: BOUNDARY DIRECTIONAL SPREAD 1

: BOUNDARY MAIN DIRECTION 2

: BOUNDARY DIRECTIONAL SPREAD 2

: BOUNDARY WEIGHTING FACTOR FOR ADF

Three angular distribution functions have been implemented and can be selected using of the keyword:
BOUNDARY ANGULAR DISTRIBUTION FUNCTION, which corresponds to the following options:

1: model in

2: model in

3: model in

Since the boundary spectrum computation procedures are similar to those for the initial spectrum, refer to
section 8.4 for further details.

8.5.3. The LIMWAC user subroutine

It should be reminded that the spectrum is discretised over both frequencies and directions and that it is a
relative spectrum, i.e. expressed in a coordinate system that moves with the current.

The subroutine LIMWAC, in its original version, allows to impose the spectrum components at each point of
a boundary with a prescribed value. The spectrum components are calculated from the parameters specified
in the CAS file (see section 7.2.1). This subroutine, however, can easily be modified to specify e.g. nonhomogeneous
(in space) boundary conditions. When such specific boundary conditions are required, these
will ideally be incorporated in the user part provided in the code of the LIMWAC procedure. The keyword
BOUNDARY SPECTRUM MODIFIED BY USER must also be set to YES.

8.6. Some useful subroutines

8.6.1. Modification of bottom topography: CORFON subroutine

The seabed levels can be modified in two different ways, as already stated in section 7.5.

The seabed levels can be modifed at the beginning of the computation using the CORFON subroutine, which
is called once at the beginning of the computation. This subroutine allows the value of the ZF variable to be
modified at each mesh point. For this purpose, a number of variables such as, for instance, the point
coordinates, the element area values, the connectivity table, etc., are provided.

By default, the CORFON subroutine performs the same number of bottom smoothing iterations as LISFON,
i.e. the same value as specified by the integer keyword BOTTOM SMOOTHINGS.

Note that the CORFON subroutine is not called in case the computation is initialized with the result of a
former TOMAWAC run (“hot start”).

This subroutine is part of the TELEMAC-2D library and is listed in APPENDIX 2.

8.6.2. Modifying the co-ordinates: CORRXY subroutine

TOMAWAC allows the mesh point co-ordinates to be modified at the beginning of the computation, so that
an up-scaling (switching from a small scale model to a full-size model), a rotation or a translation can be
performed.

Such changes are made using the CORRXY subroutine from the BIEF library, which is called in at the
beginning of the computation. This subroutine is void by default and provides, in the form of a comment, an example of programming relevant to a change of scale and origin. It is part of the TELEMAC-2D library and
is listed in APPENDIX 2.

8.6.3. Operations on vectors: OV subroutine

The BIEF library has a range of very useful subroutines including, in particular, subroutines for operations on
vectors. A number of relations have been programmed so that loops can be replaced by a mere procedure
call.

The syntax is as follows:

CALL OV(OP, X, Y, Z, C, NPOIN)

Where OP is a string of exactly 8 digits that is indicative of the operation about to be performed on the X, Y,
Z vectors and the constant C. The result is the vector X.

Example:

CALL OV('X=X+Y ', X, Y, Z, C, NPOIN)

Y is added to X, the result will be stored in X.

A full list of available operations are given in the Guide to programming in the Telemac system version 6.0.