100

Modelling Evolution of Directional Spectra of Water Waves
S.Y. Annenkov and V.I. Shrira

Keele University, UK
Abstract

Studies of the evolution of water wave spectra are usually focussed on frequency orwavenumber spectra of energy or wave action. Much less is known about the shape, or evenwidth of the angular distribution and its evolution with time or fetch. Here, we consider
the evolution of angular width of wave spectra with and without wind forcing. As a metricof the angular width, we use the second moment of the distribution at certain wavenum‐

bers relative to the wavenumber of the spectral peak k_{p} (0.5, 1, 2 and 3k_{p} ), or averaged
over all wavenumbers. Numerical experiments are performed using three different numer‐
ical approaches employing different sets of assumptions: the classic Hasselmann kinetic

equation (KE), the generalized kinetic equation (gKE) and the direct numerical simulation(DNS) algorithm based on the Zakharov equation (DNS‐ZE). The results of simulations
by each model are compared with each other and with the Tehuantepec field observations

by Romero & Melville (2010). For the case of initially narrow swell the kinetic equations
considerably overestimate the rate of the initial angular broadening , compared with the
DNS‐ZE, although at large time the integrated angular width approaches the same pow‐
erlike asymptote. The angular width of the spectral peak continues to grow slowly evenfor large time. For constant wind forcing, the models agree well on the spectral width of
the spectral tail, but show large discrepancies for small wavenumbers. The DNS‐ZE simu‐
lations are more consistent with the data of the Tehuantepec observations, demonstrating
slow growth of angular width over time at the spectral peak, in contrast to nearly constantwidth shown by the kinetic equations. Results of the study demonstrate serious shortcom‐
ings of modelling based on the Hasselmann kinetic equation and provide an insight into therole of the assumptions underlying the Hasselmann equation.

1

Introduction

Studies of oceanic wind waves are mostly concerned with the evolution of wave spectra, mostly of

one dimensional (frequency or wavenumber) spectra of energy or wave action. Spatio‐temporal
evolution of directional characteristics is less known. Observations of the spatio‐temporal evo‐lution of directional spectra are very difficult and, therefore, rare; they require precise measure‐ments of two‐dimensional spectra with a good angular resolution over a sequence of fetches.
Parameterizations of wind wave spectra usually involve a representation of the directional dis‐tribution in the form D(\theta) , where \theta is the angle, while the dependence of the function D on

101
101
frequency or wavenumber is often ignored (Holthuisen, 2007). Even in this idealised case, theshape of the function D(\theta) is not well known, usually it is assumed it having a maximum inthe downwind direction (for which it is convenient to set \theta=0 ) and a gradual decrease withincreasing absolute value of \theta . Several expressions to describe this scale averaged behaviour ofD(\theta) have been proposed (e.g. \cos 2 model, \cos^{N} model, the wrapped‐normal distribution, etc).
On summarising a large number of observations (Young et al., 1995; Ewans, 1998) concluded

that the directional width also varies with frequency, increasing approximately two‐fold between\omega_{p} and 3\omega_{p} , where \omega_{p} is the frequency of the spectral peak. To take into account the frequency

dependence of D(\theta) Mitsuyasu et al. (1975) suggested a piecewise‐constant function, with thedirectional spreading being narrowest around the spectral peak frequency and widening towards
both lower and higher frequencies. The value of the function at the spectral peak was also related
to wind speed. There were further attempts at improvement the piecewise‐constant relationship
Goda (2010). Donelan et al (1985) suggested to use sech^{2}(\beta\theta) , where \beta is a piecewise‐constantfunction of frequency \omega/\omega_{p} , independent of wind speed. Little is known about the evolution ofD(\theta) with fetch or time.

Numerical simulations of long‐term spectral evolution also often focus on the evolution of one‐dimensional (frequency or wavenumber)spectra. Badulin & Zakharov (2017) performed extensive
simulations of swell evolution with the Hasselmann kinetic equation (KE). They have found thatregardless of the initial angular width there is a tendency to a universal directional distribution
at large times, which appeared to be quasi‐stationary at least at the frequency corresponding to
the spectral peak.
Until very recently, all numerical studies of long‐term evolution of random water wave fields
were based on the KE. The KE is based on two key assumptions: the statistical closure, occuring
as a result of an asymptotic procedure based on small nonlinearity for broadband wave fields
under the presumed absence of coherent patterns, and the quasistationarity implied by the large‐

time limit (Annenkov & Shrira,2006a ;

Newell & Rumpf, 2013). To describe situations with arapid transformation of a wave field, Annenkov & Shrira (2006a) derived a generalization of the
KE called the generalized kinetic equation (gKE) . The gKE employs the same statistical closure,but is free from the quasistationarity assumption. However, the role and validity of the statisticalclosure can be verified only by direct numerical simulation (DNS). Until recently, progress inclarifying this issue was very slow, since there were no tools for studying the long‐term evolutionof wave spectra without a statistical closure. Although there were a few attempts at modellingthe spectral evolution with DNS, starting with Tanaka (2001), the algorithms used were very
demanding in computational resources and, therefore, able to trace only the initial part of theevolution. The DNS algorithm, based on the Zakharov equation (DNS‐ZE), was first suggested
and implemented by Annenkov and Shrira (2001) for simulations of deterministic evolution ofwave field, then modified and extended for modelling random wave fields in Annenkov & Shrira(2009), it was later used for the study of the spectral evolution under gusty wind (Annenkov& Shrira, 2011) and for simulations of the evolution of higher statistical moments of wave field(Annenkov & Shrira, 2013; Shrira & Annenkov, 2013). The algorithm is based on integration of
the deterministic four‐wave Zakharov equation and ensemble averaging over a sufficient numberof realisations, and is free from any statistical assumptions.

Annenkov & Shrira (2018) examined the evolution of initially narrow (in both frequency and102
angle) spectra subjected only to dissipation localised at high frequencies, using all three numericalmodels: the KE (the standard WRT algorithm was kindly provided by G. van Vledder), thegKE and the DNS‐ZE. The reason for this particular choice of the initial conditions was the vast
amount of observational and numerical data relevant to these spectra. In particular, they wereused in the MARINTEK wave tank experiments described by Onorato et al. (2009) and Toffoli
et al. (2010), so that detailed observations of the initial stages of both the spectral evolution
and the evolution of higher statistical moments were available. Since these spectra were initially
far from equilibrium, at the initial stage of evolution they undergo relatively fast broadening.

Annenkov & Shrira (2018) found a striking difference in the rate of this angular broadening,which was much larger for the gKE and the KE than for the DNS‐ZE. At the same time, resultsof the DNS‐ZE simulations were consistent with earlier short‐term DNS simulations by Xiao et al.

(2013) carried out with the well‐established high‐order spectral method, the striking discrepancy
between the kinetic equations and DNS could have been reported earlier, but apparently was notnoticed. However, since the focus of Annenkov& Shrira (2018) was on the detailed comparison of
frequency spectra, all computations were performed within the same, relatively narrow, angularsector -4\pi/9\leq\theta\leq 4\pi/9 , which could have affected the long term evolution of angular spread.
In particular, the WRT code employed for the KE simulations proved to be sensitive to the
choice of the angular segment, especially for simulations of angular spreading.
In this paper, we consider the long‐term evolution of directional width of water wave spectrawith and without wind forcing. The aim is to compare evolution of directional spectra obtainedwith different models, in order to reveal the role of different sets of assumptions under different
winds. As a metric of directional width, we use the second moment of directional distribution,either at particular wavenumbers relative to the spectral peak k_{p} , or integrated over a wide rangeof wavenumbers.

First, we revisit the evolution of initially narrow spectra without wind, considered by An‐

nenkov & Shrira (2018). For the present study, the KE simulations of evolution have been rerun
for larger time and in the full circle rather than in a sector. At relatively short times, both KE
and gKE show close evolution of directional width. However, the DNS‐ZE simulations predict
much slower rate of directional spreading of the spectral tail. In the long term, all models appear
to tend to the powerlike spreading of approximately the same rate, although the DNS spectra
have considerably larger directional width of the spectral peak.
Second, we consider the evolution of directional width under a constant wind, with windspeed equal to 3c_{p} and 5c_{p} , where c_{p} is the phase speed of the initial spectral peak. The initialconditions, identical for all three models, are taken in the form of Donelan et al (1985) spectra forthe corresponding wind speed. Third, we perform numerical simulations, with all three numericalapproaches, of the evolution of wind waves observed in the Tehuantepec experiment (Romero &Melville, 2010) and compare the results with the observations.
The paper is organized as follows. Section 2 contains basic equations and brief description of
the numerical algorithms and parameters. In section 3, results of the numerical simulations are
discussed. Section 4 provides conclusions and a brief discussion.

103
2

Basic equations and numerics

2.1

The Zakharov equation and wave kinetic equations

Our starting point is the equation of motion in the form of the “four‐wave” integrodifferential
Zakharov equation originally derived for potential gravity waves at the surface of ideal incom‐pressible fluid of infinite depth, with the accuracy up to O(\varepsilon^{3}) (Zakharov ,1994)

i\frac{\partial b_{0} {\partial t} = \omega_{0}b_{0}+\int T_{0123}b_{1}b_{2}b_ {3}\delta_{0+1-2-3}dk_{123} ,196S ;

Krasitskii,(1)

where the wave field is expressed in terms of complex canonical variables b(k) linked to the Fourier
harmonics of the surface elevation \zeta(k) and velocity potential at the surface \psi(k) through an
integral‐power series:b( k)=\frac{1}{\sqrt{2} \{\sqrt{\frac{\omega(k)}{k} \zeta(k)+i\sqrt{\frac{k} {\omega(k)} \psi(k)\}+O(\varepsilon).

(2)Compact notation is used, designating arguments by indices, e.g. \delta_{0+1-2-3}=\delta(k_{0}+k_{1}-k_{2}k_{3}), dk_{123}=dk_{1}dk_{2}dk_{3} . The interaction coefficient T_{0123} is given by an explicit but lengthy

expression (e.g. Krasitskii, 1994).
Equation (1) is the basis for the derivation of the wave kinetic equation. Consider ensemblesof random wave fields, assuming spatial homogeneity. Then the spectral density of wave action
at wavevector k_{0} is the second‐order correlator n_{0},\{b_{0}^{*}b_{1}\}=n_{0}\delta_{0-1}.
By multiplying the Zakharov equation by b_{0}^{*} , upon ensemble averaging we get\frac{\partial n_{0} {\partial t}=2{\rm Im}\int T_{0123}\langle b_{0}^{*}b_{1} ^{*}b_{2}b_{3}\rangle\delta_{0+1-2-3}dk_{123} .

(3)

In a similar way we can express time derivative of the fourth‐order correlator in terms of thesixth‐order one and so on. This infinite chain of equations can be truncated invoking the closure

hypothesis (e.g. Zakharov et al., 1992; Nazarenko, 2011). Since the statistics of a linear wave
field is gaussian, the natural assumption is that the statistics of a weakly nonlinear field will bequasi‐gaussian. Following Zakharov et al. (1992), we assume that the third‐order correlator iszero. The fourth‐order correlator can be always presented as\{b_{0}^{*}b_{1}^{*}b_{2}b_{3}\rangle=[n_{0}n_{1}(\delta_{0-2}\delta_{1-3}+ \delta_{0-3}\delta_{1-2})]+J_{0123}^{(1)}\delta_{0+1-2-3},
where J_{0123}^{(1)} is the irreducible part of the correlator (cumulant), which determines the evolutionof n_{0} . To find the evolution of J_{0123}^{(1)} we express it in terms of the sixth‐order correlator, in which
we retain just the dominant terms—the products of pair correlators( i\frac{\partial}{\partial t}+\triangle\omega)J_{0123}^{(1)}=-2T_{0123} f_{0123} ,

(4)

104
\triangle\omega=\omega_{0}+\omega_{1}-\omega_{2}-\omega_{3}, f_{0123}=n_{2}n_{3}(n_{0}+n_{1})-n_{0}n_{1}(n_{2}+n_{3}) . The derivation of theclassical KE drops \partial/\partial t from the equation for the cumulant J_{0123}^{(1)} , which leads to the approximate
solution for large time in terms of generalised functions

where

J_{0123}^{(1)}(t)=-2T_{0123}[ \frac{P}{\triangle\omega}+ i\pi\delta(\triangle\omega)]f_{0123}(t)where P is “principal value\delta is Dirac 6‐function.

,

This leads to the KE\frac{\partial n_{0} {\partial t}=4\pi\int T_{0123}^{2}f_{0123}\delta_{0+1-2- 3}\delta(\omega_{0}+\omega_{1}-\omega_{2}-\omega_{3})dk_{123} ,(5)

The gKE is derived using the same statistical closure as the KE, but without the assumption ofquasi‐stationarity and the use of large time asymptotics. Instead, the equation for the cumulant

J_{0123}^{(1)}(4) is solved exactly (Annenkov & Shrira, 2006a )J_{0123}^{(1)}(t)=-2 iT_{0123}\int_{0}^{t}e^{-i\triangle\omega(\tau-t)}f_{0123} (\tau)d\tau+J_{0123}^{()}(0)e^{i\triangle\omega t}.
Then the resulting equation (gKE) has the form

\frac{\partial n_{0} {\partial t} = 4{\rm Re}\int\{T_{0123}^{2}[\int_{0}^{t}e^ {-i\triangle\omega(\tau-t)}f_{0123}(\tau)d\tau]‐\frac{i}{2}T_{0123}J_{0123}^{(1)}(0)e^{i\triangle\omega t}\}\delta_{0+1-2-3} dk_{123} .

(6)

The gKE is nonlocal in time: evolution of the spectrum depends on the previous history of
evolution, starting from the initial moment when the value of cumulant J_{0123}^{(1)}(0) is prescribed as
the initial condition. However, the gKE can be solved iteratively. On each time step, the valueof J_{0123}^{(i)} is computed as

J_{0123}^{(1)}(t)=-2 iT_{0123}\int_{0}^{t}e^{-i\triangle\omega(\tau-t)}f_{0123} (\tau)d\tau+J_{0123}^{(1)}(0)e^{i\triangle\omega t}(7)and taken as the new initial condition, so that the‘internal’ time integration is performed over one
timestep only. Analytical results concerned with the gKE are summarised in Shrira & Annenkov(2013).
2.2

Numerics

In this study, we use three different algorithms for three different models: the classic KE (Has‐selmann), the gKE , and our DNS algorithm based on the Zakharov equation (DNS‐ZE). For
the KE, we use the standard WRT code kindly provided by Gerbrant van Vledder. Two otheralgorithms are original.
Although the gKE at first glance looks much more complicated than the KE, and requires
to take into account a large number of interactions, including those far from resonance, the105
algorithm has a simpler structure and can be parallelized with nearly perfect scalability, thus
efficiently utilizing the advantages of modern supercomputers. The nonlocality in time can betreated by accumulating all the previous history in the value of the cumulant J_{0123}^{(1)} , which istaken as the initial condition at each time step. We employ the standard Runge‐Kutta‐Fehlberg
time‐stepping algorithm with automatic step choice. Details of the algorithm can be found in

Annenkov & Shrira (2016).
The DNS algorithm is based on the efficient algorithm for the simulation of the evolutionof discrete wave systems, employing the integrodifferential Zakharov equation (Annenkov andShrira, 2001). In that study, a wave field b(k, t) was considered as an ensemble of discrete
harmonics of the formb( k, t)=\sum_{j=1}^{N}b_{j}(k_{j}, t).

(8)The DNS‐ZE algorithm is essentially an application of the original algorithm of Annenkov andShrira (2001) to the evolution of continuous random wave fields. It is based on the idea ofcoarse‐graining of a continuous wave field, retaining its fundamental properties of nonlinear
interactions. The coarse‐graining represents a wave field by a grid consisting of wave packets,
coupled through exact and approximate resonant interactions. The resulting system of discrete
dynamical equations is again solved by Runge‐Kutta‐Fehlberg time‐stepping scheme. The DNS‐ZE algorithm was described in detail by Annenkov & Shrira ( 2013, 201S) .
For all computations in this work, we use non‐regular grids with N_{\omega} logarithmically spaced
points in the range \omega_{p}/2\leq\omega\leq 3\omega_{p} and N_{\theta} uniformly spaced angles. The KE is simulated in
the full circle of the angles \theta , since this is preferable for the existing WRT code. The gKE and
DNS‐ZE algorithms can work in any sector of \theta , and a sector -7\pi/9\leq\theta\leq 7\pi/9 , sufficiently
wide for the purposes of the present study, was used. The frequency resolution was chosen as
N_{\omega}=160 for the DNS‐ZE, and N_{\omega}=100 for both kinetic equations. Angular resolution was
N_{\theta}=70 within -7\pi/9\leq\theta\leq 7\pi/9 for the DNS, N_{\theta}=30 within the same sector for the gKE,and N_{\theta}=50 in full circle for the KE. In some computations which involved initially narrowspectra (§3.1 below), angular resolution for both kinetic equations was refined, and results wereverified to be non‐dependent on any further refinement.

3

Results

Here we consider results of simulations for several different cases of spectral evolution with orwithout wind forcing, focussing on the evolution of directional spreading. As a measure of
directional spreading, we use the second moment of the normalized angular distribution function

D(k, \theta) (Hwang et al., 2000, eg)\theta_{2}(k)=(\int_{0}^{\pi/2}\theta^{2}D(k, \theta)d\theta)^{1/2}(\int_{0}^{ \pi/2}D(k, \theta)d\theta)^{-1/2}(9)Often, directional spreading is characterized by the mean directional spread \theta_{m} , defined as theaverage of \theta_{2} overk

(Xiao et al., 2013; Annenkov & Shrira,201S ).

Being a scale averaged

106(D^{\Xi}smle3tau^{\cir}wedgbot{D
0

1000

2000

3000

4000

0

2000

4000

t (periods)

t

6000

8000

10000

(periods)

Figure 1: Evolution of averaged second moment of directional distribution \theta_{m} without wind
forcing. Initial conditions are JONSWAP spectra with \gamma=6 and narrow directional spreading,corresponding to (a) N=850(b)N=24 in the \cos^{N} directional model.

characteristics of the angular spectrum, \theta_{m} provides only the most basic information about its
width. Since D(k, \theta) , angular distributions at all wavenumbers play equal role, regardless ofthe spectral amplitude, diminishing the relative contribution from the spectral peak. In orderto understand the evolution of the angular distributions at key parts of the spectrum, at the
spectral peak, slope and front separately, along with \theta_{m} we consider the evolution of theta_{2} forat selected specific values of k , relative to the spectral peak k_{p} ( 0.5k_{p}, k_{p}, 2k_{p} and 3k_{p} ).
In the rest of this section, we consider three spectral evolution scenarios: evolution of initiallynarrow (both in frequency and angle) spectra without wind forcing, evolution of directionallywide spectra under constant wind, and simulation of a realistic observation‐based wind wave
evolution, with comparison with the experimental data.

3.1

Evolution of initially narrow spectra without wind

The laboratory experiment by Onorato et al. (2009), who studied short‐term (about 30 wave‐
length) evolution of wave spectra having initially JONSWAP form and different directional dis‐tributions provided inspiration for a few studies that aimed at modelling the same evolution

numerically. Xiao et al. (2013) performed numerical modelling using high‐order spectral methodwith high resolution, taking as initial conditions the spectra generated mechanically in the ex‐periment and simulating their evolution for about 150 wave periods. Two of these spectra had
JONSWAP form with peakedness \gamma=6 , significant wave height H_{s}=0.08m , the peak period
T_{p}=1s , and different directional distributions D(\theta) , given by cosine squareD(\theta)=\{ begin{ar y}{l \frac{2}\Theta}\cos^{2}(\frac{\pi\theta}{\Theta}) for|\theta|\leq\Theta/2 0 for|\theta|>\Theta/2 \end{ar y}

(10)

107

(D^{\Xi}smle3tau^{\cir}wedgbot{D10^{I}

10^{2}t (periods)

10^{3}

10’

10^{2}

10^{3}t (periods)

10^{4}

Figure 2: As in figure 1, in logarithmic coordinates, and power fits for the KE results for large
timewhere \Theta is the directional spreading width in radians. The first spectrum was very narrow inangle with \Theta=\pi/15 , approximately corresponding to N=840 in the \cos^{N} directional model.
The second spectrum was wider, with \Theta=\pi/2.9 , corresponding to N=24 . This angular width
can be considered as typical width of swell.

Recently, Annenkov & Shrira (2018) performed numerical experiments with the same initialconditions, using three numerical models (the KE, the gKE , and the DNS‐ZE), focussing mostly
on one‐dimensional spectral evolution. Using the advantages of the DNS‐ZE, these simulations
were performed for thousands of wave periods. The validation of the DNS‐ZE algorithm wasprovided by quantitative comparisons with the results of short‐term DNS simulations by Xiao

et al. (2013)) and tank observations by Onorato et al. (2009). Long‐term simulations allowed
Annenkov& Shrira (2018) to reveal fundamental differences between the predictions of the kinetic
equations and DNS. These differences were attributed to the role of the statistical closure.
One of these differences was in the rate of directional spreading at the initial stage of theevolution, which was found to be much faster for the kinetic equations than for the DNS‐ZE,

although the latter was in close agreement with the results by Xiao et al. (2013). The long‐term angular width evolution appeared to tend to the same limiting angle for all three models.

However, since the focus of Annenkov & Shrira (2018) was on one‐dimensional evolution, allcomputations were performed for a relatively narrow range of angles -4\pi/9\leq\theta\leq 4\pi/9 , identicalfor all models to facilitate the comparison between them.
For this work, all simulations involving the KE have been recomputed in the full circle, withthe base resolution of N_{\theta}=50, -\pi\leq\theta\leq\pi . For the narrow N=840 case, this angularresolution was found to be insuffient at the initial stage of the evolution, so N_{\theta} was increased
to 130. Simulations with the gKE and the DNS‐ZE were found to have sufficient resolution andto be unaffected by the angular limits for the time span of the simulations. In this work, thesesimulations are used unchanged.
Figure la, b shows the long‐term evolution of the averaged second moment \theta_{m} for both initial108\hat{CD^N}smile\trangf^{vj)}b
0

1000

2000
t (periods)

3000

4000

0

2000

4000
6000
t (periods)

8000

10000

Figure 3: Evolution of second moment of directional distribution \theta_{2}(k) at wavenumbers k=k_{p},0.5k_{p}, 2k_{p} and 3k_{p} , for two initially JONSWAP spectra with \gamma=6 and narrow directional

spreading, corresponding to (a) N=850(b)N=24 in the \cos^{N} directional model

spectra, obtained with the KE, the gKE and the DNS‐ZE. Computations for different models
were performed for different times, due to the numerical instability of the gKE in the highfrequencies, which is particularly effective without wind forcing (Annenkov & Shrira, 2018),
and high computational cost of the DNS‐ZE simulations, which involved averaging over 100

realizations. Case (b), with wider initial angular spectrum, was traced for larger time to get amore complete picture of long‐term behaviour. In the short term, both kinetic equations givenearly identical results and demonstrate much faster angular expansion than the DNS‐ZE. In
the long term, the least‐squares fit for the KE shows slow powerlike growth of \theta. , and all models
appear to tend to it (figure 2a,b ). The exponents for the powerlike growth are slightly different
for the two initial spectra, apparently due to a large difference in angular resolution. While \theta_{m}is an averaged characteristic, a more comprehensive picture of the directional spreading can beobtained from the evolution of the second moment \theta_{2}(k, t) at different k . This evolution, for a
few values of k relative to the spectral peak k_{p} , is shown in figure 3a,b . All three models nearlycoincide for \theta 2(0.5k_{p}, t) , but at the spectral peak and for higher k the DNS‐ZE demonstrates
a larger angular width at large times than the kinetic equations, although with close rates ofgrowth. However, the gradient of the angular width with respect to wavenumber is higher for
the KE, so that the difference in the averaged \theta_{m} becomes small at large times due to the roleof the highest wavenumbers.

Badulin & Zakharov (2017), having performed the long‐term simulations with the KE of theevolution of swell of different initial angular width, found that the angular width at the spectral
peak at large times tends to a certain universal value, although they provided only a qualitative
demonstration of the approach to this limit. Our simulation with the KE shows that the growthat the spectral peak continues even for large times ( \gg 10^{4} periods), although it becomes very
slow.

109

\hat{Pi^X}prod_\smle0{ut}do_3
0

500

1000

t

1500

(penods)

2000

2500

0

200

400

600

800

1000

t (periods)

Figure 4: Evolution of averaged second moment of directional distribution \theta_{m} for wind waves,with wund speed (a) U/c_{p}=3(b)U/c_{p}=5 , where c_{p} is the phase speed of the initial spectral
peak

3.2

Evolution of wind wave spectra with constant wind forcing

Here, we consider the evolution of wind wave spectra under the action of constant wind, for
two different wind speeds, corresponding to either moderate or strong wind. Initial conditions,
including the initial directional distribution, are specified in the form of the empirical spectra by

Donelan et al (1985) for U/c_{p}=3 (moderate wind) or 5 (strong wind), where U is the wind speed
and c_{p} is the initial phase speed of the spectral peak. The subsequent wind forcing, for the samewind speed, modelled according to Hsiao and Shemdin (1983), applied to frequencies \omega\leq 2.5\omega_{p},where \omega_{p} is the initial spectral peak frequency. Dissipation is applied to \omega>2.8\omega_{p} . Simulations
are performed with the same frequency discretizations as before. The angular discretization
involves 50 uniformly spaced angles for the KE, 30 angles within -7\pi/9\leq\theta\leq 7\pi/9 for the gKE
and 70 angles within the same sector for the DNS‐ZE. The DNS‐ZE algorithm includes a simple
parameterization of wave breaking: in all realizations all harmonics within the forcing domain
are subjected to strong dissipation instead of forcing if they exceed a certain level of nonlinearity.

Only a small number of harmonics is affected by this breaking parameterization (typically, 1−10
for moderate wind and 10‐20 for strong wind, out of the total number of 11431). Averaging forthe DNS‐ZE is over 20 realizations.

Figure 4 shows the evolution of the averaged second moment \theta_{m} for both cases. Since initial
spectra are already directionally wide, the change of angular width during evolution is muchsmaller than in the previous section. All spectra initially adapt to the constant wind forcing
over the first few hundred periods, after which the evolution of \theta_{m} is very similar for all models,
showing a gradual increase of angular width. For strong wind, there is a relatively small sys‐
tematic difference between the DNS‐ZE and both kinetic equations, which give almost identicalresults between themselves. This difference is likely to be due to a different shape of directional
distribution in the case of the DNS‐ZE, which requires a special study.

110\hat{Pi^N}smle\urcon30dt{}\mega2k_{p}

0 500 1000 1500 2000 2500 0 200 400 600 800 1000t (periods)t(periods)

Figure 5: Evolution of second moment of directional distribution \theta_{2}(k) at wavenumbers k=k_{p},0.5k_{p}, 2k_{p} and 3k_{p} for wind waves, with wind speed (a) U/c_{p}=3(b)U/c_{p}=5 , where c_{p} is thephase speed of the initial spectral peak

The details of the angular distribution, as shown in figure 5, reveal a close agreement forangular width at higher wavenumbers between the three models, but a large discrepancy forlarge scales, where the DNS‐ZE gives a considerably narrower angular distribution. For small
scales, the angular width continues to increase at large times, while at the spectral peak itsaturates.

3.3

Evolution of wind wave spectra in natural conditions

In the Tehuantepec experiment, Romero & Melville (2010) performed high quality airborne mea‐
surements of two‐dimensional wind wave spectra with simultaneous measurements of wind, overa few hundred kilometers of fetch. Here we model one particular example of this evolution (corre‐
sponding to research flight 10 of the experiment) numerically, taking as the initial condition thetwo‐dimensional spectrum measured at fetch 100 km offshore, where wind waves start dominat‐
ing over opposing swell. During the evolution, wind is gradually weakening and turning clockwise
by about 40 degrees, the wind data being input into the numerical models using, as before, theparamaterization by Hsiao and Shemdin (1983). The same numerical grid as in the previous
section was used, the initial condition and wind data renormalized accordingly. Averaging forthe DNS‐ZE was performed over 30 realizations. The opposing swell and its effect on wind waves
are not modelled.

Figure 6 shows the evolution of \theta_{m} for three numerical models, and comparison with \theta_{m}obtained from the data. It demonstrates a rough consistency of all the models with the obser‐
vations. Figure 7, which details the evolution of the angular width for different wavenumbers,
reveals a complex picture. For long waves, the angular width is identical for all the models, butquite different from the data, apparently due to the presence of swell. For the spectral peak andshorter waves, the DNS‐ZE shows good agreement with the data, demonstrating slow growth ofangular width over time at the spectral peak (in contrast to nearly constant width shown by the

111
111

\hat{Pi^X}B0J\smiledot{D}
0

ı000

2000

3000

4000

5000

6000

7000

t (periods)

Figure 6: Evolution of averaged second moment of directional distribution \theta_{m} for wind waves in

the Tehuantepec experiment (Romero & Melville, 2010), data and numerical modelling
kinetic equations) and weak dependence of angular width on wavenumber.

4

Concluding remarks

Better description of the evolution of the directional spectra is key for probing the adequacy ofour present understanding of the underpinning physical mechanisms. In this work, we examineda few representative examples of the evolution of directional wave spectra (with and without wind
forcing) obtained numerically using simulations with three different models, based on different
sets of assumptions. By comparing the results of simulations with different models and against
the Tehuantepec observations we have got important insights into the ability of existing numerical
tools to capture the angular evolution of wave spectra.
As a quantitative characteristics of wave field directionality we confined our attention to thedirectional width, as given by the second moment of the angular distribution function, calcu‐
lated at a few “ characteristic”’ wavenumbers relative to the spectral peak, or averaged over allwavenumbers. For the case of swell, initially narrow in angle, we have found that although the
kinetic equations, compared with the DNS‐ZE, considerably overestimate the rate of the initial
angular broadening, the integrated angular width at large time approaches the same powerlike
asymptote, corresponding to slow broadening with time. However, the DNS‐ZE shows a weakerdependence of the angular width on wavenumber. Within all models, the angular width of thespectral peak continues to grow slowly even for large time. For constant wind forcing, the mod‐els agree well on the spectral width of the spectral tail, but show large discrepancies for small
wavenumbers. For the case of the Tehuantepec experiment data, the DNS‐ZE is consistent withthe data, demonstrating slow growth of angular width over time at the spectral peak, in contrast

112\hat{Pi^nfy}smle\aprox^{0ust}O\dIJ
t (periods)

Figure 7: Evolution of second moment of directional distribution \theta_{2}(k) at wavenumbers k=k_{p},0.5k_{p}, 2k_{p} and 3k_{p} for wind waves in the Tehuantepec experiment, data and numerical modellingto nearly constant width shown by the kinetic equations, and weak dependence of angular width
on wavenumber. The overall conclusion is that the modelling based upon the Hasselmann kinetic
equation or gKE does not capture well the evolution of the angular width, DNS‐ZE does a muchbetter job in this respect.
The results presented have important, although not necessarily short term, implications.
The knowledge of the directional distribution of water waves is essential for various practicalapplications. Accurate measurements of two‐dimensional spectra are rare, and most of theexisting understanding of directional properties comes from modelling, which is almost exclusively
based on the Hasselmann kinetic equation. Understanding the role of the assumptions underlying
the KE is necessary for improving wave prediction models. From the theoretical perspective,
revealing the discrepancies between the kinetic equations and the DNS helps understand the roleof the statistical closure in kinetic models, which is important both within and beyond the water
wave context.

We are grateful to G. van Vledder for providing his WRT code, and to L. Romero and W.
K. Melville for the access to Tehuantepec experiment data. The work was supported by UK
NERC grant NE/I01229X/1 . Computations were performed on the CUDA High Performance
Computing cluster at Keele University and on the ECMWF Supercomputing facility within theSpecial Project SPGBSHRI.

References
Annenkov, S. Y. and Shrira, V. I. 2001 Numerical modelling of water waves evolution based on
the Zakharov equation. J. Fluid Mech. 449, 341‐371.

113
Annenkov, S. Y. & Shrira, V. I. 2006 Role of non‐resonant interactions in the evolution of
nonlinear random water wave fields. J. Fluid Mech. 561, 181‐207.

Annenkov, S. Y. & Shrira, V. I. 2009 “Fast” nonlinear evolution in wave turbulence. Phys. Rev.
Lett. 102, 024502.

Annenkov, S. Y. & Shrira, V. I. 2011 Evolution of wave turbulence under “gusty” forcing. Phys.
Rev. Lett. 107, 114502.

Annenkov, S. Y. & Shrira, V. I. 2013 Large‐time evolution of statistical moments of a wind wave
field. J. Fluid Mech. 726, 517‐546.
Annenkov, S. Y. & Shrira, V. I. 2016. Modelling transient sea states with the generalised kinetic
equation. In: Rogue and Shock Waves in Nonlinear Dispersive Media. Springer, pp. 159‐178.
Annenkov, S. Y. & Shrira, V. I. 2018 Spectral evolution of weakly nonlinear random waves:kinetic description versus direct numerical simulations. J. Fluid Mech. 844, 766‐795.
Badulin, S. I. & Zakharov, V. E. 2017 Ocean swell within the kinetic equation for water waves.
Nonlin. Proc. Geophys. 24, 237‐253.
Donelan, M. A., J. Hamilton & Hui, W. H. 1985 Directional spectra of wind‐generated waves.
Phil. Trans. R. Soc. London A315, 509‐562.
Ewans, K. C. 1998 Observations of the directional spectrum of fetch‐limited waves. J. Phys.
Oceanogr.28, 495‐512.
Goda, Y. 2010 Random seas and design of maritime structures. World Scientific.
Holthuisen, L. H. Waves in oceanic and coastal waters. Cambridge University Press.

Hwang, P. A., Wang, D. W., Walsh, E. J., Krabill, W. B. & Swift, R. N. 2000 Airborne measure‐ments of the wavenumber spectra of ocean surface waves. Part II: Directional distribution. J.
Phys. Oceanogr. 30, 2768‐2787.
Krasitskii, V. P. 1994 On reduced Hamiltonian equations in the nonlinear theory of water surface
waves. J. Fluid Mech. 272, 1‐20.

Mitsuyasu, H., Tasai, F., Suhara, T., Mizuno, S., Ohkusu, M., Honda, T. & Rikiishi, K. 1975
Observations of the directional spectrum of ocean WavesUsing a cloverleaf buoy. J. Phys.
Oceanogr. 5, 750‐760.
Nazarenko, S. V. 2011 Wave Turbulence. Springer.
Newell, A. C. & Rumpf, B. 2013 Wave turbulence: A story far from over. In: Advances in Wave
Turbulence, World Scientific Series on Nonlinear Science, 83, pp. 1‐51.

114
Onorato, M., Cavaleri, L., Fouques, S., Gramstad, O., Janssen, P.A.E.M., Monbaliu, J., Osborne,
A.R., Pakozdi, C., Serio, M., Stansberg, C.T. & Toffoli, A. 2009 Statistical properties ofmechanically generated surface gravity waves: a laboratory experiment in a three‐dimensional
wave basin. J. Fluid Mech. 627, 235‐257.

Hsiao, S. V. & Shemdin, O. H. 1983 Measurements of wind velocity and pressure with a wave
follower during MARSEN. J. Geophys. Res. 88, 9841‐9849.
Romero, L. & Melville, W. K. 2010 Airborne observations of fetch‐limited waves in the Gulf ofTehuantepec. J. Phys. Oceanogr. 40, 441‐465.
Shrira, V. I. & Annenkov, S. Y. 2013 Towards a new picture of wave turbulence. In: Advances
in Wave Turbulence, World Scientific Series on Nonlinear Science, 83, pp. 239‐281.
Tanaka, M. 2001 Verification of Hasselmann’s energy transfer among surface gravity waves bydirect numerical simulations of primitive equations. J. Fluid Mech. 444, 199‐221.

Toffoli, A., Gramstad, O., Trulsen, K., Monbaliu, J., Bitner‐Gregersen, E. & Onorato, M. 2010
Evolution of weakly nonlinear random directional waves: laboratory experiments and numer‐
ical simulations. J. Fluid Mech. 664, 313‐336.
Xiao, W., Liu, Y., Wu, G. & Yue, D. K. 2013 Rogue wave occurrence and dynamics by direct
simulations of nonlinear wave‐field evolution. J. Fluid Mech. 720, 357‐392.
Young, I. R., Verhagen, L. A. & Banner, M. L. 1995 A note on the bimodal directional spreadingof fetch‐limited wind waves. J. Geophys. Res. 100, 773‐778.
Zakharov, V. E.

196S

Stability of periodic waves of finite amplitude on the surface of a deep

fluid. J. Appl. Mech. Tech. Phys. (USSR) 9, 86‐94.
Zakharov, V. E., L’vov, V. S. & Falkovich, G. 1992 Kolmogorov Spectra of Turbulence I: Wave
Turbulence. Springer.