# Influence of varying bathymetry in rogue wave occurrence within unidirectional and directional sea-states

Guillaume Ducrozet, M. Gouin

► To cite this version:

Guillaume Ducrozet, M. Gouin. Influence of varying bathymetry in rogue wave occurrence within unidirectional and directional sea-states. Journal of Ocean Engineering and Marine Energy, 2017, 10.1007/s40722-017-0086-6. hal-01563227

HAL Id: hal-01563227

https://hal.science/hal-01563227v1

Submitted on 27 May 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
1

2

3

4

5

6

7

8

9

10

11

# Influence of varying bathymetry in rogue waves occurrence within unidirectional and directional sea-states

Guillaume Ducrozet · Maïté Gouin

Received: date / Accepted: date

**Abstract** Recent experimental and numerical studies have demonstrated an increased rogue waves' activity during the propagation of a wave field over a sloping bottom, from a deeper to a shallower domain. These studies have shown the influence of several parameters (wave steepness, amplitude of depth variation, slope profile, etc.) but were limited to unidirectional sea states. In this work, we focus on the effect of the directional spreading on the wave

G. Ducrozet and M. Gouin

Ecole Centrale Nantes, LHEEA Lab., UMR CNRS 6598

1, rue de la Noë - 44321 Nantes - France

Tel.: +33 240 371 645

Fax: +33 240 372 523

E-mail: guillaume.ducrozet@ec-nantes.fr

M. Gouin

Institut de Recherche Technologique Jules Verne

Chemin du Chaffault - 44340 Bouguenais - France
A highly nonlinear potential flow solver based on the High-Order Spectral method is used. We demonstrate that the enhancement of the extreme waves' occurrence observed close to the shallower side of the slope is reduced when considering the directionality of the sea state. We can state that the underlying physics is different between a real configuration and simplified unidirectional simulations. It is consequently essential to include directional effects in the context of rogue waves to have an accurate estimation of their probabilities of occurrence.

**Keywords** Rogue Waves · Varying bathymetry · Directional sea states · High-Order Spectral method

PACS 92.10.Sx · 92.10.Hm · 47.85.Dh · 47.35.Bb · 47.15.km · 47.11.Kb

## 1 Introduction

The understanding and description of extreme ocean waves have drawn a wide attention in the past years, see e.g. [Kharif et al., 2008, Osborne, 2010]. This is mainly due to their possible disastrous effects on structures at sea. Several physical phenomena are known to be possible candidates to explain such features [Kharif and Pelinovsky, 2003, Dysthe et al., 2008]. For instance, the non-linear wave interactions are gathering a wide community of scientists toward their in-depth understanding. The modulational instability of weakly nonlinear waves [Benjamin and Feir, 1967, Zakharov and Gelash, 2013] is indeed of major scientific interest due to its occurrence in a variety of different physical contexts [Onorato et al., 2013]. However, this is one among other physical
mechanisms and its pertinence to explain everything related to rogue waves is still an ongoing debate [Fedele et al., 2016]. It is well known that the formation of such waves may be attributed to focusing either due to dispersion (frequency focusing), to spatial location (directional focusing) or to the nonlinear wave-wave interactions. These three phenomena may be seen as competitive elements during the formation of extreme waves. For instance, for unidirectional waves and deep water, the linear dispersion may lead the evolution for broad spectra [Johannessen and Swan, 2003], but narrow-banded conditions are governed by the modulational instability [Onorato et al., 2013]. Nonlinear self-focusing also appears essential in the evolution of unidirectional wave fields in intermediate water depths [Katsardi and Swan, 2011]. The same applies when studying multidirectional sea states, which dynamics are substantially different for long and short-crested seas. The possible energy directional spreading induces a decrease of the nonlinear focusing due to modulational instability [Katsardi and Swan, 2006, Onorato et al., 2009, Toffoli et al., 2010, Fedele, 2015].

50

Recently, the propagation of gravity waves over variable bathymetry profiles has been studied as a possible configuration enhancing the occurrence of large waves. Different studies have described the statistical properties of gravity waves in this configuration both experimentally [Trulsen et al., 2012] and with different numerical methods ranging from Korteweg-de Vries models [Pelinovsky and Sergeeva, 2006, Sergeeva et al., 2011] through Non-Linear
Schrödinger (NLS) [Zeng and Trulsen, 2012] and Boussinesq models [Gramstad et al., 2013]
up to fully nonlinear potential flow solvers [Viotti and Dias, 2014].

57
58

In the present study, we take as a starting point the experimental work
of [Trulsen et al., 2012]. These laboratory experiments were conducted in a
wave basin at MARIN facility [Bunnik, 2010]. The propagation of long-crested
waves normally incident over a non-uniform bathymetry is studied. Different
irregular sea states are propagated over a 1 : 20 slope from a water of con-
stant depth 0.60 m to a shallower constant depth domain at 0.30 m. Specific
attention is paid to the spatial evolution of the statistical properties of the
wave field. The exact set-up is reminded in Sec. 4.1. Different behaviors have
been exhibited with respect to the propagation of irregular sea states in such
configuration. These physical features are dependent on depth variations, rel-
ative water depths in both domains, bottom slope, etc. However, it has been
demonstrated that there can be a local maximum of kurtosis and skewness
close to the shallower side of the slope, which is associated to a local maxi-
mum of probability of large wave envelope. This behavior can be interpreted
[Viotti and Dias, 2014] as a transition between two different statistical equi-
libria. As a result, non-equilibrium dynamics are induced in the wave field
which are responsible for the possible increase of rogue waves occurrence. Sim-
ilar behavior may be observed due to the initialisation process [Fedele, 2015],
change of wind direction [Annenkov and Shrira, 2009], wave-current interac-
tions [Toffoli et al., 2011], etc.

59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79

80

81

82

83

84

85

86

87

88

89

90

91

92

93

94

95

96

97

98

99

100

101

Nevertheless, those studies are limited to uni-directional sea states and it is expected that enhancing the complexity to realistic directional wave fields may change the observed results in 2D. As stated previously, the physical processes are usually different in directional sea-states and non-equilibrium dynamics may be affected, as it is the case for modulational instability [Onorato et al., 2009, Toffoli et al., 2010, Fedele, 2015].

The objective of this paper is consequently to assess the possible effects of directional spreading on the statistical properties of a wave field propagating from deeper to shallower depth. As a first step, we will consider in this study a bathymetric profile with infinite transversal length. Consequently, we do not consider **depth profiles resulting in focal regions, the so-called refractive focusing**, which may be of major importance for the wave statistics [Janssen and Herbers, 2009].

We choose to use numerical simulations for this study for their easiness in setting-up accurately controlled directional sea states as well as varying physical parameters, in comparison to experimental ocean wave basins. In terms of numerical methods, as noted in [Viotti and Dias, 2014], it appears important to consider a highly nonlinear numerical model in order to ensure the accurate simulation of those extreme waves, which by definition possibly include high degree of nonlinearity due to their large amplitudes.

For practical use, it is important to have an efficient and accurate numerical model. To this end, we will use a recently developed High-Order
Spectral (HOS) model adapted to varying bathymetry [Gouin et al., 2016, 102 Gouin et al., 2017]. It has been carefully validated in different configurations 103 and have shown its efficiency for the solution of the wave propagation of di- 104 rectional sea-states over varying bottom. 105

106

This paper is divided as follows. In a first part the numerical model used, 107 based on the HOS method [West et al., 1987, Dommermuth and Yue, 1987] is 108 briefly presented together with the chosen initial conditions and the important 109 numerical parameters. Then, the statistical moments used to characterize the 110 wave field evolution are presented, together with the method used for their 111 evaluation and some theoretical predictions. Next, the third part presents a 112 validation of the numerical model with unidirectional waves, comparing to 113 the experiments of [Trulsen et al., 2012]. Finally, the last section is dedicated 114 to the presentation of the results obtained with directional sea-states, with 115 comparisons to the unidirectional wave field. 116

117

## 2 Numerical method

2.1 High-Order Spectral method 118

### 2.1.1 Hypothesis and formulation of the problem 119

A 3D rectangular fluid domain and a Cartesian coordinate system with the 120 origin O located at one corner of the domain are considered, as presented in 121 Fig.1. $z = \eta(x, y, t)$ represents the free surface elevation and $z = 0$ the still 122
water level. The total depth $h(x, y)$ is decomposed into a mean depth $h_0$ and the bottom variation $\beta(x, y)$, such as $-h(x, y) = -h_0 + \beta(x, y)$, where $\beta$ is assumed to evolve in space but not in time. Thus, the considered domain becomes: $-h_0 + \beta(x, y) \le z \le \eta(x, y, t)$.

<figure><img src="image_2.png" /><figcaption>Fig. 1 Description of the computational domain with a mean water depth $h_0$ and a bottom variation $\beta(x, y)$.</figcaption></figure>

For this wave propagation problem, we use the potential flow theory. Given this hypothesis, the continuity equation becomes the Laplace equation in the fluid domain D. For the boundary conditions, we first consider periodic boundary conditions in the horizontal plane, assuming we have an infinite domain. The free-slip bottom boundary condition reads ($\tilde{\nabla}$ standing for the surface gradient):

$$ \tilde{\nabla}\phi \cdot \tilde{\nabla}\beta - \frac{\partial\phi}{\partial z} = 0 \quad \text{at } z = -h_0 + \beta(x, y) \qquad (1) $$

Following [Zakharov, 1968], both kinematic and dynamic nonlinear free-surface boundary conditions (FSBC) are written in terms of surface quantities $\eta$ and $\tilde{\phi}$ ($\tilde{\phi}(x, y, t) = \phi(x, y, z = \eta, t)$), and expressed at the exact free-surface position $z = \eta(x, y, t)$:

$$ \frac{\partial \eta}{\partial t} = \left( 1 + |\tilde{\nabla} \eta|^2 \right) \frac{\partial \phi}{\partial z} - \tilde{\nabla} \phi \cdot \tilde{\nabla} \eta \quad \text{on } z = \eta(x, y, t) \qquad (2) $$
$$ \frac{\partial \tilde{\phi}}{\partial t} = -g\eta - \frac{1}{2} |\tilde{\nabla} \tilde{\phi}|^2 + \frac{1}{2} \left( 1 + |\tilde{\nabla} \eta|^2 \right) \left( \frac{\partial \phi}{\partial z} \right)^2 \quad \text{on } z = \eta(x, y, t) \qquad (3) $$

To account for the time evolution of the quantities of interest $\eta$ and $\tilde{\phi}$ one needs to evaluate the vertical velocity at the free surface $W(x, y, t) = \frac{\partial \phi}{\partial z}(x, y, z = \eta(x, y, t), t)$ from the known quantities at the free-surface $\eta$ and $\tilde{\phi}$. This is the purpose of the HOS method described briefly in the following.

### 2.1.2 HOS method with variable bottom

The HOS model is a pseudo-spectral method with an arbitrary order of nonlinearity $M$, initially developed by [West et al., 1987] and [Dommermuth and Yue, 1987] for a flat bottom. We take as basis the open-source HOS-ocean solver [Ducrozet et al., 2016] available in such configuration. The first step to obtain the vertical velocity $W$ is to express the potential as a truncated power series of components $\phi^{(m)}$ for $m = 1$ to $M$: $\phi = \sum_{m=1}^{M} \phi^{(m)}$. $M$ stands for the order of nonlinearity of the method, which is arbitrary. The potential evaluated at the free surface is then expanded as a Taylor series with respect to the still water level $z = 0$.

When extended to variable bottom (we refer to [Gouin et al., 2016, Gouin et al., 2017] for details), the potential $\phi^{(m)}$ (solution of the problem) is expressed as the sum of two components: $\phi_{h_0}^{(m)}$ and $\phi_{\beta}^{(m)}$. $\phi_{h_0}$ satisfies the boundary condition on the flat bottom $z = -h_0$ and $\phi_{\beta}$ allows the definition of the correct bottom boundary condition (Eq.(1)) and satisfies a Dirichlet condition on $z = 0$. The latter condition allows to separate the free-surface problem from the varying bottom problem in the total solution. Both potentials are expressed in terms of basis functions taking into account previous boundary conditions and the
periodicity of the domain:

$$ \phi_{h_0}^{(m)}(x, y, z, t) = \sum_{p,q} A_{pq}^{(m)}(t) \frac{\cosh(k_{pq}(z+h_0))}{\cosh(k_{pq}h_0)} e^{ik_{x_p}x+ik_{y_q}y} \quad (4) $$

$$ \phi_{\beta}^{(m)}(x, y, z, t) = \sum_{p,q} B_{pq}^{(m)}(t) \frac{\sinh(k_{pq}z)}{\cosh(k_{pq}h_0)} e^{ik_{x_p}x+ik_{y_q}y} \quad (5) $$

with $k_{x_p} = p \frac{2\pi}{L_x}$, $k_{y_q} = q \frac{2\pi}{L_y}$ and $k_{pq} = \sqrt{k_{x_p}^2 + k_{y_q}^2}$. $A_{pq}^{(m)}(t)$ and $B_{pq}^{(m)}(t)$ are the so-called modal amplitudes of $\phi_{h_0}^{(m)}$ and $\phi_{\beta}^{(m)}$ respectively.

With this formulation and the use of the bottom boundary condition, one can compute the modal amplitudes $B_{pq}^{(m)}(t)$ as function of the $A_{pq}^{(m)}(t)$ with one of the two methods introduced in [Gouin et al., 2017]. With the first method, called original method, the same order of truncation $M$ is required at the free-surface ($\phi_{h_0}$) and at the bottom ($\phi_{\beta}$). For a possible enhanced accuracy and efficiency, the second method allows an increased flexibility with two truncation orders: $M$ for the free-surface, and $M_b$ for the bottom.

The last step of the HOS scheme is to evaluate the vertical velocity at the free surface $W$. In the original formulation [West et al., 1987, Dommermuth and Yue, 1987], a similar series expansion than the one done on the velocity potential is applied to the vertical velocity ($W = \sum_{m=1}^{M} W^{(m)}$), leading to another triangular system for $W^{(m)}$ solved iteratively. The same process is followed when studying varying bottom, with the previous decomposition of the velocity potential, solution of the problem, into two components. This way, once the modal amplitudes $B_{pq}^{(m)}(t)$ and $A_{pq}^{(m)}(t)$ have been computed with one of the two previous methods, the vertical velocity $W$ at the free surface can be obtained.
Then, the free surface conditions (Eqs. (2) and (3)) are used to march
in time $\eta$ and $\tilde{\phi}$ using the classical fourth-order Runge-Kutta scheme with
an adaptive time step [Ducrozet et al., 2016]. A time-tolerance parameter is
introduced in this concern. Note that a careful dealiasing is performed as
explained in [Bonnefoy et al., 2009] to ensure the method's convergence and
accuracy, even for waves with high steepness.

The initial HOS method on flat bottom is very efficient thanks to the use
of Fast Fourier Transforms (FFTs). The extension proposed still allows the
use of FFTs, preserving the numerical efficiency of the original HOS scheme
[Gouin et al., 2017]. Moreover, it has been demonstrated that the HOS method
with varying bottom conserves its convergence properties, i.e. exponential con-
vergence rate with respect to the number of discretization points $N_x$, $N_y$ and
the orders of nonlinearity $M$ and $M_b$ [Gouin et al., 2016, Gouin et al., 2017].
As a conclusion, the enhancement of the HOS-ocean model to take into account
variable bottom leads to an efficient and accurate model, which gives access
to the modeling of nonlinear directional wave fields with reasonable compu-
tational effort. Note that in the present study, the so-called original method
is used to take care of the varying depth in the HOS model. Consequently, a
unique order of nonlinearity $M$ is defined in the following, which stands for the
free surface nonlinearities as well as the ones induced by the bottom variation.
## 2.2 Initial conditions and relaxations zones

Any nonlinear wave model relies on the definition of appropriate initial conditions. In our case, we consider an initial linear irregular sea state defined from a wave spectrum, which may include directional spreading. To prevent from numerical instabilities, this initial condition is relaxed in time to allow the progressive introduction of the nonlinearities in the computation as well as the influence of the non-flat bathymetry [Dommermuth, 2000, Gouin et al., 2016].

In the following, the frequency spectrum $S_J(f)$ is defined as a JONSWAP spectrum specifying the peak period $T_p = 1/f_p$, the significant wave height $H_s$ and the shape factor $\gamma$:

$$S_J(\omega) = \frac{\alpha_J}{2\pi} H_s^2 f_p^4 f^{-5} \exp \left[ -\frac{5}{4} \left( \frac{f}{f_p} \right)^{-4} \right] \gamma^{\exp \left[ -\frac{(f-f_p)^2}{2\sigma^2 f_p^2} \right]} \quad (6)$$

with $\sigma = \begin{cases} 0.07 & \text{for } f < f_p \\ 0.09 & \text{for } f \ge f_p \end{cases}$ and $\alpha_J$ chosen to obtain the correct significant wave height.

Different irregular unidirectional wave fields will be tested, following the experimental work of [Trulsen et al., 2012]. Table 1 summarizes the parameters of these sea states. Note that the significant wave height $H_s$ has been adjusted to correspond exactly to the measurements available in [Trulsen et al., 2012]. The standard deviation measured in the deeper side of the domain is used for this minor modification. The value used in the following is consequently slightly different from the target one, which was $H_s = 0.06 m$ in the experiments.
<table><caption>**Table 1** Parameters of the wave-field tested</caption><thead><tr><th>Case</th><th>H<sub>s</sub> [m]</th><th>T<sub>p</sub> [s]</th><th>(k<sub>p</sub>h)<sub>deep</sub></th><th>(k<sub>p</sub>h)<sub>shallow</sub></th></tr></thead><tbody><tr><td>1</td><td>0.0624</td><td>1.273</td><td>1.6</td><td>0.99</td></tr><tr><td>2</td><td>0.0624</td><td>1.697</td><td>1.1</td><td>0.70</td></tr><tr><td>3</td><td>0.0624</td><td>2.121</td><td>0.81</td><td>0.54</td></tr></tbody></table>

For directional sea-states, the directional spectrum is defined as:

$$S(f, \theta) = S_J(f)D(\theta) \qquad (7)$$

The directional spreading at use is the wrapped-normal distribution suggested in [Mardia, 1972]:

$$D(\theta) = \frac{1}{2\pi} + \frac{1}{\pi} \sum_{j=1}^{N} \exp \left\{ -\frac{(j\sigma_{\theta})^2}{2} \cos (j[\theta - \theta_0]) \right\} \qquad (8)$$

with $\sigma_{\theta}$ the directional spreading standard deviation and $\theta_0$ the mean wave direction. N is the upper bound of the sum, which is a parameter fixed to N = 100 to ensure the convergence of the series. As a matter of completeness and possible comparisons to other studies, note that this directional distribution can be approximated as a more typical $cos^{2s}$ function [Mitsuyasu et al., 1975]

$$D_{2s}(\theta) = \frac{1}{\pi} 2^{2s-1} \frac{\Gamma^2(s+1)}{\Gamma(2s+1)} \cos^{2s} \left( \frac{\theta - \theta_0}{2} \right) \qquad (9)$$

where $\Gamma(\cdot)$ denotes the Gamma function and s is the spreading parameter. In the present study one narrow directional spreading, chosen as $\sigma_{\theta} = 11^{\circ}$ in Eq. (8) will be studied. This corresponds to a spreading parameter (Eq. (9)) s = 50. The directional functions are depicted in Fig. 2. This narrow distribution is expected to be representative of swells, especially when arriving in
limited water depth. Referring to Tab. 1, the relative water depth of the three cases studied is sufficiently small in the deeper side of the domain so that wave refraction have narrowed the directional spreading [Goda, 2000].

<figure><img src="image_3.png" /><figcaption>Fig. 2 Directional spreading function at use in the numerical tests.</figcaption></figure>

As a reminder, we have seen previously that the unknowns are expressed on basis functions that ensures periodic conditions in the horizontal plane (see Eqs. (4)-(5)). In the case of variable bathymetry, this is a possible issue if one wants to control the wave field that is propagating over a given bottom profile. Relaxation zones are used to overcome this problem. One generating zone is located at the beginning of the computational domain and one absorbing zone at the other end. The numerical solution is relaxed in space towards a specific solution (the linear one obtained from the wave spectrum defined above). The respective lengths of the relaxation zones is adjusted so that:

i) the generated wave field accurately satisfies the target wave spectrum ($H_s$ and $T_p$) and ii) no reflection occurs at the end of the domain. The reader can
refer to [Gouin et al., 2016] where the use of relaxation zones with the HOS method has been detailed.

## 3 Statistical moments

The statistical properties of interest are intended to characterize the extreme waves occurrence. We will consequently focus on the third and fourth-order moments of the free surface elevation, namely the skewness and the kurtosis defined as:

$$ \text{skewness} = \lambda_3 = \frac{\langle (\eta - \langle \eta \rangle^3) \rangle}{\sigma^3} \qquad (10) $$

$$ \text{kurtosis} = \lambda_4 = \frac{\langle (\eta - \langle \eta \rangle^4) \rangle}{\sigma^4} \qquad (11) $$

$\langle \cdot \rangle$ stands for the average over time and $\sigma$ is the standard deviation of $\eta$, directly related to the significant wave height ($H_s = 4\sigma$). The skewness characterizes the asymmetry of the distribution with respect to the mean while the kurtosis measures the importance of the tails. We recall that for a Gaussian distribution (i.e. linear waves), skewness = 0 and kurtosis = 3. Note also that those moments are evaluated locally in space from the temporal signals at every locations.

### 3.1 Numerical evaluation

These statistical properties are evaluated numerically by: i) averaging in time on a given window chosen as $T_w = 150 T_p$ once all wave components have propagated through the whole domain and ii) ensemble-averaged on different
random runs $n_{run} = 100$. For an accurate estimation of these quantities, it is important to assess the convergence of the ensemble-averaged data with respect to the number of runs. Figure 3 presents an example of evolution of the ensemble-averaged statistical moments for Case3. It is found that 100 random runs is sufficient to obtain a reliable estimate. The same applies to each case tested in the following.

<figure><img src="image_4.png" /><figcaption>Fig. 3 Evolution of the ensemble-averaged statistical moments evaluated at the beginning of the bar as a function of the number of runs (Case3).</figcaption></figure>
## 3.2 Theoretical predictions on flat bottom

A lot of efforts have been carried out in the past years to evaluate theoretically the skewness and kurtosis of a given sea states, see e.g. [Janssen and Onorato, 2007, Janssen, 2009, Janssen and Bidlot, 2009, Fedele, 2015, Fedele et al., 2016]. One of the main objective is to be able to assess if specific wave conditions are prone to extreme waves occurrence. Under different set of assumptions, it is possible to describe the effect of second and third-order nonlinearities on those statistical moments.

We are interested in theoretical predictions valid for intermediate to small water depths. To exhibit simple analytic expressions, we assume here that the sea state is homogeneous, weakly nonlinear and narrow-banded. We adopt the notations of [Fedele et al., 2016] and introduce different quantities that will be evaluated using the target initial wave spectrum (Eqs. (6)-(8)). The spectral moments are defined as $m_j = \int \omega^j S(\omega) d\omega$, allowing the definition of the mean angular frequency $\omega_m = m_1/m_0$ and the associated wavenumber $k_m$. The characteristic wave steepness is defined as $\mu_m = k_m \sigma$ and following quantities are also useful: $q_m = k_m h$, $Q_m = \tanh q_m$, $c_m = \omega_m / k_m$ the phase velocity, $c_s = \sqrt{gh}$ the shallow water phase velocity and $c_g = c_m [1 + 2q_m / \sinh(2q_m)] / 2$. Following [Serio et al., 2005, Janssen and Bidlot, 2009], the spectral bandwidth is evaluated as $\nu = 1 / (Q_p \sqrt{\pi})$ with $Q_p = 2 \int \omega S^2(\omega) d\omega / m_0^2$. The angular spreading is evaluated as $\sigma_\theta = \sqrt{2(1-a)/m_0}$ with $a = \iint \cos(\theta) S(\omega, \theta) d\omega d\theta$.
The first physical feature to take into account is the effect of bound modes. These are the dominant process for the skewness $\lambda_{3,NB}$ and one of the contribution to the kurtosis, denoted $\lambda_{4,NB}^b$

$$ \lambda_{3,NB} = \lambda_{3,NB}^b = 6\mu_m (\alpha + \Delta) \quad (12) $$

$$ \lambda_{4,NB}^b = 3 + \frac{4}{3}\lambda_{3,NB}^2 \left( 1 + \frac{\beta + \gamma}{2(\alpha + \Delta)^2} \right) \quad (13) $$

with the subscript NB standing for narrow-band approximation, the superscript b for the bound modes contribution and the introduction of the following parameters

$$ \alpha = \frac{3 - Q_m^2}{4Q_m^3}, \quad \Delta = -\frac{1}{4} \frac{c_s^2}{c_s^2 - c_g^2} \left[ 2 \frac{1 - Q_m^2}{Q_m} + \frac{1}{q_m} \right] \quad (14) $$

$$ \beta = \frac{24 + 3(1 - Q_m^2)^3}{64Q_m^6}, \quad \gamma = -\frac{\alpha^2}{2} \quad (15) $$

For the kurtosis, the dynamic contribution representing the effect of resonant interactions, may be important and as a consequence the complete kurtosis is written as the sum of two terms accounting for both the dynamic and the static contributions, namely $\lambda_{4,NB}^{dyn}$ and $\lambda_{4,NB}^b$ respectively:

$$ \lambda_{4,NB} = \lambda_{4,NB}^{dyn} + \lambda_{4,NB}^{b} \quad (16) $$

For the evaluation of this dynamic kurtosis in the context of deep-water directional waves, we introduce the Benjamin-Feir Index $BFI = \frac{\sqrt{2\mu_m}}{\nu}$ [Janssen, 2003] as well as the ratio between directional and frequency width $R = \frac{1}{2} \frac{\delta_g^2}{\nu^2}$. Then,
the dynamic kurtosis is expressed (for $0 \le R \le 1$) as a function of non-dimensional time $\tau = \nu^2\omega_0t$

$$ \lambda_{4,NB}^{dyn}(\tau, R) = BFI^2 J(\tau, R) \qquad (17) $$

$$ J(\tau, R) = 2Im \int_{0}^{\tau} \frac{1}{\sqrt{1 - 2i\alpha + 3\alpha^2\sqrt{1 + 2iR\alpha + 3R^2\alpha^2}}} d\alpha \quad (18) $$

This equation describes accurately the time-evolution of the kurtosis excess induced by resonant interactions occurring in a sea state [Fedele, 2015]. In the general case $0 < R < 1$ (focusing regime), this dynamic kurtosis reaches a maximum, which amplitude and time of occurrence is dependent on $R$ (and BFI): a smaller $R$ (i.e. smaller directional spreading) induces a larger maximum happening after longer time. The limiting case $R = 0$ describes the well-known unidirectional behavior where the kurtosis monotonically increases with time up to an asymptotic value. Note that for $R \ge 1$ (defocusing regime), it is possible to use the same equation with the following change of variable [Janssen and Bidlot, 2009]

$$ \lambda_{4,NB}^{dyn}(\tau, R) = -\frac{1}{R}\lambda_{4,NB}^{dyn}\left(\tau, \frac{1}{R}\right) \quad \text{for } R \ge 1 \qquad (19) $$

The evolution is then inverted: the dynamic kurtosis is negative with a minimum amplitude and location dependent on $R$.

Then, it is interesting to possibly extend the previous theoretical results to shallow water. It has been shown that the previous set of equations can be used with a modified Benjamin Feir Index valid for such conditions, namely
B_S [Janssen and Bidlot, 2009]. It is expressed as a function of the deep water

BFI:

$$B_S^2 = -BFI^2 \left(\frac{c_{g0}}{c_0}\right)^2 \frac{gX_{nl}}{k_0\omega_0\omega_0''} \quad (20)$$

with $c_0$ and $c_{g0}$ the phase and group velocities and (introducing $T_0 = \tanh(k_0h)$)

$$\omega_0'' = -\frac{g}{4\omega_0 k_0 T_0} \left\{ [T_0 - k_0 h (1 - T_0^2)]^2 + 4(k_0 h)^2 T_0^2 (1 - T_0^2) \right\} \quad (21)$$

$$X_{nl} = \frac{9T_0^4 - 10T_0^2 + 9}{8T_0^3} - \frac{1}{k_0h} \left[ \frac{(2c_{g0} - c_0/2)^2}{c_s^2 - c_g^2} + 1 \right] \quad (22)$$

Interestingly, we have now two important parameters: $B_S$ and $R$. For $k_0h < 1.363$, also referred as the defocusing regime, it is known that mod- ulational instability is canceled. This is reflected with a negative value for $B_S^2$ [Janssen and Onorato, 2007]. The dynamic kurtosis consequently becomes negative (with monotonous decrease in unidirectional wave fields). However, it is worth to note that this simple equation also reflects that for broad direc- tional spreading $R > 1$ and shallow water, it is expected to have an increase of extreme waves occurrence. This is consistent with some experimental mea- surements [Toffoli et al., 2009].

The previous expressions are valid in a limited framework, namely narrow- banded weakly nonlinear sea states evolving over constant depth. It is also important to recall that the wave spectrum is assumed to keep its shape during the evolution of the considered wave field. This is a strong hypoth-
349

350

351

352

353

354

355

356

357

358

esis, which has possibly large effect on the statistical moments of interest
[Mori et al., 2011]. In addition, we are presently interested in the evolution of
such quantities in an inhomogeneous medium (with varying depth). We will
consequently use those theoretical results outside of their application range,
assuming the behavior can be approximated to the one we have in two differ-
ent regions of constant depth. This is for sure a strong approximation. Note
that the change from time to space variables in Eq. (17) is achieved using the
group velocity.

However, we expect these theoretical predictions to provide possible useful
information on the processes at play in the configuration studied here.

359

4 Unidirectional wave-field

360

361

362

363

364

365

366

367

368

369

370

In this section, the previous HOS model is applied to the study of the evolu-
tion of the statistical properties of unidirectional gravity waves propagating
over a variable bathymetry. This follows different experimental and numer-
ical studies in which waves propagate from a deeper to a shallower domain
[Trulsen et al., 2012, Pelinovsky and Sergeeva, 2006, Sergeeva et al., 2011, Zeng and Trulsen, 2012,
Gramstad et al., 2013, Katsardi et al., 2013, Viotti and Dias, 2014]. We will take
as reference the experimental results of [Trulsen et al., 2012] that we intend to
reproduce numerically in order to validate our highly nonlinear HOS model in
this 2D configuration. These results are also useful in the view of assessing the
effect of directionality on the statistics of interest, which is treated in details
in Sec. 5.
## 4.1 Configuration

The geometry of the computational domain is defined to reproduce the experiments of [Trulsen et al., 2012]. We consequently refer to this publication for more details. In terms of geometry, the only difference with the wave tank configuration is associated to the presence of relaxation zones for the generation and absorption of waves as well as the necessary periodicity in the HOS model, see Sec. 2.2. As a consequence, the computational domain is defined as depicted in Fig. 4. The grey zones indicate the location of the relaxation zones, waves propagating from the left to the right. The stars stand for the probe locations in the experiments. We refer to [Trulsen et al., 2012] for the exact locations.

<figure><img src="image_5.png" /><figcaption>Fig. 4 Sketch of the computational domain</figcaption></figure>

In the experiments, the shallow side extends up to the absorbing zone. Note that the length of the shallow region is sufficient to ensure that it has no effect at the probe locations (where comparisons are to be made). Without
the relaxation zones, the computational domain is of length $L_x = 200$ m. In addition to the geometry of the computational domain, it is necessary to define carefully the numerical parameters for an accurate description of the evolution of the sea state of interest. A convergence study has been carried out, leading to a horizontal discretization $N_x = 1536$, an order of nonlinearity $M = 8$ and a time tolerance fixed to $T_{\text{ol}} = 10^{-6}$.

As an example of convergence study, Fig. 5 presents the influence of the nonlinear order $M$ on the free surface elevation at two locations. The first one is located in the deeper side, while the second one is located after the bottom slope, in the shallower domain. It is reminded that the so-called original method is used in this study for the treatment of bottom variation in the HOS model. Consequently, the same nonlinear order $M$ is used for the free surface and the bottom.

It is clear that for in the deeper side, with waves propagating over constant depth, the results obtained with $M = 5$ are converged. However, the behavior appears to be different after the bottom slope, in the shallower side, where small discrepancies are observed between the computations with $M = 5$ and $M = 8$ in some of the crests and troughs. As expected for typical wave propagation over a flat bottom, $M = 5$ is sufficient to capture all nonlinear free surface features. In the present configuration, one has to choose a higher value for $M$ in order to capture correctly the depth influence on the sea state evolution. Similar studies have been conducted with respect to the horizontal
<figure><img src="image_6.png" /><figcaption>**Fig. 5** Convergence study with respect to the HOS nonlinear order *M* (Case3). Temporal evolution of the free surface at two locations: first probe at deeper side (top) and last probe in the shallower side (bottom)</figcaption></figure>

discretization, leading to the choice of numerical parameters given previously.

## 4.2 Validation of numerical solutions

Figure 6 presents for Case3 the comparisons between the HOS simulations and the experimental data for the standard deviation, the skewness and the kurtosis along the tank. The theoretical predictions assuming flat bottom (Eqs. (12) & (16)) are also added for comparison.

The comparison between experiments and nonlinear simulations is very accurate in this configuration, validating the proposed numerical approach. In
<figure><img src="image_7.png" /><figcaption>Fig. 6 Standard deviation, skewness and kurtosis along the tank for Case3. Comparison between numerical simulations (solid line), experiments [Trulsen et al., 2012] (circles) and theoretical predictions (dashed line).</figcaption></figure>

more details, starting with the standard deviation, the comparison indicates that the incident wave field, before the step, has a similar amount of energy. As stated earlier, we adjusted the incident significant wave height so that it corresponds to the measurement (see Tab. 1). During the propagation on the slope, the behavior is slightly different, even if the increase observed experi-
mentally is captured. The last experimental point indicates a large decrease of the standard deviation which is a probable trace of physical dissipation, which is not taken into account in our potential flow solver.

Regarding the skewness and kurtosis, the increase and the local maxima after the end of the slope are correctly captured. An increased probability of occurrence of extreme waves is expected in this region. As indicated in [Viotti and Dias, 2014], there is a clear enhancement of the numerical predictions when using a highly nonlinear model, compared to simpler models such as the weakly nonlinear and weakly dispersive one of [Gramstad et al., 2013]. Note that the model is also able to capture correctly the negative excess kurtosis taking place in small relative water depths (below the theoretical threshold for modulational instability $kh < 1.363$) [Janssen and Onorato, 2007].

From the theoretical predictions, even if they suffer from strong assumptions (see Sec. 3.2) some interesting features are observed. The skewness in the deeper side of the tank is perfectly modeled while the effect of varying depth appears clearly from the evolution over the slope and in the shallower side. The theoretical model assumes constant depth at both side of the slope and the observed departure with experiments and numerical simulations is consequently a signature of the effects of the inhomogeneity in the medium. It is expected that on a longer spatial scale, the HOS numerical simulations will slowly converge towards this value.
Regarding the evolution of the kurtosis, the theoretical prediction is interesting to assess the effect of the bound and free waves components. Similarly to the skewness prediction (only due to bound waves) $\lambda_{4,NB}^b$ is positive and constant at both sides of slope with an increase in the shallower domain. The effect of resonant interactions, which is the difference between the two theoretical curves presented, appears to be very important. In the deeper side of the domain, no spatial evolution is observed since the length from the wave generation system is sufficiently long for the interactions to be settled. As expected, resonant interactions induces a negative contribution in this configuration, reminding that $k_p h = 0.81$. In the shallower side of the domain, the theory predicts a very large negative contribution of the free waves to the kurtosis. This takes place over a quite long distance, which appears consistent with the HOS nonlinear simulations. Surprisingly, these simple theoretical predictions are quite accurate in predicting the overall evolution of the kurtosis in the present configurations. Indeed, we remind that in this theory the wave spectrum is assumed of constant shape and narrow-banded, which is for sure not the case here.

The results for the evolution of the skewness and the kurtosis for the Cases 1 and 2 are reported in Fig. 7. These cases exhibit smaller rogue waves’ activity but are interesting to present as a matter of completeness with respect to the available experimental results [Trulsen et al., 2012].
<figure><img src="image_8.png" /><figcaption>**Fig. 7** Skewness and kurtosis along the tank for Case2 (left) and Case1 (right). Comparison between numerical simulations (solid line) and experiments [Trulsen et al., 2012] (circles) and theoretical predictions (dashed line).</figcaption></figure>

Conclusions are similar to Case3, even if for Case 1 some discrepancies appear, especially on the evaluation of the skewness, which is underestimated in the simulations. Note that a similar trend is observed in the computations of [Viotti and Dias, 2014]. The latter numerical results are difficult to analyze in more details because of the large oscillations they exhibit.

The theoretical predictions are also providing insightful information in those two configurations. At first, the asymptotic behavior of the skewness is accurately predicted in both cases. Regarding the kurtosis, the evolution in the deeper side presents some differences with previous case, associated to the relative water depth $k_p h = 1.1$ and $k_p h = 1.6$ for Case2 and Case1 respectively. The dynamic contribution to the kurtosis is consequently: i) close to zero for Case2 (we are close to the limit $k_p h = 1.363$ for modulational instability) and
ii) positive and dominant for Case1, which is in the focusing regime. After the slope, the spatial evolution is consistent with nonlinear simulations for Case2 and some discrepancies are observed for Case1, which exhibits limited contribution of the dynamic kurtosis. The latter is once again associated to the relative water depth $k_p h = 0.99$ in this configuration. It is expected that the changes in the wave spectrum may have significant influence, changing the *BFI* experienced at the beginning of the shallower side and consequently possibly explaining the behavior observed in the HOS nonlinear simulations and the experiments.

As a conclusion, we consider that the numerical method as well as the current set-up is validated for the proposed configuration. As demonstrated in [Gouin et al., 2017], the nonlinear model at use is very efficient, allowing a possible extension of the study to 3D directional sea states.

## 5 Multidirectional sea states

In the previous unidirectional framework, the enhanced probability of rogue waves is clear for waves propagating from deep to shallow water depths. During the wave propagation, the medium characteristics are changed, inducing non-equilibrium dynamics resulting in this feature. However, it is well established that the physics in the water wave propagation are essentially different in 3D. The extra degree of freedom offered in directional sea states (i.e. the possibility to spread in direction) usually reduces these non-equilibrium dynamics,
as seen for modulational instability [Onorato et al., 2009, Toffoli et al., 2010, Fedele, 2015]. This section intends to assess if this process is also at play when the out of equilibrium is induced by water depth variation.

## 5.1 Numerical set-up

The configuration is similar to the one presented for the unidirectional case (Sec. 4). The main difference is that we consider a directional sea state with a depth profile having an infinite transversal length. The waves are still generated and absorbed thanks to relaxation zones at both ends of the domain. We simply consider a transverse direction to Fig. 4. As an example, Fig. 8 presents for a specific case the free surface profile in the 3D computational domain, together with the bottom profile.

As detailed in Sec 2.2 we limit the present study to one medium directional spreading, chosen as $\sigma_{\theta} = 11^{\circ}$ in Eq. (8). The lateral size of the domain is chosen so that the wave spectrum is sufficiently discretized. This is consequently dependent on the chosen directional spreading. In the case considered, the transversal length is fixed to $L_y = 125$ m. Then, a convergence study allowed us to evaluate the necessary discretization in the transverse direction. Figure 9 presents an example of initial directional spectrum for the Case2. Summarizing, the numerical parameters chosen in the HOS simulations presented in this section are: $N_x \times N_y = 1536 \times 128$ discretization points (free of aliasing errors), a HOS nonlinear order $M = 8$ and a time tolerance fixed to $T_{ol} = 10^{-6}$.
<figure><img src="image_9.png" /><figcaption>**Fig. 8** Example of free surface elevation with directional sea states, Case 2 and medium spreading $\sigma_{\theta} = 11^{\circ}$ ($x-y$ ratio conserved).</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Fig. 9** Initial directional spectrum, Case2, medium directional spreading $\sigma_{\theta} = 11^{\circ}$. Color scale refers to the modal amplitude of the free surface elevation.</figcaption></figure>
Similarly to the unidirectional configurations, the statistical properties of
these directional sea state are obtained with ensemble-averaged quantities on
100 random runs. The average is achieved over the last 150 $T_p$ of the simulation.

## 5.2 Results and discussion

It has been demonstrated in [Ducrozet et al., 2017] that over a flat bottom, the limits of applicability of the HOS-ocean solver are associated to the presence of breaking waves in the simulations. These events are presently not approximated nor filtered out in order to perfectly control the accuracy of the numerical solution. In this concern, it has been shown that the unidirectional and directional configurations may have very different behaviors, especially for small relative water depths. This was a confirmation of the possible importance of the directional effects, when studying ocean wave properties.

As a consequence, it is observed in the present configuration that the Case3, which exhibits for unidirectional seas the highest rogue waves activity, is not accessible to the conservative HOS-ocean solver. A possible explanation resides in the fact that for small relative water depths, the phase velocities of the frequency components with most of the energy are similar. The directional spreading allows the superposition of waves traveling with different directions and the inherent possible wave breaking. This effect is completely absent from a unidirectional simulation where all those frequency components propagate with similar velocities and low probability of superposition. Some methods may be used to take care of wave breaking in potential flow solvers such as
HOS [Perlin et al., 2013, Seiffert and Ducrozet, 2016] but are not used in the present study. As a consequence, the following results are restricted to Case1 & Case2, which do not exhibit any problem of stability during the wave propagation.

Figure 10 presents, for Case1, the results for the standard deviation, skewness and kurtosis of the directional simulations together with the previous 2D results and unidirectional experiments. The theoretical predictions of skewness and kurtosis for this directional sea-state are also included.

As seen previously with the unidirectional simulations and experiments, the change of water depths in this case induces limited effects on both third and fourth moments. This is associated to the range of relative water depth involved in this configuration. As studied in details in [Gramstad et al., 2013, Viotti and Dias, 2014], for a given configuration, the decrease of the shallower depth will result in an increase of the deviations induced by nonlinearities. As a general comment, the simulations performed with a directional sea state also exhibit small influence of the change of bathymetry on the statistical properties of the wave field for Case1.

However, after noticing the exact correspondence of the configurations with respect to the energy content, some differences can be pointed out. At first, slight differences may be seen on the skewness, especially after the propagation over the slope. It is observed that the skewness experiences an increase when considering the directional spreading in the wave field. The effect of direction-
<figure><img src="image_11.png" /><figcaption>Fig. 10 Standard deviation, skewness and kurtosis along the tank for Case1. Comparisons between 2D experiments, unidirectional (HOS-2D) and directional (HOS-3D) computations and theoretical predictions.</figcaption></figure>

ality on the skewness is actually dependent on the relative water depth. For in-
termediate or deep water depths, the directional spreading induces a reduction
of the skewness [Onorato et al., 2009]. In smaller water depths ($k_p h < 1.36$) it
has been observed in experiments and numerical simulations that short-crested
seas exhibits larger values of skewness than long-crested ones [Forristall, 2000,
Toffoli et al., 2006, Toffoli et al., 2009]. The latter feature is actually not taken into account in the theoretical predictions for skewness (noticing that Eq. (12) is independent of directional spreading). The apparent better agreement with the 3D computations is consequently fortuitous. It is well known that during the propagation over a flat bottom, the deviation of skewness is mainly dominated by the bound modes, even though the free waves' dynamics can also contribute, see e.g. [Onorato et al., 2005]. The latter may explain the observed differences, keeping in mind that the deviations are weak in this configuration.

For the kurtosis, the differences between unidirectional and spreading waves are more obvious. The directionality appears especially sensible in the deeper flat bottom region. The possible kurtosis enhancement over the flat region due to nonlinearities is noticeably reduced in the case of directional wave field. The present results are in concordance with existing literature (see [Onorato et al., 2009, Toffoli et al., 2010]) and theoretical predictions. The kurtosis excess is now dominated by bound modes (the two theoretical curves are perfectly superposed), contrary to the 2D configuration presented in Fig. 7.

Then, the propagation in the varying water depth region does not induce any significant change on the tail of the probability distribution of the free surface elevation in both 2D and 3D configurations. Extreme waves activity is consequently reduced in this set-up. The discrepancies observed with the theoretical predictions in the shallower side of the domain is attributed to the change in the wave spectrum during its propagation. This leads to a change of the BFI as well as the directional to frequency width ratio *R* in Eq. (17).
This has been shown to possibly have large influence in the prediction of the kurtosis evolution [Mori et al., 2011].

We now focus on Case2, which exhibits more pronounced effects of the varying bottom on the wave statistics. Figure 11 presents the evolution of the different wave statistics along the computational domain.

Regarding the skewness, we observe that the effect of directionality is very limited in this configuration. The asymptotic value in deeper and shallower side of the slope is close in the two presented configurations and is the result of bound harmonics in the wave field. This is consistent with previous experimental observations [Onorato et al., 2009] and with the present choice of directional spreading $s = 50$ (equivalent to $N = 100$ in [Onorato et al., 2009]). When propagating in varying water depth, the directional wave field exhibits a slightly lower value for the skewness (and especially its maximal value).

In addition, the directional spreading has a strong influence on the kurtosis distribution throughout the whole domain. Similarly to Case1 and to existing literature, the kurtosis is reduced in the flat region when adding directionality. Some discrepancies are observed with the theoretical evolution of $\lambda_{4,NB}$ in this region, which is dominated by bound waves when directionality is added. As previously, this may be attributed to the spatial changes in the frequency content of the wave spectrum, not taken into account. The evolution of the kurtosis in the nonlinear simulations suggests that the dynamic contribution plays an active role in this configuration.
<figure><img src="image_12.png" /><figcaption>**Fig. 11** Standard deviation, skewness and kurtosis along the tank for Case2. Comparisons between 2D experiments and unidirectional (HOS-2D) and directional (HOS-3D) computations</figcaption></figure>

The major finding here is that the non-equilibrium dynamics induced by the depth change is greatly influenced by the directional spreading of the sea state considered. The enhancement of the extreme waves' activity is reduced in a real configuration compared to unidirectional simulations or experiments.
As a consequence, it is of major importance, in the context of rogue waves, to
include directional effects to have an accurate estimation of their probabilities of occurrence. It is also worth to indicate that with 3D simulations, negative kurtosis is also observed in the shallow side of the domain, as it is in 2D. Finally, kurtosis appears to tend toward different asymptotic values after the step, with a larger value in directional sea states than unidirectional ones.

As a final comparison, Fig. 12 presents different probability density func-
tions (pdf) of the surface elevation. The unidirectional and directional results
are compared and different locations are studied: i) before the slope in the
deeper region ($x = x_{bar} - 5.4m$), ii) just after the slope where statistical mo-
ments are maximum ($x = x_{bar} + 1.1m$) and iii) further away in the shallower
side ($x = x_{bar} + 6.6m$). We focus on Case2, which exhibits the larger deviations
in terms of statistical moments.

This figure illustrates the previous results on statistical moments with the
probability density functions. The space evolution of the pdf is clear in both
configurations with the increased asymmetry of the distribution just after the
slope (characterized through skewness). The importance of large events is more
complicated to accurately assess with those pdf, kurtosis being an interesting
parameter to have a clearer view about this point. However, it is still possible
to demonstrate the influence of directional spreading on the probabilities of
occurrence of extreme waves with the pdf, leading to the same conclusions
than previously.
<figure><img src="image_13.png" /><figcaption>Fig. 12 Probability density functions of surface elevation at different locations with unidirectional and directional sea states, Case2.</figcaption></figure>

## 6 Conclusion

This paper studies the propagation of wave fields over a variable bathymetry consisting in a deeper and a shallower regions connected thanks to a sloping bottom, following the experiments presented in [Trulsen et al., 2012]. The main objective is to assess the effect of directional spreading on the wave statistics in such configuration. Previous experimental and numerical stud-
ies were limited to unidirectional sea-states with varying parameters (rel-
ative water depths of the two regions, bottom slope and profile, etc.), see
[Trulsen et al., 2012, Pelinovsky and Sergeeva, 2006, Sergeeva et al., 2011, Zeng and Trulsen, 2012,
Gramstad et al., 2013, Viotti and Dias, 2014].

To this end, we focus on the simple configuration of directional sea states
encountering a submerged step of infinite transversal extent. We do not con-
sider depth profiles resulting in focal regions such as bathymetric lenses that
can induce a local focusing of the sea state behind it. Specific features have
been observed in this context, such as the uncoupling between the local wave
heights and the kurtosis [Janssen and Herbers, 2009]. This configuration has
already been widely studied experimentally (see e.g. [Berkhoff et al., 1982,
Vincent and Briggs, 1989, Bobinski et al., 2015]) and numerically (see [Belibassakis et al., 2001,
Griffiths and Porter, 2012, Gouin et al., 2017] among others).

At first, the numerical model that is based on the extension of the HOS-
ocean solver [Ducrozet et al., 2016] to a variable bottom [Gouin et al., 2016,
Gouin et al., 2017] is presented and validated against existing experimental
results with unidirectional wave fields. The accuracy of the model is demon-
strated and thanks to its efficiency, the simulation of directional wave fields is
accessible in reasonable computational time.

Then, the effective influence of the directionality in a wave field is inves-
tigated through the use of a narrow spreading. The comparisons between 2D
and 3D simulations are achieved with particular attention to the statistical
moments of the free surface elevation, namely the skewness and the kurtosis.
The latter information is intended to characterize accurately the deviations
from Gaussian properties in the considered sea-states. We highlight that for
the skewness, which is mainly dominated by the bound modes, the effect of
directionality is limited and dependent on the tested configuration.

As a major point, we demonstrate in this study that for the kurtosis the
non-equilibrium dynamics induced by the depth change is greatly influenced
by the directional spreading of the sea state considered. The resulting enhance-
ment of the extreme waves' activity is reduced in a real configuration compared
to simplified unidirectional simulations. Thus, it is essential to include direc-
tional effects in the context of rogue waves to have an accurate estimation of
their probabilities of occurrence.

As a consequence, this work includes the depth variation as another feature
inducing out-of-equilibrium dynamics, similar to unstable initial conditions,
change of wind or current conditions, etc. These dynamics are expected to be
in all configurations an essentially 3D phenomena in which the unidirectional
assumption leads to fundamental different physics. We hope that this work
will motivate the community to pursue the research for the effects of depth
variation on wave propagation in this direction. For instance, more system-
atic studies can be conducted in order to study in greater details the effects of
the different parameters of interest: i) respective values of kh<sub>deep</sub> and kh<sub>shallow</sub>
which conduct for unidirectional wave fields to different behaviors of the statis-
tical moments [Trulsen et al., 2012], ii) directional spreading, iii) peakedness
of the spectrum, iv) level of nonlinearities, etc. The study of the long distance evolution of the sea states after the slope is also of great interest, together with a more detailed analysis on the respective contributions of bound and free waves.

In addition, this study demonstrates the ability of this highly nonlinear model to simulate this complex problem in an efficient way, opening the path to numerous interesting studies to characterize, for instance, the extreme waves' occurrence. As a final point, the main limitation of the model is associated to the breaking events which cannot be modeled presently. Inclusion of a wave breaking model to be able to treat the configurations exhibiting the larger occurrences of extreme waves is also a direction where research efforts are currently carried on [Seiffert and Ducrozet, 2016].

**Acknowledgements** The authors would like to thank the two anonymous referees for their constructive comments.

## References

[Annenkov and Shrira, 2009] Annenkov, S. Y. and Shrira, V. (2009). Evolution of kurtosis for wind waves. *Geophysical Research Letters*, 36(13).

[Belibassakis et al., 2001] Belibassakis, K., Athanassoulis, G., and Gerostathis, T. P. (2001). A coupled-mode model for the refraction-diffraction of linear waves over steep three-dimensional bathymetry. *Applied Ocean Research*, 23(6):319-336.

[Benjamin and Feir, 1967] Benjamin, T. B. and Feir, J. (1967). The disintegration of wave trains on deep water. Part 1. Theory. *Journal of Fluid Mechanics*, 27(03):417-430.
[Berkhoff et al., 1982] Berkhoff, J., Booy, N., and Radder, A. C. (1982). Verification of numerical wave propagation models for simple harmonic linear water waves. *Coastal Engineering*, 6(3):255-279.

[Bobinski et al., 2015] Bobinski, T., Eddi, A., Petitjeans, P., Maurel, A., and Pagneux, V. (2015). Experimental demonstration of epsilon-near-zero water waves focusing. *Applied Physics Letters*, 107(1):014101.

[Bonnefoy et al., 2009] Bonnefoy, F., Ducrozet, G., Le Touzé, D., and Ferrant, P. (2009). Time-domain simulation of nonlinear water waves using spectral methods. *Advances in Numerical Simulation of Nonlinear Water Waves*, 11:129-164.

[Bunnik, 2010] Bunnik, T. (2010). Benchmark workshop on numerical wave modeling - description of test cases. Technical Report 70022-1-RD, MARIN, The Netherlands.

[Dommermuth, 2000] Dommermuth, D. (2000). The initialization of nonlinear waves using an adjustment scheme. *Wave Motion*, 32(4):307-317.

[Dommermuth and Yue, 1987] Dommermuth, D. and Yue, D. (1987). A high-order spectral method for the study of nonlinear gravity waves. *Journal of Fluid Mechanics*, 184:267-288.

[Ducrozet et al., 2016] Ducrozet, G., Bonnefoy, F., Le Touzé, D., and Ferrant, P. (2016). HOS-Ocean: Open-source solver for nonlinear waves in open ocean based on High-Order Spectral method. *Computer Physics Communications*.

[Ducrozet et al., 2017] Ducrozet, G., Bonnefoy, F., and Perignon, Y. (2017). Applicability and limitations of highly non-linear potential flow solvers in the context of water waves. *Ocean Engineering*, submitted.

[Dysthe et al., 2008] Dysthe, K., Krogstad, H. E., and Müller, P. (2008). Oceanic rogue waves. *Annu. Rev. Fluid Mech.*, 40:287-310.

[Fedele, 2015] Fedele, F. (2015). On the kurtosis of deep-water gravity waves. *Journal of Fluid Mechanics*, 782:25-36.

[Fedele et al., 2016] Fedele, F., Brennan, J., De León, S. P., Dudley, J., and Dias, F. (2016). Real world ocean rogue waves explained without the modulational instability. *Scientific reports*, 6.
747

Forristall, 2000 [Forristall, G. Z. (2000). Wave crest distributions: Observations and second-order theory. *Journal of physical oceanography*, 30(8):1931–1943.

748

749

[Goda, 2000] Goda, Y. (2000). *Random seas and design of maritime structures*. Advanced Series on Ocean Engineering: Vol. 15. World Scientific.

750

751

[Gouin et al., 2016] Gouin, M., Ducrozet, G., and Ferrant, P. (2016). Development and validation of a non-linear spectral model for water waves over variable depth. *European Journal of Mechanics-B/Fluids*, 57:115–128.

752

753

754

[Gouin et al., 2017] Gouin, M., Ducrozet, G., and Ferrant, P. (2017). Propagation of 3d nonlinear waves over an elliptical mound with a high-order spectral method. *European Journal of Mechanics-B/Fluids*, 63:9–24.

755

756

757

[Gramstad et al., 2013] Gramstad, O., Zeng, H., Trulsen, K., and Pedersen, G. (2013). Freak waves in weakly nonlinear unidirectional wave trains over a sloping bottom in shallow water. *Physics of Fluids*, 25(12):122103.

758

759

760

[Griffiths and Porter, 2012] Griffiths, L. and Porter, R. (2012). Focusing of surface waves by variable bathymetry. *Applied Ocean Research*, 34:150–163.

761

762

[Janssen and Bidlot, 2009] Janssen, P. and Bidlot, J.-R. (2009). On the extension of the freak wave warning system and its verification. Technical Memorandum 588, ECMWF.

763

764

[Janssen, 2003] Janssen, P. A. (2003). Nonlinear four-wave interactions and freak waves. *Journal of Physical Oceanography*, 33(4):863–884.

765

766

[Janssen, 2009] Janssen, P. A. (2009). On some consequences of the canonical transformation in the hamiltonian theory of water waves. *Journal of Fluid Mechanics*, 637:1–44.

767

768

[Janssen and Onorato, 2007] Janssen, P. A. and Onorato, M. (2007). The intermediate water depth limit of the zakharov equation and consequences for wave prediction. *Journal of Physical Oceanography*, 37(10):2389–2400.

769

770

771

[Janssen and Herbers, 2009] Janssen, T. T. and Herbers, T. (2009). Nonlinear wave statistics in a focal zone. *Journal of Physical Oceanography*, 39(8):1948–1964.

772

773

[Johannessen and Swan, 2003] Johannessen, T. B. and Swan, C. (2003). On the nonlinear dynamics of wave groups produced by the focusing of surface-water waves. *Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences*, 459(2032):1021–1052.

774

775

776
777

778

779

[Katsardi et al., 2013] Katsardi, V., De Lutio, L., and Swan, C. (2013). An experimental study of large waves in intermediate and shallow water depths. part i: wave height and crest height statistics. *Coastal Engineering*, 73:43-57.

780

781

782

783

[Katsardi and Swan, 2006] Katsardi, V. and Swan, C. (2006). The importance of directionality in the evolution of large waves in intermediate and shallow water. In *Proc. of the 30th Intl. Conf. Coastal Engng.*, San Diego, USA. ASCE AMERICAN SOCIETY OF CIVIL ENGINEERS.

784

785

786

787

[Katsardi and Swan, 2011] Katsardi, V. and Swan, C. (2011). The evolution of large non-breaking waves in intermediate and shallow water. i. numerical calculations of uni-directional seas. *Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences*, 467(2127):778-805.

788

789

[Kharif and Pelinovsky, 2003] Kharif, C. and Pelinovsky, E. (2003). Physical mechanisms of the rogue wave phenomenon. *European Journal of Mechanics-B/Fluids*, 22(6):603-634.

790

791

[Kharif et al., 2008] Kharif, C., Pelinovsky, E., and Slunyaev, A. (2008). *Rogue Waves in the Ocean*. Springer Science & Business Media.

792

793

[Mardia, 1972] Mardia, K. (1972). Probability and mathematical statistics: statistics of directional data. *London: Academic*.

794

795

796

[Mitsuyasu et al., 1975] Mitsuyasu, H., Tasai, F., Suhara, T., Mizuno, S., Ohkusu, M., Honda, T., and Rikiishi, K. (1975). Observations of the directional spectrum of ocean waves using a cloverleaf buoy. *Journal of Physical Oceanography*, 5(4):750-760.

797

798

799

[Mori et al., 2011] Mori, N., Onorato, M., and Janssen, P. (2011). On the estimation of kurtosis in directional sea states for freak waves forecasting. *Journal of Physical Oceanography*, 41:1484-1496.

800

801

802

803

[Onorato et al., 2009] Onorato, M., Cavaleri, L., Fouques, S., Gramstad, O., Janssen, P., Monbaliu, J., Osborne, A., Pakozdi, C., Serio, M., Stansberg, C., et al. (2009). Statistical properties of mechanically generated surface gravity waves: a laboratory experiment in a three-dimensional wave basin. *Journal of Fluid Mechanics*, 627:235-257.

804

805

[Onorato et al., 2005] Onorato, M., Osborne, A., and Serio, M. (2005). On deviations from gaussian statistics for surface gravity waves. In *Rogue Waves: Proceedings of the Hawaiian*
Winter Workshop.

[Onorato et al., 2013] Onorato, M., Residori, S., Bortolozzo, U., Montina, A., and Arecchi, F. (2013). Rogue waves and their generating mechanisms in different physical contexts. *Physics Reports*, 528(2):47–89.

[Osborne, 2010] Osborne, A. (2010). *Nonlinear ocean waves and the Inverse Scattering Transform*, volume 97 of *International Geophysics Series*. Academic Press.

[Pelinovsky and Sergeeva, 2006] Pelinovsky, E. and Sergeeva, A. (2006). Numerical modeling of the kdv random wave field. *European Journal of Mechanics-B/Fluids*, 25(4):425–434.

[Perlin et al., 2013] Perlin, M., Choi, W., and Tian, Z. (2013). Breaking waves in deep and intermediate waters. *Annual Review of Fluid Mechanics*, 45:115–145.

[Seiffert and Ducrozet, 2016] Seiffert, B. and Ducrozet, G. (2016). Deep water wave-breaking in a High-Order Spectral model. USA.

[Sergeeva et al., 2011] Sergeeva, A., Pelinovsky, E., and Talipova, T. (2011). Nonlinear random wave field in shallow water: variable korteweg-de vries framework. *Natural Hazards and Earth System Sciences*, 11(2):323–330.

[Serio et al., 2005] Serio, M., Onorato, M., Osborne, A., and Janssen, P. (2005). On the computation of the benjamin-feir index. *Il Nuovo Cimento C*, 28(6):893–903.

[Toffoli et al., 2009] Toffoli, A., Benoit, M., Onorato, M., and Bitner-Gregersen, E. (2009). The effect of third-order nonlinearity on statistical properties of random directional waves in finite depth. *Nonlinear Processes in Geophysics*, 16(1):131–139.

[Toffoli et al., 2011] Toffoli, A., Cavaleri, L., Babanin, A., Benoit, M., Bitner-Gregersen, E., Monbaliu, J., Onorato, M., Osborne, A., and Stansberg, C. (2011). Occurrence of extreme waves in three-dimensional mechanically generated wave fields propagating over an oblique current. *Natural Hazards and Earth System Sciences*, 11(3):895–903.

[Toffoli et al., 2010] Toffoli, A., Gramstad, O., Trulsen, K., Monbaliu, J., Bitner-Gregersen, E., and Onorato, M. (2010). Evolution of weakly nonlinear random directional waves: laboratory experiments and numerical simulations. *Journal of Fluid Mechanics*, 664:313–336.
835

836

837

[Toffoli et al., 2006] Toffoli, A., Onorato, M., and Monbaliu, J. (2006). Wave statistics in unimodal and bimodal seas from a second-order model. *European Journal of Mechanics-B/Fluids*, 25(5):649-661.

838

839

[Trulsen et al., 2012] Trulsen, K., Zeng, H., and Gramstad, O. (2012). Laboratory evidence of freak waves provoked by non-uniform bathymetry. *Physics of Fluids*, 24(9):097101.

840

841

842

[Vincent and Briggs, 1989] Vincent, C. L. and Briggs, M. J. (1989). Refraction?diffraction of irregular waves over a mound. *Journal of Waterway, Port, Coastal, and Ocean Engineering*, 115(2):269-284.

843

844

[Viotti and Dias, 2014] Viotti, C. and Dias, F. (2014). Extreme waves induced by strong depth transitions: Fully nonlinear results. *Physics of Fluids*, 26(5):051705.

845

846

[West et al., 1987] West, B. J., Brueckner, K. A., Janda, R. S., Milder, M., and Milton, R. L. (1987). A new numerical method for surface hydrodynamics. 92:11803-11824.

847

848

[Zakharov, 1968] Zakharov, V. (1968). Stability of periodic waves of finite amplitude on the surface of a deep fluid. *J. Appl. Mech. Tech. Phys.*, pages 190-194.

849

850

[Zakharov and Gelash, 2013] Zakharov, V. E. and Gelash, A. (2013). Nonlinear stage of modulation instability. *Physical review letters*, 111(5):054101.

851

852

853

[Zeng and Trulsen, 2012] Zeng, H. and Trulsen, K. (2012). Evolution of skewness and kurtosis of weakly nonlinear unidirectional waves over a sloping bottom. *Natural Hazards and Earth System Sciences*, 12(3):631.