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.

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.

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…

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”).

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).

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.

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 .

It is a function (in Joule.Hz-1.rad-1) that depends on:

• wave frequency f (in Hertz), conventionally only positive (ranging from 0 to )
• propagation direction , ranging within a length interval.

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.

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:

The n-order (n = 0, 1, 2,…) moments of the variance density directional spectrum are defined as:

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]).

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 quantity, noted as , will remain constant (without considering the source and sink terms) even though the propagation medium is neither homogeneous nor steady [Komen et al., 1994] [Willebrand, 1975] [Phillips, 1977] [Bretherton, 1969].

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

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:

= (absolute frequency; propagation direction)

= (relative frequency; propagation direction)

= (wave number; propagation direction)

= (wave number vector)

For the numerical resolution of equations, the model TOMAWAC uses the pair = (relative frequency; propagation direction)

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

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.

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:

• bathymetry and sea bottom geometry (bottom friction, refraction, surf-breaking, non-linear effects of interactions with the bottom, sand rippling…)
• atmospheric circulation (wind and pressure effects)
• tide pattern (variation of currents and water heights),
• three-dimensional oceanic circulation currents,
• over/underelevations caused by exceptional weather events, resulting in sea levels variations up to several meters (storm, surges).

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.

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.

• bottom friction-induced dissipation, mainly occurring in shallow water (bottom grain size distribution, ripples, percolation…)

• dissipation through bathymetric breaking. As the waves come near the coast, they swell due to shoaling until they break when they become to steep.

• dissipation through wave blocking due to strong opposing currents.

—> Non-linear energy transfer conservative processes:

• non-linear resonant quadruplet interactions, which is the exchange process prevailing at great depths.

• non-linear triad interactions, which become the prevailing process at small depths.

-—> Wave propagation-related processes:

• wave propagation due to the wave group velocity and, in case, to the velocity of the medium in which it propagates (sea currents).

• depth-induced refraction which, at small depths, modifies the directions of the wave-ray and then implies an energy transfer over the propagation directions.

• shoaling: wave height variation process as the water depth decreases, due to the reduced wavelength and variation of energy propagation velocity.

• current-induced refraction which also causes a deviation of the wave-ray and an energy transfer over the propagation directions.

• interactions with unsteady currents, inducing frequency transfers (e.g. as regards tidal seas).

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.

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):

The transfer rates are given by the linear wave theory [Chaloin, 1989] [Komen et al., 1994] [Mei, 1983] [Tolman, 1991]:

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).
• Discretization of angular spectrum uses the pair (,) = (relative frequency ; propagation direction).

The following conventions are adopted for writing the equations :

• 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.

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.

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.

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:

wherein :

: wind-driven wave generation : whitecapping-induced energy dissipation : non-linear quadruplet interactions : bottom friction-induced energy dissipation : bathymetric breaking-induced energy dissipation : non-linear triad interactions

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.

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:

• 0 : no wind input (default value)
• 1 : Janssen’s model [Janssen, 1989] [Janssen, 1991] (WAM cycle 4) (see in paragraph 4.2.3.2.1).
• 2 : Snyder et al. model [Snyder et al., 1981] (see in paragraph 4.2.3.2.2).
• 3 : Yan's model [Yan, 1987] (see in paragraph 4.2.3.2.3)

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).

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.

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.

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 .

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.

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:

0 no whitecapping-induced dissipation (default value)

1 Komen et al. [Komen et al., 1984] and Janssen’s [Janssen, 1991] dissipation model.

2 Westhuysen et al. [Westhuysen et al., 2007] dissipation model

For a more detailed description of the issues related to the whitecapping dissipation modelling and of the recent advances in this field, reference can be made to [WISE group, 2007].

In deep water, that term is written as follows in TOMAWAC:

With a finite water height, TOMAWAC uses the following formulation:

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.

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.

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.

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.

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)

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:

• « DEPTH-INDUCED BREAKING 2 (TG) WEIGHTING FUNCTION Two values are possible:

  1. weighting function 1 (see in (4.42.a))

  2. weighting function 2 (see in (4.42.b)) (default value)

• « DEPTH-INDUCED BREAKING 2 (TG) CHARACTERISTIC FREQUENCY. Six values are possible:

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

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

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

  4. discrete peak frequency:

  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.

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.
• « DEPTH-INDUCED BREAKING 3 (RO) EXPONENT WEIGHTING FUNCTION, corresponding to the N variable. Its default value in the model is 10.
• « DEPTH-INDUCED BREAKING 3 (RO) CHARACTERISTIC FREQUENCY Six values are possible:

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

  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:

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.
• « DEPTH-INDUCED BREAKING 4 (IH) CHARACTERISTIC FREQUENCY. Six values are possible:

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

  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:

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:

0 no non-linear quadruplet interaction (default value)

1 DIA method (Discrete Interaction Approximation) of Hasselmann et al. [Hasselmann et al., 1985] which is a discrete parameterization of the exact computation operator as proposed by Hasselmann [Hasselmann, 1962] [Hasselmann, 1962].

2 MDIA method (Multiple DIA) as proposed by Tolman [Tolman, 2004]

3 Quasi exact GQM method (Gaussian Quadrature Method) as introduced by Lavrenov [Lavrenov, 2001] and implemented by Gagnaire-Renou et al. [Gagnaire-Renou et al., 2010].

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.

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 .

Two constants can me modified through keywords:

• constant corresponding to the keyword TRIADS 1 (LTA) COEFFICIENT ALPHA. Its default value is
• constant corresponding to the keyword TRIADS 1 (LTA) COEFFICIENT RFMLTA. Its default value is

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.

Three constants can be modified through keywords:

• constant K corresponding to the keyword 'TRIADS 2 (SPB) COEFFICIENT K. Its default value is K = 0.34
• the lower and upper directional markers corresponding to the keywords TRIADS 2 (SPB) LOWER DIRECTIONAL BOUNDARY and TRIADS 2 (SPB) UPPER DIRECTIONAL BOUNDARY. Their respective default values are 0 and 360

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).

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.

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)

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 !

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.

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.

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.

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.

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.

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].

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.

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.

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).

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.

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.

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.

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).

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).

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).

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.

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.

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.

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.

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).

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.

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.

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.

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.

The following libraries are used during the link editing operation:

• TOMAWAC library: contains the specific subroutines of the TOMAWAC computation model.
• telemac library: contains the specific subroutines of the TELEMAC-2D computation model.
• utile library: contains a number of utility subroutines, such as the file read and write routines.
• damocles library: contains the subroutines handling the steering file reading.
• 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).

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'.

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.

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.

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

The computation case title is specified by the keyword TITLE

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).

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.

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.

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.

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.

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”.

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.

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.

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.

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.

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.

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).

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.

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.

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.

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

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.

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

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.

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:

: INITIAL SIGNIFICANT WAVE HEIGHT
: INITIAL PEAK FREQUENCY
: INITIAL PEAK FACTOR
: INITIAL VALUE OF SIGMA-A FOR SPECTRUM
: INITIAL VALUE OF SIGMA-B FOR SPECTRUM
: INITIAL PHILLIPS CONSTANT
: INITIAL MEAN FETCH VALUE
: INITIAL MAXIMUM PEAK FREQUENCY
: INITIAL MAIN DIRECTION 1
: INITIAL DIRECTIONAL SPREAD 1
: INITIAL MAIN DIRECTION 2
: INITIAL DIRECTIONAL SPREAD 2
: INITIAL WEIGHTING FACTOR FOR ADF

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.

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.

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.

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.

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.

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.

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.

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.