Table of Contents

TOMAWAC User Manual links from User Manuals, TOMAWAC


1.Introduction


1.1 Generals

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

TELEMAC-based Operational Model Addressing Wave Action Computation

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

1.2 Implementing TOMAWAC

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

1.3 TOMAWAC general purposes

TOMAWAC can be used for three types of applications:

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:


2.Representing waves in TOMAWAC


2.1 General definition of waves

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

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

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

2.2 Plane monochromatic waves

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

\begin{equation*} \eta(x,y,t) = a.cos \bigg(k(x sin \theta + y cos \theta) - \omega.t+ \varphi \bigg) \end{equation*}

$a$ 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.

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

$k$ is the wave number (in rad/m). The wavelength (in meters): $L = 2\pi/k$ 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:

$\omega^2= g.k.tanh(k.d)$

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

$\varphi$ 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:

$E = 1/2 \rho g a^2 = 1/8 \rho g H^2$

wherein:

g is the gravity acceleration (g ≈ 9.81 $m/s^2$)
$\rho$ is the water density (in $kg/m^3$) ($\rho \approx 1025 kg/m^3$ for seawater).

2.3 Random multidirectional waves

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

\begin{equation*} \eta(x,y,t)=\sum_{m=1}^M \eta_m(x,y,t)=\sum_{m=1}^M a_m.cos\bigg[k_m(x sin \theta_m +y cos \theta_m )-\omega_m.t+\varphi_m\bigg] \end{equation*}

A major point in the above expression concerns the phase distribution $\varphi_m$ of elementary wave components. The approach used in the TOMAWAC model assumes that these phases are randomly distributed over the [0; $2\pi$] 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:

\begin{equation*} \ E=\sum_{m=1}^M 1/2.\rho g a_m^2 \end{equation*}

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

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

2.4 Sea state directional power spectrum

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

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

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

\begin{equation*} \sum_f^{f+df}\sum_{\theta}^{\theta+d\theta} 1/2.\rho g a_m^2 = E(f,\theta)df.d\theta \end{equation*}

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

2.5 Directional spectrum of sea state variance

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

$F(f,\theta)= \frac{E(f,\theta)}{\rho g}$

Then, in particular, we have:

\begin{equation*} \sum_f^{f+df}\sum_{\theta}^{\theta+d\theta} 1/2.\rho g a_m^2 = E(f,\theta)df.d\theta \end{equation*}

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

\begin{equation*} \eta(x,y,t)=\int_{f=0}^{\infty}\int_{\theta=0}^{2\pi}\sqrt{2F(f,\theta)df.d\theta} .cos \bigg[k(xcos\theta+ysin\theta)-\omega.t+\varphi\bigg] \end{equation*}

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:

\begin{equation*} \ a_m = \sqrt{2F(f,\theta)df.d\theta} \end{equation*}

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

\begin{equation*} \ m_n\int_{f=0}^{\infty}\int_{\theta=0}^{2\pi}f^nF(f,\theta)df.d\theta \end{equation*}

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

\begin{equation*} \ <\eta^2> = \lim_{t_0 \rightarrow \infty} \frac{1}{t_0}\int_0^{t_0}\eta^2(t)dt = m_0=\int_{f=0}^{\infty}\int_{\theta=0}^{2\pi} F(f,\theta)df.d\theta \end{equation*}

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:

$h_{m0}=4\sqrt{m_0}$

The average frequencies $f_{01}$ and $f_{02}$ and $f_{-10}=f_e$ are also used and computed as follows:

\begin{equation*} \ f_{01}=\frac{m_1}{m_2}\qquad f_{02}=\sqrt{\frac{m_1}{m_2}}\qquad f_e=\frac{m_0}{m_{-1}} \end{equation*}

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

2.6 Sea state directional spectrum of wave action

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

That quantity, noted as $N(f,\theta)$, 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:

$N = F/\sigma $

wherein $\sigma$ 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 $\overrightarrow{U}$

$\Omega(\overrightarrow{k},\overrightarrow{x},t)=\omega=\sigma+\overrightarrow{k}.\overrightarrow{U}$

2.7 Selecting the directional spectrum discretization variables

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

$(f_a,\theta)$ = (absolute frequency; propagation direction)
$(fr,\theta)$ = (relative frequency; propagation direction)
$(k,\theta)$ = (wave number; propagation direction)
$(kx,ky) = (k.sin \theta ; k.cos \theta)$ = (wave number vector)
For the numerical resolution of equations, the model TOMAWAC uses the pair $(f_r,\theta)$ = (relative frequency; propagation direction)

The directional spectra output by TOMAWAC, however, are always expressed in $(f_a,\theta)$. The equations solved by TOMAWAC are thoroughly reviewed in section 4.


3. Hypotheses and application domain of TOMAWAC


3.1. Application domain of the model TOMAWAC

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

The application domain of TOMAWAC includes:

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

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

3.2. Wave interactions with other physical factors

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

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

The physical processes modelled in TOMAWAC

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

—> Energy source/dissipation processes:

—> Non-linear energy transfer conservative processes:

-—> Wave propagation-related processes:

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)


4. Mathematical modelling procedures used by TOMAWAC


4.1. Scope of sea state modelling

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

$N(\overrightarrow{x},\overrightarrow{k},t)=N(x,y,k_x,k_y,t)$

using, as discretization variables:

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

\begin{equation*} \frac{\partial N}{\partial t}+\frac{\partial (\dot{x}N)}{\partial x}+\frac{\partial(\dot{y}N)}{\partial y}+\frac{\partial (\dot{k_x}N)}{\partial k_x}+\frac{\partial (\dot{k_y}N)}{\partial k_y}= Q(k_x,k_y,x,y,t) \end{equation*}

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 :

\begin{equation*} \dot{g}=\frac{dg}{dt}=\frac{\partial g}{\partial t}+\frac{\partial x}{\partial t}\frac{\partial g}{\partial x}+\frac{\partial y}{\partial t}\frac{\partial g}{\partial y} \end{equation*}

In that form (conservative writing in the form of a flux), equation (4.1) can be transposed to other coordinate systems and, for instance, $(k, \theta)$, $(f_a,\theta)$ or else $(f_r,\theta)$ 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):

\begin{equation*} \dot{x}=\frac{\partial \Omega}{\partial k_x}\qquad\text{and}\qquad \dot{y}=\frac{\partial \Omega}{\partial k_y} \end{equation*}

\begin{equation*} \dot{k_x}=\frac{-\partial \Omega}{\partial x}\qquad\text{and}\qquad \dot{k_y}=\frac{-\partial \Omega}{\partial y} \end{equation*}

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

$\Omega(\overrightarrow{k},\overrightarrow{x},t)= \omega = \sigma + \vec{k}.\vec{U}$

wherein:

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

$f_a = \omega/(2\pi)$ is named absolute frequency.

$\vec{U}$ denotes the current velocity (depth-integrated).

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

$\sigma^2 = g.k.tanh(k.d)$

$fr = \sigma/(2\pi)$ 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:

\begin{equation*} \frac{\partial (\dot{x})}{\partial x}+\frac{\partial(\dot{y})}{\partial y}+\frac{\partial (\dot{k_x})}{\partial k_x}+\frac{\partial (\dot{k_y})}{\partial k_y}= 0 \end{equation*}

or $div(\vec{V}) = 0$ when defining : $V = (\dot{x},\dot{y},\dot{k_x},\dot{k_y})$

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

\begin{equation*} \frac{\partial N}{\partial t}+ \dot{x}\frac{\partial N}{\partial x}+\dot{y}\frac{\partial N}{\partial y}+\dot{k_x}\frac{\partial N}{\partial k_x}+\dot{k_y}\frac{\partial N}{\partial k_y}= Q(k_x,k_y,x,y,t) \end{equation*}

\begin{equation*} \frac{\partial N}{\partial t}+\vec{V}.grad_{\vec{x},\vec{k}}(N)= Q \end{equation*}

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

\begin{equation*} \begin{align} \dot{x} & = C_g \frac{k_x}{k}+U_x \\ \dot{y} & = C_g \frac{k_y}{k}+U_y \\ \dot{k_x} & = -\frac{\partial \sigma}{\partial d}\frac{\partial d}{\partial x}-\vec{k}.\frac{\partial\vec{U}}{\partial x} \\ \dot{k_y} & = -\frac{\partial \sigma}{\partial d}\frac{\partial d}{\partial y}-\vec{k}.\frac{\partial\vec{U}}{\partial y} \\ \end{align} \end{equation*}

$C_g$ 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:

\begin{equation*} \ C_g = \frac{\partial \sigma}{\partial k} = n\frac{\sigma}{k}\qquad \text{with}\qquad n= \frac{1}{2}\bigg(1+\frac{2kd}{sinh(2kd)}\bigg) \end{equation*}

The relative (or intrinsic) phase velocity C of waves is also introduced : $C=\frac{\sigma}{k}$

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 :

The following conventions are adopted for writing the equations :

coordinate system) is assumed to be horizontal, directed to the right, whereas the y-axis (in the Cartesian coordinate system) or the $\varphi$-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 the clockwise direction.

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


4.2. Equation solved

4.2.1. Equations solved in a Cartesian spatial coordinate system

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

\begin{equation*} \begin{align} \ N(x,y,k_x,k_y,t)= \frac{C.C_g}{2\pi\sigma}\tilde{N}(x,y,f_r,\theta,t)=\tilde{B}.\tilde{F}(x,y,f_r,\theta,t)\\ \text{putting :}\qquad\tilde{B}=\frac{C.C_g}{2\pi\sigma^2}=\frac{C_g}{(2\pi)^2k f_r}\\ \end{equation*}

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

\begin{equation*} \frac{\partial(\tilde{B}\tilde{F})}{\partial t}+\dot{x}\frac{\partial(\tilde{B}\tilde{F})}{\partial x}+\dot{y}\frac{\partial(\tilde{B}\tilde{F})}{\partial y}+\dot{\theta}\frac{\partial(\tilde{B}\tilde{F})}{\partial \theta}+\dot{f_r}\frac{\partial(\tilde{B}\tilde{F})}{\partial f_r}= \tilde{B}.\tilde{F}(x,y,f_r,\theta,t) \end{equation*}

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

\begin{equation*} \begin{align} \dot{x} & = C_g.sin\theta + U_x\\ \dot{y} & = C_g.sin\theta + U_y\\ \dot{\theta} & = -\frac{1}{k}\frac{\partial\sigma}{\partial d}\tild{G}_n(d)-\frac{\vec{k}}{k}.\tilde{G}_n(\vec{U})\\ \dot{f_r} & = \frac{1}{2\pi}\bigg[\frac{\partial\sigma}{\partial d}\bigg(\frac{\partial d}{\partial t}+\vec{U}.\vec{\nabla}d\bigg)- C_g \vec{k}.\vec{\tild{G}_n}(\vec{U})\bigg]\\ \end{align} \end{equation*}

The operators $\tilde{G}_n$ and $\tilde{G}_t$ refer to the computation of a function gradient in directions that are respectively normal and tangential to the characteristic curve with the direction $\theta$ :

\begin{equation*} \begin{align} \tilde{G}_n(g) & = \vec{n}.\vect{nabla}g = cos \theta\frac{\partial g}{\partial x} - sin\theta\frac{\partial g}{\partial y} \tilde{G}_t(g) & = \vec{t}.\vect{nabla}g = sin \theta\frac{\partial g}{\partial x} + cos\theta\frac{\partial g}{\partial y} \end{align} \end{equation*}

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

\begin{equation*} \frac{\partial\sigma}{\partial d} = \frac{\sigma k}{sinh(2kd)} \end{equation*}

The spatial transfer rates $\dot{x}$ and $\dot{y}$ (equations 4.12.a and 4.12.b) model the spatial wave propagation and the shoaling. The directional transfer rate $\dot{\theta}$ (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 $\dot{f_r}$ (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 $\tilde{B}$ allows to switch from the term Q to a term $\tilde{Q}$ that is directly expressed in terms of the directional variance spectrum with a variance $\tilde{F}(f_r,\theta)$. The content of that term is explained in paragraph 4.2.3.

4.2.2. Equations solved in a spherical spatial coordinate system

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

\begin{equation*} \N(x,y,k_x,k_y,t)=\frac{C.C_g}{2\pi\sigma^2 R^2 cos\varphi}=\frac{C_g}{(2\pi)^2 k f_r R^2 cos\varphi} \end{equation*}

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

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

\begin{equation*} \frac{\partial(\hat{B}\hat{F})}{\partial t}+\dot{\lambda}\frac{\partial(\hat{B}\hat{F})}{\partial \lambda}+\dot{\varphi}\frac{\partial(\hat{B}\hat{F})}{\partial \varphi}+\dot{\theta}\frac{\partial(\hat{B}\hat{F})}{\partial \theta}+\dot{f_r}\frac{\partial(\hat{B}\hat{F})}{\partial f_r}=\hat{B}.\hat{Q}(\lambda,\varphi,f_r,\theta,t) \end{equation*}

with the following transfer rates:

\begin{equation*} \begin{align} \dot{\lambda} & = \frac{1}{Rcos\varphi}(C_g.sin\theta + U_{\lambda})\\ \dot{\varphi} & = \frac{1}{R}(C_g.sin\theta + U_{\varphi})\\ \dot{\theta} & = \frac{1}{R}\bugg[C_g.sin\theta.tan\varphi-\frac{1}{k}\frac{\partial\sigma}{\partial d}\hat{G}_n(d)-\frac{\vec{k}}{k}.\hat{G}_n(\vec{U})\\ \dot{f_r} & = \frac{1}{2\pi R}\bigg[\frac{\partial\sigma}{\partial d}\bigg(\frac{\partial d}{\partial t}+\frac{U_{\lambda}}{cos\varphi}\frac{\partial d}{\partial \lambda} + U_{\varphi}\frac{\partial d}{\partial \varphi}\bigg)- C_g \vec{k}.\hat{G}_t(\vec{U})\bigg]\\ \end{align} \end{equation*}

As in the previous case, the operators $\hat{G}_n$ nd $\hat{G}_t$ refer to the computation of a function gradient in directions that are respectively normal and tangential to the characteristic curve with the direction $\theta$ :

\begin{equation*} \begin{align} \hat{G}_n(g) & =\frac{cos\theta}{cos\varphi}\frac{\partial d}{\partial \lambda}-sin\theta \frac{\partial g}{\partial \varphi}\\ \hat{G}_n(g) & =\frac{sin\theta}{cos\varphi}\frac{\partial d}{\partial \lambda}-cos\theta \frac{\partial g}{\partial \varphi}\\ \end{align} \end{equation*}

As previously, the spatial transfer rates $\dot{\lambda}$ and $\dot{\varphi}$ (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 $\dot{\theta}$ 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 $\dot{\theta}$ (equation 4.18.c) model the refraction caused respectively by bathymetry and currents. The relative frequency transfer rate $\dot{f_r}$ (equation 4.18.d) models the changes of relative frequency resulting from variations of the sea level or of the current in both space and time. It is noteworthy that this last term is zero in the case of zero current and of no variation of the sea level in time: the advection equation is then reduced to a three-dimensional equation.

4.2.3. TOMAWAC source and sink terms

4.2.3.1. Generals

The source and sink terms that compose $\tilde{Q}$ and $\hat{Q}$ 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:

\begin{equation*} \ Q = Q_{in}+ Q_{ds}+ Q_{nl}+ Q_{bf}+ Q_{br}+ Q_{tr} \end{equation*}

wherein :

$Q_{in}$: wind-driven wave generation $Q_{ds}$: whitecapping-induced energy dissipation $Q_{nl}$: non-linear quadruplet interactions $Q_{bf}$: bottom friction-induced energy dissipation $Q_{br}$: bathymetric breaking-induced energy dissipation $Q_{tr}$: 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.

4.2.3.2. Wind input (term Qin)

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

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

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

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

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

\begin{equation*} \ Q_{in} = \sigma\epsilon\beta\bigg((\frac{u*}{C}+ z_{\alpha})max(cos (\theta-\theta_w);0)\bigg)^2 F \end{equation*}

with the following notations:

$u*=\sqrt{\frac{\tau_s}{\rho_{air}}}$

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

It is written as $\beta = \frac{\beta_m}{\kappa}\mu ln^4 \mu$

where

\begin{equation*} \mu = min\bigg[\frac{g.z_0}{C^2}.exp\bigg(\frac{\kappa}{\Big[\frac{u*}{C}+z_{alpha}\Big]cos(\theta-\theta_w)}\bigg);1\bigg] = min \bigg[\Omega.A^2exp\bigg(\frac{\kappa}{\Big[A + z_{\alpha}\Big]cos(\theta-\theta_w)}\bigg);1\bigg] \end{equation*}

The Janssen’s model [Janssen, 1989] [Janssen, 1991] is characterized by the method it uses for computing $u*$ and $z_0$. The surface stress $\tau_s$ 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 $\tau_s$. It is obtained by solving the following system of equations:

\begin{equation*} \begin{align} \ U_{10} & = \frac{u*}{\kappa}ln\big(\frac{10+z_0+\tilde{z_0}}{z_0}\big)\approx\frac{u*}{\kappa}ln\big(\frac{10}{z_0}\big)\\ \ z_0 & = \frac{\tilde{z_0}}{\sqrt{1-\tau_w/\tau_s}}\\ \tilde{z_0} & = \alpha\frac{u*^2}{g}\\ \ u* & = \sqrt{\frac{\tau_s}{\rho_{air}}} \end{align} \end{equation*}

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

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

$u* = \sqrt{C_D}U_{10}$ where: $C_D = 1.2875.10^{-3}$ by default.

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

\begin{equation*} \tau_w = \int\int\rho_{water}\sigma Q_{in}(f_r,\theta)df_r d\theta \end{equation*}

That integral is numerically computed over the discretized portion of the spectrum and a parametrization, based upon a decrement of variance in $f^{-n}$, is used for the high frequencies portion of the spectrum.

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

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

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

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

\begin{equation*} \ Q_{in} = \beta F \qquad\text{where:}\qquad \beta = max\bigg[0;0.25\frac{\rho_{air}}{\rho_{water}}\Big(28\frac{u*}{C} cos(\theta - \theta_w)-1\Big)\bigg]\sigma \end{equation*}

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

$u* = \sqrt{C_D}U_{10}$ where:

\begin{equation*} \begin{align} \ C_D = 6.5.10^{-5}U_{10}+8.10^{-4} \qquad\text{if}\qquad U_{10}>7.5m/s\\ \ C_D = 1.28875.10^{-3} \qquad\text{if}\qquad U_{10}<7.5m/s\\ \end{align} \end{equation*}

That source term only uses two parameters, namely:

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

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

\begin{equation*} \ Q_{in} = \beta F \qquad\text{where:}\qquad \beta = \bigg[D\Big(\frac{u*}{C}\Big)^{2} cos(\theta - \theta_w)+E\Big(\frac{u*}{C}\Big)cos(\theta - \theta_w)+F.cos(\theta - \theta_w)+H\bigg]\sigma \end{equation*}

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

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

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

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

\begin{equation*} \ Q_m(f,\theta) = \alpha(f,\theta) = 1.5.10^{-3}g^{-2}\Big[u_* max(0,\theta - \theta_w)\Big]^4.exp\bigg[-\Big(\frac{f}{f_{PM}}\Big)^{-4}\bigg] \end{equation*}

where $u_*$ is the friction wind velocity, $\theta_w$ the wind direction and $f_{PM}$ is a peak frequency called Pierson- Moskowitz frequency [Pierson & Moskowitz, 1964], defined as:

\begin{equation*} \ f_{PM} = \frac{1}{2\pi}\frac{g}{28u_*} \end{equation*}

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

4.2.3.3. Whitecapping-induced dissipations (term Qds)

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

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

4.2.3.3.1. option 1 for whitecapping: Komen and Janssen dissipation model

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

\begin{equation*} \ Q_{ds} = -\frac{1}{2g^4} C_{dis}\bar{\sigma}^9 m_0^2\Big(\delta\frac{\sigma}{\bar{\sigma}}+(1-\delta)\big(\frac{\sigma}{\bar{\sigma}}\big)^2\Big)F \end{equation*}

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

\begin{equation*} \ Q_{ds} = -\frac{1}{2} C_{dis}\bar{\sigma}\bar{k} m_0^2\Big(\delta\frac{k}{\bar{k}}+(1-\delta)\big(\frac{k}{\bar{k}}\big)^2\Big)F \end{equation*}

$m_0$ denotes the total variance, $\bar{\sigma}$ denotes the average intrinsic frequency and $\bar{\sigma}$ denotes the average wave number; they are respectively computed as followings:

\begin{equation*} \ m_0= \int_{f_r=0}^{\infty} \int_{\theta=0}^{2\pi} F(f_r,\theta)df_rd\theta \end{equation*}

\begin{equation*} \bar{\sigma}= \bigg(\frac{1}{m_0}\int_{f_r=0}^{\infty} \int_{\theta=0}^{2\pi} \frac{1}{\sigma}F(f_r,\theta)df_rd\theta\bigg)^{-1} \end{equation*}

\begin{equation*} \bar{k}= \bigg(\frac{1}{m_0}\int_{f_r=0}^{\infty} \int_{\theta=0}^{2\pi} \frac{1}{\sqrt{k}}F(f_r,\theta)df_rd\theta\bigg)^{-2} \end{equation*}

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 +$\infty$) considering a decreasing variance in $f^{-n}$.

That source term has two parameters:

4.2.3.3.2. Option 2 for whitecapping: Westhuysen dissipation model

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

The expression proposed by Westhuysen is:

\begin{equation*} \ Q_{ds} = -C_{dis,break}\bigg(\frac{B(k)}{B_r}\bigg)^{\frac{p_0}{2}} g^{\frac{1}{2}}k^{\frac{1}{2}} F(f,\theta) \end{equation*}

where:

\begin{equation*} \ B(k) = \frac{1}{2\pi}\int_{0}^{2\pi} C_g k^3 F(f,\theta)d\theta = C_g k^3\frac{E(f)}{2\pi} \end{equation*}

and

\begin{equation*} \ p_0\frac{u_*}{C} = 3 + tanh \Big[w\big(\frac{u_*}{C}-0.1\big)\Big] \end{equation*}

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:

\begin{equation*} \ Q_{ds} = f_{br}(f).Q_{ds}^W + \big(1-f_{br}(f)\big).Q_{ds}^k \end{equation*}

\begin{equation*} \ f_{br} = \frac{1}{2}+\frac{1}{2}tanh\Bigg(10\bigg[\bigg(\frac{B(k)}{B_r}\bigg)^{1/2}-1\bigg]\Bigg) \end{equation*}

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

4.2.3.4. Bottom friction-induced dissipations (term Qbf)

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

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

\begin{equation*} \ Q_{bf} = -\Gamma\bigg(\frac{\sigma}{g.sinh(kd)}\bigg)^2F \end{equation*}

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

\begin{equation*} \ Q_{bf} = -\Gamma\frac{2k}{g.sinh(2kd)}F \end{equation*}

That source term has a single parameter:

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

4.2.3.5. Bathymetric breaking-induced dissipations (term Qbr)

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

0 No breaking-induced dissipation (default value)
1 Battjes and Janssen’s model [Battjes, 1978]
2 Thornton and Guza’s model [Thornton, 1983]
3 Roelvink’s model [Roelvink, 1993]
4 Izumiya and Horikawa’s model [Izumiya, 1984]

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

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

4.2.3.5.1.Battjes and Janssen’s model (1978)

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

\begin{equation*} \ D_{br} = -\frac{\alpha Q_b f_c H_m^2}{4} \end{equation*}

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

\begin{equation*} \ H_m = \gamma_2 d \end{equation*}

or through a relation derived from the Miche’s criterion

\begin{equation*} \ H_m = \frac{\gamma_1}{k_c}.tanh\bigg(\frac{\gamma_2 k_c d}{\gamma_1}\bigg) \end{equation*}

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

\begin{equation*} \frac{1-Q_b}{ln Q_b} = -\frac{H_{m0}}{2H_m^2} \end{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:

\begin{equation*} \ b= \frac{H_{m0}}{\sqrt{2H_m^2}} \end{equation*}

version 1:

\begin{equation*} \begin{align} \ Q_b & = 0 \qquad & \text{if}\qquad b < C_b , (Cb = 0.5)\\ \ Q_b & = \bigg(\frac{b-C_b}{1-C_b}\bigg)^2 \qquad & \text{if}\qquad b > C_b\\ \end{align} \end{equation*}

version 2:

\begin{equation*} \begin{align} \ q_0 & = (2b-1)^2 \qquad & \text{if}\qquad 0.5 < b < 1\\ \ q_0 & = 0 \qquad & \text{if}\qquad b \le 0.5\\ \ q_1 & = q_0 - b^2 \Bigg(\frac{q_0 - e^{\frac{q_0-1}{b^2}}}{b^2 - e^{\frac{q_0-1}{b^2}}}\Bigg)\\ \ Q_b & = 0 \qquad & \text{if}\qquad b \le C_b , (Cb = 0.3)\\ \ Q_b & = q_1 \qquad & \text{if}\qquad C_b < b < 0.9\\ \ Q_b & = q_0 \qquad & \text{if}\qquad 0.5 \le b \le 1.0\\ \end{align} \end{equation*}

version 2:

\begin{equation*} \ Q_b = 2.4 b^7 \end{equation*}

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:

\begin{equation*} \ Q_{br}(f,\theta) = -\frac{\alpha Q_b f_c H_m^2}{4}\frac{F(f,\theta)}{m_0} \end{equation*}

Three constants can be modified using keywords:

The following keywords are for selecting the model options:

4.2.3.5.2. Thornton and Guza’s model (1983)

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

function 1:

\begin{equation*} \ Q_{br1}(f,\theta) = -48\sqrt{\pi}B^3 f_c \frac{(2m_0)^{5/2}}{H_m^4 d} F(f,\theta) \end{equation*}

function 2:

\begin{equation*} \ Q_{br2}(f,\theta) = -12\sqrt{\pi}B^3 f_c \frac{(2m_0)^{3/2}}{H_m^4 d}\Bigg[1-\bigg(1+\bigg(\frac{8m_0}{H_m^2}\bigg)\bigg)^{-5/2}\Bigg]F(f,\theta) \end{equation*}

$f_c$ is the characteristic wave frequency (average frequency, $f_{01}$ , $f_{02}$ 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, $H_m$ , is governed by the parameter $\gamma (H_m = \gamma d )$.

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

4.2.3.5.3.Roelvink’s model (1993)

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

Weibull’s distribution:

\begin{equation*} \ Q_{br1}(f,\theta) = -\alpha f_c mA\sqrt{\frac{2}{m_0}}F(f,\theta)\int_0^{\infty}\bigg(\frac{H}{\sqrt{8m_0}}\bigg)^{2m+1}exp\Bigg[-A\bigg(\frac{H}{\sqrt{8m_0}}\bigg)^{2m}\Bigg]\Bigg[1-exp\bigg(-\bigg(\frac{H}{\gamma d}\bigg)\bigg)^N\Bigg].dH \end{equation*}

\begin{equation*} \ A = \bigg[\Gamma\bigg(1+\frac{1}{m}\bigg)\bigg]^m \qquad\tex{with}\qquad m = 1+0.7 tan^2\bigg(\frac{\pi}{2}\frac{1}{\gamma_2}\frac{\sqrt{8m_0}}{d}\bigg) \end{equation*}

The coefficient $\gamma_2$ is usually set to 0.65.

Rayleigh’s distribution:

\begin{equation*} \ Q_{br2}(f,\theta) = -\alpha f_c mA\sqrt{\frac{2}{m_0}}F(f,\theta)\int_0^{\infty}\bigg(\frac{H}{\sqrt{8m_0}}\bigg)^{3}exp\Bigg[-A\bigg(\frac{H}{\sqrt{8m_0}}\bigg)^{2}\Bigg]\Bigg[1-exp\bigg(-\bigg(\frac{H}{\gamma d}\bigg)\bigg)^N\Bigg].dH \end{equation*}

$f_c$ denotes the characteristic wave frequency (average frequency, $f_{01}$ , $f_{02}$ or peak frequency), $\alpha$ is a numerical constant of order 1, $\gamma$ is the proportional control factor between the allowable wave height and the water depth (by default, $\gamma = 0.54$ ) 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:

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

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

\begin{equation*} \frac{d}{dx}(EC_g) = -\alpha \frac{E^{3/2}}{rho^{1/2}d^{3/2}}\bigg(\frac{2C_g}{c}-1\bigg)^{1/2} \end{equation*}

E denotes the total wave energy, Cg and c are respectively group and phase velocities associated to the characteristic wave frequency $f_c$ (average frequency $f_{01}$ , $f_{02}$ or peak frequency), $\alpha$ 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 $\alpha$ in terms of a deviation from a steady state:

\begin{equation*} \alpha = \beta_0(M_*^2 - M_{*S}^2)^{1/2} \end{equation*}

where $M_*$ is a dimensionless quantity in the form:

\begin{equation*} \ M_*=\frac{C_g}{c}\frac{E}{\rho g d^2} \end{equation*}

From laboratory data, Izumiya and Horikawa set $M_{*S}^2$ to 9 10-3 and $\beta_0$ to 1.8.

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

\begin{equation*} \ Q_{br}(f,\theta) = -\beta_0\bigg(\frac{C_g}{c}\frac{m_0}{d^2} -M_{*S}^2\bigg)^{1/2}\frac{g^{1/2}m_0^{1/2}}{d^{3/2}}\bigg(\frac{2C_g}{c}-1\bigg)^{1/2}F(f,\theta) \end{equation*}

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

4.2.3.6. Non-linear quadruplet interactions (term Qnl)

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

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].
4.2.3.6.1.Option 1 for non-linear quadruplet interactions: DIA method

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

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

\begin{equation*} \begin{align} \ Q_{nl}^{exact} & = \int \int \int \sigma_4.G.\delta (\vec{k_1}+\vec{k_2}-\vec{k_3}-\vec{k_4})\delta(\sigma_1+\sigma_2-\sigma_3-\sigma_4)\\ \ & \bigg[\frac{F(\vec{k_1})}{\sigma_1}\frac{F(\vec{k_2})}{\sigma_2}\bigg(\frac{F(\vec{k_3})}{\sigma_3}+\frac{F(\vec{k_4})}{\sigma_4}\bigg)-\frac{F(\vec{k_3})}{\sigma_3}\frac{F(\vec{k_4})}{\sigma_4}\bigg(\frac{F(\vec{k_1})}{\sigma_1}+\frac{F(\vec{k_2})}{\sigma_2}\bigg)\bigg]d\vec{k_1}.d\vec{k_2}.d\vec{k_3}\\ \end{align} \end{equation*}

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

\begin{equation*} \sigma_1+\sigma_2 = \sigma_3 +\sigma_4 \qquad\text{and}\qquad \vec{k_1}+\vec{k_2}= \vec{k_3}+\vec{k_4} \end{equation*}

as evidenced by the two Dirac functions $\delta$ in the integral.

G denotes the value of the coupling term for the resonant quadruplet interactions $(\vec{k_1},\vec{k_2},\vec{k_3},\vec{k_4})$

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:

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

\begin{equation*} \begin{bmatrix} Q_{nl}\\ Q_{nl}^-\\ Q_{nl}^+\\ \end{bmatrix} = \begin{bmatrix} -2\\ 1\\ 1\\ \end{bmatrix} \Pi g^{-4}f_r^{11}\Bigg[F^2\bigg(\frac{F^+}{(1+\lambda)^4}+\frac{F^-}{(1-\lambda)^4}\bigg)-\frac{FF^+F^-}{(1-\lambda^2)^4}\Bigg] \end{bmatrix} \end{equation*}

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:

\begin{equation*} \ Q_{nl}(d)= R.Q_{nl}(d= \infty) \end{equation*}

Coefficient R is a function of the normalized water height $\bar{k}.d$ and is expressed as follows:

\begin{equation*} \ R(\chi) = 1+ \frac{5.5}{\chi}\big(1-\frac{5}{6}\chi\big)exp\big(-\frac{5}{4}\chi\big)\qquad\text{where}\qquad \chi = \frac{3}{4}\bar{k}.d \end{equation*}

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

That source-term has a single parameter:

4.2.3.6.2.Option 2 for non-linear quadruplet interactions: MDIA method

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:

\begin{equation*} \begin{align} \vec{k_0} + \vec{k_1} & = \vec{k_2} +\vec{k_3}= 2\vec{k}\\ \sigma_0 & = (1+\mu)\sigma\\ \sigma_1 & = (1-\mu)\sigma\\ \sigma_2 & = (1+\lambda)\sigma\\ \sigma_3 & = (1-\lambda)\sigma\\ \end{align} \end{equation*}

The $\lambda$ and $\mu$ values proposed by Tolman that allow to best estimate the $Q_{nl4}$ 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.

4.2.3.6.3. Option 3 for non-linear quadruplets interactions: GQM method

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

\begin{equation*} \begin{align} \sigma_a & = \sigma_0 + \sigma_1\\ \vec{k_a} & = \vec{k_0} + \vec{k_1}\\ \epsilon_a & = \frac{2gk_a}{\sigma_a^2} \end{align} \end{equation*}

and depending only on ($\sigma_0$ et $\theta_0$) and ($\sigma_1$ et $\theta_1$), are used as well.

\begin{equation*} \ Q_{nl4}= \int_{\sigma_1=0}^{+\infty}\int_{\theta_1=0}^{2\pi}\int_{\sigma_2=0}^{\sigma_a / 2}2\frac{\sigma_a^4 T}{\sigma_1\sigma_2\sigma_3}\frac{F_2F_3(F_0\sigma_1^4+F_1\sigma_0^4)-F_0F_1(F_2\sigma_3^4+F_3\sigma_2^4)}{\sqrt{\tilde{B_0}(\epsilon_a,s_2)\tilde{B_1}(\epsilon_a,s_2)\tilde{B_2}(\epsilon_a,s_2)}}d\sigma_1 d\theta_1 d\sigma_2 \end{equation*}

where $s_2$ is defined as $s2= \sigma_2 / \sigma_a$.

Equation 4.55 is then integrated using different quadrature methods:

Three different GQM method resolutions have been tested:

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:

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:

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


4.2.3.7.1.LTA (Lumped Triad Approximation) model

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

\begin{equation*} \begin{align} \ Q_{LTA}(f,\theta) & = Q_{LTA}^+(f,\theta)+Q_{LTA}^-(f,\theta)\\ \ Q_{LTA}^+(f,\theta) & = \alpha_{LTA} cC_g g^2 \frac{R^2_{(f/2,f/2)}}{S_f^2} sin |\beta_{f/2,f/2}|\bigg[F^2(f/2,0)-2F(f,\theta)F(f/2,0)\bigg]\\ \ Q_{LTA}^-(f,\theta) & = -2Q_{LTA}^+(2f,\theta)\\ \end{align} \end{equation*}

$\alpha_{LTA}$ is the model adjustment coefficient; c and Cg denote the phase and group velocities, respectively.

R is the self-interaction coefficient

\begin{equation*} \ R_{f,f} = (2k)^2\bigg[\frac{1}{2}+\frac{(2\pi f)^2}{gdk^2}\bigg] \end{equation*}

\begin{equation*} \ S_f = -2k\bigg[gd + 2Bgd^3k^2-(B+\frac{1}{3}d^2)(2\pi f)^2\bigg] \end{equation*}

The biphase $\beta$ is given by the relation

\begin{equation*} \beta(f,f) = -\frac{\pi}{2}+\frac{\pi}{2} tanh \bigg(\frac{0.2}{Ur}\bigg) \end{equation*}

where Ur denotes the Ursell's number

\begin{equation*} \ Ur = \frac{g}{8\pi^2\sqrt{2}}\frac{H_{m0}T_m^2}{d^2} \end{equation*}

with $H_{m0}$ being the significant spectral height and $T_m$ being the average wave time. $Q_{LTA}^{\pm}$ denotes the negative and positive contributions of the self-interactions. Since $Q_{LTA}^+$ denotes the positive contributions to the first upper harmonic, it should be positive. The negative values of $Q_{LTA}^+$ 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 $R_{fm} = 2.5$.

Two constants can me modified through keywords:

4.2.3.7.2. SPB model

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

\begin{equation*} \begin{align} \ Q(f,\theta) = & \frac{B'g}{2S_{1,f}}\int_0^{f}\int_0^{2\pi}\int_0^{f}\int_0^{2\pi} df_1 df_2 d\theta_1 d\theta_2 T_{f_1,f_2} \delta(\theta_{\vec{k}}-\theta_{\vec{k_1}+\vec{k_2}})\delta(f-f_1-f_2)\\ \ & +\frac{B'g}{S_{1,f}}\int_0^{\infty}\int_0^{2\pi}\int_0^{\infty}\int_0^{2\pi} df_1 df_2 d\theta_1 d\theta_2 T_{-f_2,f_1} \delta(\theta_{\vec{k_1}}-\theta_{\vec{k}+\vec{k_2}})\delta(f_1-f-f_2)\\ \end{align} \end{equation*}

\begin{equation*} \begin{align} \ T_{f_1,f_2} & = \frac{gK}{K^2+\Delta k^2}R_{f_1,f_2}\bigg[-\frac{R_{-f_2,f}}{S_{2,f_1}k_1}FF_2-\frac{R_{-f_1,f}}{S_{2,f_2}k_2}FF_1+\frac{R_{f_1,f_2}}{S_{2,f}k}F_1F_2\bigg]\\ \ B' & = \frac{Cg}{2\pi k}\\ \ R_{f_1,f_2} & = (k_1+k_2)^2\bigg[\frac{1}{2}+\frac{(2\pi)^2 f_1f_2}{gdk_1k_2}\bigg]\\ \ S_f & = -2k\bigg(gd + 2Bgd^3k^2 - (B+\frac{1}{3})d^2(2\pi f)^2\bigg)\\ \end{align} \end{equation*}

F denotes the variance spectrum in terms of frequencies and directions, $F(f ,\theta)$. $T_{f_1,f_2}$ and $T_{-f_2,f_1}$ 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:


5. Discretizations used in TOMAWAC


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

5.1. Spatial discretization

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

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

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

5.2. Spectro-angular discretization

5.2.1. Frequency discretization

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

$f_n = f_1.q^{n-1}$ with n ranging from 1 to NF

The minimum frequency is then f1 and the maximum frequency is $f_1.q^{NF-1}$.

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

5.2.2. Directional discretization

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

$\theta_m = (m-1).360/ND$ 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:

5.2.3. Spectro-angular grid

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

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

5.3. Temporal discretization

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

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


6. Numerical methods used in TOMAWAC


6.1. General solution algorithm

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

\begin{equation*} \frac{\partial(BF)}{\partial t} + \vec{V}.\vec{\nabla}(BF) = BQ \end{equation*}

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

\begin{equation*} \vec{V} = \begin{bmatrix} \dot{x}\\ \dot{y}\\ \dot{\theta}\\ \dot{f_r}\\ \end{bmatrix} \end{equation*}

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 $t = n.\Delta t$, in which the variance spectrum $F^n$ is known in all points:

\begin{equation*} \frac{\partial(BF)}{\partial t} + \vec{V}.\vec{\nabla}(BF) = 0 \end{equation*}

discretized as follows:

\begin{equation*} \frac{(BF)^*-(BF)^n}{\Delta t} = [\vec{V}.grad(BF)]^n \end{equation*}

from which a value of $(B.F)^*$, then of $F^*$, intermediate after the convection step, is derived

\begin{equation*} \frac{\partial F}{\partial t} = Q \end{equation*}

discretized as follows:

\begin{equation*} \frac{F^{n+1}-F^*}{\Delta t} = \frac{Q^{n+1}-Q^*}{2} \end{equation*}

since coefficient B is time independent.

The variance density spectrum $F^{n+1}$ for a time step (time $t = (n+1).\Delta t$) is then obtained. That operation is then repeated for the next time step and as many times as necessary for covering the simulation period being considered.

6.2. Processing the propagation step

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

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

\begin{equation*} \frac{(BF)^*-(BF)^n}{\Delta tp} = [\vec{V}.grad(BF)]^n \end{equation*}

The convector field $\vec{V}$, 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:

\begin{equation*} \frac{BF^*-BF^n}{\Delta tp} = [\vec{V}.grad(BF)]^n \end{equation*}

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 $(x_Q, y_Q, \theta_Q, fr_Q)$ of the discretized spatial and spectro-angular variables, the characteristic curve is traced back to the time step $\Delta tp$ and the “arrival” point $(x_P, y_P, \theta_P, fr_P)$, 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:

$(B.F)* (x_Q, y_Q, \theta_Q, fr_Q) = (B.F)^n(x_P, y_P, \theta_P, fr_P)$

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

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

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

6.3. Processing the source term integration step

6.3.1. Source term integration numerical scheme

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

\begin{equation*} \frac{F^{n+1}-F^*}{\Delta t} = \frac{Q^{n+1}-Q^*}{2} \end{equation*}

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: $\Delta F = F^{n+1} – F^*$ and the source and sink terms are classified as linear and non-linear terms in F:

Q = Qnl + Ql
$Ql^{n+1} = \beta^{n+1}.F^{n+1} = \beta^{n+1}.F^* + \beta^{n+1}.\Delta F$

\begin{equation*} \ Q_{nl}^{n+1} \approx Q_{nl}^* + \frac{\partial Q_{nl}^*}{\partial F} \Delta F \end{equation*}

$\frac{\partial Q_{nl}^*}{\partial F}$ is a matrix of differential increments that is broken down into a diagonal portion [L*] and an extradiagonal portion $[N^*]$:

\begin{equation*} \frac{\partial Q_{nl}^*}{\partial F} = [M^*] = [\Lambda^*] = [N^*] \end{equation*}

substituting into the expression of $Q_{nl}^{n+1}$, this yields:

\begin{equation*} \ Q_{nl}^{n+1} \approx Q_{nl}^* + \bigg([\Lambda^*]+[N^*]\bigg)\Delta F \end{equation*}

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

\begin{equation*} \begin{align} \ Q^* & = \beta^*F^* + Q_{nl}^*\\ \ Q^{n+1} & = Q_{l}^{n+1}+Q_{nl}^{n+1} = \beta^{n+1}F^*+\beta^{n+1}\Delta F+Q_{nl}^* + \bigg([\Lambda^*]+[N^*]\bigg)\Delta F\\ \end{align} \end{equation*}

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

\begin{equation*} \Delta F = F^{n+1}-F^* = \frac{\Delta t}{2}(Q^{n+1}+Q^*) \end{equation*}

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

\begin{equation*} \Delta F = \frac{\Delta t}{2}\bigg[\beta^*F^*+ Q_{nl}^* + \beta^{n+1}F^*+\beta^{n+1}\Delta F+Q_{nl}^* + \bigg([\Lambda^*]+[N^*]\bigg)\Delta F\bigg] \end{equation*}

\begin{equation*} \Delta F \bigg(1+\frac{\Delta t}{2}\bigg[\beta^{n+1}+ \Big([\Lambda^*]+[N^*]\Big)\bigg]\bigg) = \Delta t \bigg( \frac{\beta^*+\beta^{n+1}}{2}F^*+Q_{nl}^*\bigg) \end{equation*}

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:

\begin{equation*} \Delta F = \frac{ \Delta t \bigg( \frac{\beta^*+\beta^{n+1}}{2}F^*+Q_{nl}^*\bigg)}{1+\frac{\Delta t}{2}\bigg[\beta^{n+1}+[\Lambda^*]\bigg]} \end{equation*}

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

\begin{equation*} \Delta F = \frac{\Delta t.Q_{TOT}}{1+\frac{\Delta t}{2} QDER} \end{equation*}

where :

\begin{equation*} \ Q_{TOT} = \frac{\beta^*+\beta^{n+1}}{2}F^*+Q_{nl}^* \end{equation*}

denotes a total source term, and

\begin{equation*} \ QDER = \beta^{n+1}+[\Lambda^*] \end{equation*}

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:

\begin{equation*} \frac{F^{n+1}-F^*}{\Delta t_2}=Q^*\qquad\text{i.e.}\qquad \Delta F=\Delta t_2 Q^* \end{equation*}

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

6.3.2. Monitoring the growth of the wave spectrum

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

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

\begin{equation*} \Delta F_{lim} = 0.62 10^{-4} \frac{\Delta}{1200} f^{-5} \end{equation*}

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

6.4. Processing the boundaries – Boundary conditions


6.4.1. Spatial grid

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

6.4.2. Spectro-angular grid

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

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


7. Inputs-outputs


7.1. Preliminary remark

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

The input files are:

The output files are:

7.2. Preliminary remark

7.2.1. The steering (or CAS) file

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

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

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

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

The write rules are as follows:

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:

7.2.2. The geometry file

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

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

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

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

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

7.2.3. The boundary conditions file

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

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

This file is described in detail in 8.5.1.

7.2.4. The currents file

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

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

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

7.2.5. The tidal water level file

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

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

7.2.6. The winds file

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

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

7.2.7. The previous computation file

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

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

7.2.8. The global results file

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

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

7.2.9. The 2D results file

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

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

1- all the data about the mesh geometry;
2- the names of the variables being stored;
3- for each time step, the time and the values of the variables are given for each point of the 2D mesh.

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

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

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

7.2.10.The punctual or spectrum results file

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

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

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

7.2.11.The printout listing

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

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

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

7.2.12.The User FORTRAN file

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

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

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

7.2.13.The auxiliary files

Other input/output files may be used by TOMAWAC.

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

7.2.14.The dictionary file

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

7.2.15.The libraries

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

The following libraries are used during the link editing operation:

7.3. Binary files

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

The following keywords can be used:

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

7.4. Files standard

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

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

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

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

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

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

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

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

7.5. Bathymetry data

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

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

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


8. Controlling the simulation


8.1. General parameterisation

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

8.1.1. Computation title

The computation case title is specified by the keyword TITLE

8.1.2. Computation time

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

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

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

8.1.3. Spectral discretisation

The spectral discretisation is defined by the following 5 keywords:

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:

$f_k = f_0 r^{k-1}$ (k = 1,NF)

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

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

8.1.4. Release

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

8.1.5. Environment

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

8.2. Computation options

8.2.1. Co-ordinate system

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

8.2.2. Finite depth

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

8.2.3. Taking a stationary current into account

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

The current can be specified in various ways:

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

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

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

8.2.4. Taking a wind into account

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

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

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

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

8.2.5. Recovering a TELEMAC data item

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

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

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

8.2.6. Taking the tide into account

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

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

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

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

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

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

8.2.7. Waves-current interactions: direct coupling with Telemac (2D or 3D) flow simulation

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:

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:

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

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

8.2.8. Convection step

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

8.3. Parameterising the source term integration step

8.3.1. Introduction

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

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

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

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

8.3.2.1. Wind input

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

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

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

AIR DENSITY,
WATER DENSITY,
WIND GENERATION COEFFICIENT,
VON KARMAN CONSTANT,
CHARNOCK CONSTANT,
SHIFT GROWING CURVE DUE TO WIND,
WIND DRAG COEFFICIENT,
WIND MEASUREMENT LEVEL.

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

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

AIR DENSITY,
WATER DENSITY,

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

YAN GENERATION COEFFICIENT D
YAN GENERATION COEFFICIENT E
YAN GENERATION COEFFICIENT F
YAN GENERATION COEFFICIENT H

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

8.3.2.2. White capping dissipation

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

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

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

WHITE CAPPING DISSIPATION COEFFICIENT
WHITE CAPPING WEIGHTING COEFFICIENT.

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

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

WESTHUYSEN DISSIPATION COEFFICIENT'
SATURATION THRESHOLD FOR THE DISSIPATION
WESTHUYSEN WHITE CAPPING DISSIPATION
WESTHUYSEN WEIGHTING COEFFICIENT'

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

8.3.2.3. Bottom friction dissipation

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

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

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

8.3.2.4. Non-linear transfers between quadruplets

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

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

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

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

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

SETTING FOR INTEGRATION ON OMEGA1
SETTING FOR INTEGRATION ON THETA1
SETTING FOR INTEGRATION ON OMEGA2
THRESHOLD0 FOR CONFIGURATIONS ELIMINATION
THRESHOLD1 FOR CONFIGURATIONS ELIMINATION
THRESHOLD2 FOR CONFIGURATIONS ELIMINATION

8.3.3. Source/sink terms in shallow water depth

8.3.3.1. Time integration scheme and time step

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

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

$\delta t_{i+1} = q \delta ti$

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 $H_{m0}/D$ ratio (D is the water depth) to 1. This ratio can be modified by means of the keyword MAXIMUM VALUE OF THE RATIO HM0 TO D. However, this is generally not advisable.

8.3.3.2. Wave breaking dissipation

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

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

1: Battjes and Janssen's model (1978)

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

2: Thornton and Guza's model (1983)

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

3: Roelvink's model (1993)

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

4: Izumiya and Horikawa's model (1984)

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

8.3.3.3. Triad interactions

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

Two formulations (see section 4.2.3.7) have been implemented.

1: LTA model

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

2: SPB model

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

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

8.4. Prescribing the initial conditions

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

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

$F(f,\theta)= E(f)D(\theta)$

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

\begin{equation*} \ E(f)= \alpha_{phil}H_{m0}^2 \frac{f_p^4}{f^5}.exp\bigg[-\frac{5}{4}\bigg(\frac{f_p}{f}\bigg)^4\bigg]\gamma^{exp\bigg[\frac{(f-f_p)^2}{2\sigma^2 f_p^2}\bigg]} \end{equation*}

where :

\begin{equation*} \begin{align} \sigma & = \sigma_a\qquad\text{for}\qquad f<f_p\\ \sigma & = \sigma_b\qquad\text{for}\qquad f>f_p\\ \alpha_{phil} & = \frac{0.0624}{0.23+0.0336\gamma - \frac{0.185}{1.9+\gamma}} \end{align} \end{equation*}

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 ( $\theta_1$ and $\theta_2$) as well as the weighting factor ($\lambda$ ) between the two power peaks:

\begin{equation*} \ D(\theta) = \frac{\lambda}{\Delta_1}(\theta-\theta_1)+\frac{1-\lambda}{\Delta_2} \end{equation*}

must be specified.

$\Delta_1$ and $\Delta_2$, 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 $cos^2(\theta-\theta_0);\theta\in [\theta_0 - \frac{\pi}{2};\theta_0+\frac{\pi}{2}]$
2. model in $exp\bigg(-0.5\Big(\frac{\theta-\theta_0}{s}\Big)^2\bigg);\theta\in [\theta_0 - \frac{\pi}{2};\theta_0+\frac{\pi}{2}]$
3. model in $cos^2\bigg(\frac{\theta-\theta_0}{2}\bigg);\theta\in [\theta_0 - \pi;\theta_0+\pi]$

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

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

$H_{m0}$: INITIAL SIGNIFICANT WAVE HEIGHT
$f_p$: INITIAL PEAK FREQUENCY
$\gamma$: INITIAL PEAK FACTOR
$\sigma_a$: INITIAL VALUE OF SIGMA-A FOR SPECTRUM
$\sigma_b$: INITIAL VALUE OF SIGMA-B FOR SPECTRUM
$\alpha_{phil}$: INITIAL PHILLIPS CONSTANT
$Fetch$: INITIAL MEAN FETCH VALUE
$f_{pmax}$: INITIAL MAXIMUM PEAK FREQUENCY
$\theta_1$: INITIAL MAIN DIRECTION 1
$s_1$: INITIAL DIRECTIONAL SPREAD 1
$\theta_2$: INITIAL MAIN DIRECTION 2
$s_2$: INITIAL DIRECTIONAL SPREAD 2
$\lambda$: 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.

8.5. Prescribing the boundary conditions

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

Only two kinds of boundary conditions are available in TOMAWAC.

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

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

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

8.5.1. The boundary conditions file

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

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

8.5.2. Prescribing the boundary conditions in the CAS file

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

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

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

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:

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

8.5.3. The LIMWAC user subroutine

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

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

8.6. Some useful subroutines


8.6.1. Modification of bottom topography: CORFON subroutine

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

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

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

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

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

8.6.2. Modifying the co-ordinates: CORRXY subroutine

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

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

8.6.3. Operations on vectors: OV subroutine

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

The syntax is as follows:

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

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

Example:

CALL OV('X=X+Y ', X, Y, Z, C, NPOIN)
Y is added to X, the result will be stored in X.

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