# Simulation of breaking waves using the High-Order Spectral method with laboratory experiments: Wave-breaking onset

Betsy Seiffert, Guillaume Ducrozet, Félicien Bonnefoy

► To cite this version:

Betsy Seiffert, Guillaume Ducrozet, Félicien Bonnefoy. Simulation of breaking waves using the High-Order Spectral method with laboratory experiments: Wave-breaking onset. Ocean Modelling, 2017, 10.1016/j.ocemod.2017.09.006. hal-01612035

HAL Id: hal-01612035

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

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
Simulation of breaking waves using the High-Order
Spectral method with laboratory experiments:
Wave-breaking onset

Betsy R. Seiffert<sup>a,b,c,*</sup>, Guillaume Ducrozet<sup>b</sup>, Félicien Bonnefoy<sup>b</sup>

Abstract

This study investigates a wave-breaking onset criteria to be implemented in the
non-linear potential flow solver HOS-NWT. The model is a computationally ef-
ficient, open source code, which solves for the free surface in a numerical wave
tank using the High-Order Spectral (HOS) method. The goal of this study is to
determine the best method to identify the onset of random single and multiple
breaking waves over a large domain at the exact time they occur. To identify
breaking waves, a breaking onset criteria based on the ratio of local energy
flux velocity to the local crest velocity, introduced by Barthelemy et al. [1], is
selected. The breaking parameter is uniquely applied in the numerical model
in that calculations of the breaking onset criteria ratio are not made only at
the location of the wave crest, but at every point in the domain and at every
time step. This allows the model to calculate the onset of a breaking wave the
moment it happens, and without knowing anything about the wave a priori.
The application of the breaking criteria at every point in the domain and at
every time step requires the phase velocity to be calculated instantaneously ev-
erywhere in the domain and at every time step. This is achieved by calculating
the instantaneous phase velocity using the Hilbert transform and dispersion re-

<footnote><sup>a</sup>Florida Atlantic University, Department of Ocean & Mechanical Engineering, Boca Raton, FL 33431-0991, USA</footnote>

<footnote><sup>b</sup> École Centrale de Nantes, LHEEA, UMR CNRS 6598, Nantes, France</footnote>

<footnote><sup>c</sup>Institut de Recherche Technologique Jules Verne, Bouguenais, France</footnote>

<footnote>*Corresponding author.</footnote>

<footnote>Email address: eseiffert@fau.edu (Betsy R. Seiffert)</footnote>
lation. A comparison between more traditional crest-tracking techniques shows
the calculation of phase velocity using Hilbert transform at the location of the
breaking wave crest provides a good approximation of crest velocity.

The ability of the selected wave breaking criteria to predict single and mul-
tiple breaking events in two dimensions is validated by a series of large-scale
experiments. Breaking waves are generated by energy focusing and modula-
tional instability methods, with a wide range of primary frequencies. Steep
irregular waves which lead to breaking waves, and irregular waves with an en-
ergy focusing wave superimposed are also generated. This set of waves provides
a wide range of breaking-wave strengths, types and scales for validation of the
model. A comparison of calculations made using HOS-NWT with experimental
measurements show that the model is successful at predicting the occurrence
of wave breaking, as well as accurately calculating breaking onset time and
location.

Although the current study is limited to a unidirectional wave field, the success of the wave-breaking model presented provides the basis for application of the model in a multidirectional wave field. By including wave breaking onset with the addition of an appropriate energy dissipation model into HOS-NWT, we can increase the application range of the model, as well as decrease the occurrence of numerical instabilities that are associated with breaking waves in a potential flow solver. An accurate description of the wave field is useful for predicting the dynamic response of offshore vessels and marine renewable energy devices, predicting loads on marine structures and the general physics of ocean waves, for example.

*Keywords:* Wave-breaking, High-Order Spectral method, Laboratory experiments, Numerical wave tank, Wave focusing

1. Introduction

2 The nonlinear and irregular processes that occur just before a wave breaks

3 make predicting onset difficult. Breaking develops quickly, reportedly as fast
as half of a wave period (Babanin et al. [2], Chalikov and Babanin [3]), and
deep water breaking waves can range from a highly energetic plunging breaker
to a mildly sloping spilling breaker. Grue and Jensen [4] observed "backward
breaking" characterized by a steep negative slope on the back side of the wave
and an asymmetric velocity field behind the wave crest. Therefore, even the
definition of wave-breaking onset itself can be somewhat subjective given these
challenges.

Wave-breaking onset criteria is traditionally classified into one of three categories: geometric, kinematic or dynamic. Perlin et al. [5] provides a recent review of each of the three wave breaking criteria and the strengths a limitations of each.

Geometric criteria use a limiting steepness parameter to determine wave
breaking onset. Babanin et al. [6] claims that waves break not because of a
particular physical mechanism, but because they have reached a limiting steep-
ness and the water surface becomes unstable, regardless of the process behind
steepness growth. Babanin et al. [7] investigated the causes of wave break-
ing in initially monochromatic, two-dimensional waves using both numerical
and experimental methods. Measurements of steepness, skewness and asym-
metry were made for initially monochromatic waves that evolve to breaking
due to Benjamin-Feir instability. The authors conclude that the Stokes’ limit-
ing steepness is the single parameter which determines whether the wave will
break. However laboratory experiments such as Rapp and Melville [8], Tulin
and Waseda [9], Tian et al. [10], and Tian et al. [11] demonstrate that the
limiting steepness is smaller than the Stokes’ limiting steepness and can vary
widely according to breaker type and the method used to generate the breaking
event. Through the numerical calculations of Barthelemy et al. [1] and exper-
imental measurements of Saket et al. [12], it was shown that the variation in
wave steepness at breaking made it unsuitable as a breaking onset parameter.
Furthermore, due to the nature of the wave crests generated during irregular
wave propagation, determining an accurate and repeatable breaking parameter
based on wave steepness that would be applicable to irregular waves presents
many difficulties, and is impractical.

Kinematic breaking criteria is based on the relationship between the hori-
zontal wave crest particle velocity $U$ and the wave phase speed $C$, in that wave
breaking occurs when $U/C \ge 1$. Difficulty arises in this criteria when determin-
ing $U$ and $C$ for highly unsteady and rapidly changing waves. Stansell and Mac-
Farlane [13] experimentally investigated the kinematic criteria using PIV mea-
surements and using three different methods of determining $C$. Even so, they
found $U/C$ to be a maximum of 0.81 for plunging breakers and 0.95 for spilling
breakers. This suggests that the kinematic criteria may underestimate wave
breaking. Wu and Nepf [14] examined kinematic criteria in three-dimensional
experiments. They reported $U/C \ge 1$ successfully predicted breaking waves.
Additionally, they found that for 3D plunging breakers, $U/C \ge 1.5$ and for
spilling breakers $U/C \ge 1$. For modulational wave breaking, Tulin and Lan-
drini [15] suggest the wave will eventually break if the orbital velocity exceeds
the phase velocity.

Researchers such as Melville and Matusov [16], Banner and Peirson [17], Tian et al. [18], and Shemer [19] have used PIV to investigate kinematic criterion with results showing that a wave will break when the ratio of horizontal fluid velocity at the crest to crest speed $U/C$ is in the range between 0.75 and 0.95. In their computations for deep water irregular waves with breaking, Kurnia and van Groesen [20] compared criterion of $U/C = 0.8$ and $U/C = 1$, and found the overall differences in post breaking water surface to be small.

Song and Banner [21] developed a dynamic wave breaking onset threshold
parameter for two-dimensional nonlinear waves based on the growth rate of
the mean local wave energy. The experimental results of Banner and Peirson
[17], Tian et al. [10] and Tian et al. [18] further validate the use of the the
dynamic breaking parameter introduced by Song and Banner [21] but with some
modifications of the choice of local wave number. Recently, Derakhti and Kirby
[22], Derakhti et al. [23] and Derakhti et al. [24] have proposed a modified wave
breaking onset threshold parameter based on the breaking parameter of Song
and Banner [21] for use in their non-hydrostatic numerical model, NHWAVE.
In their numerical calculations, Barthelemy et al. [1] concluded that the breaking parameter of Song and Banner [21] was not applicable to 3D and have proposed a new wave breaking threshold parameter based on the ratio of local energy flux velocity to the local crest velocity. Barthelemy et al. [1] argues that wave-breaking occurs when a waveform is unable to remain stable due to an excess of local wave energy flux, which can be caused by wave energy focusing, wind input or wave-current interaction, for example. The local energy flux velocity is calculated as the ratio of the local energy flux vector to local energy density. Some researchers, Banner et al. [25] and Barthelemy et al. [1], for example, have observed a slowdown in crest velocity just before breaking onset, and based on this observation, Barthelemy et al. [1] defines a breaking criterion as $B_x = F_x/(Ec_x)$ where $F_x$ and $c_x$ are the projections in the direction of wave propagation (in this case, the x-axis) of energy flux and crest velocity, respectively, and $E$ is the local energy density. Barthelemy et al. [1] performed numerical calculations using the WSIM Boundary Element Method code (Grilli et al. [26], Grilli et al. [27]) and determined that wave breaking occurs in the narrow range of $0.85 < B_x < 0.86$ for both 2D and 3D, and subsurface analysis of this parameter found the maximum value for $B_x$ at the free surface. When comparing the parameter $B_x$ with local wave steepness ($S_c$) just before breaking, $B_x$ was found to have a much smaller threshold range when it came to distinguishing between a breaking versus a non-breaking wave than steepness. Saket et al. [12] investigated the breaking criterion of Barthelemy et al. [1] experimentally and found $B_x$ to be $0.84 \pm 0.016$, which is in good agreement with Barthelemy et al. [1]. Additionally, Saket et al. [12] found this threshold parameter to be robust for different classes of wave groups and wave steepness.

It is noted that at the free surface where pressure can be approximated as zero, the energy flux velocity reduces to the water particle velocity, reducing this dynamic breaking criterion to a kinematic criterion. This is different from the traditional kinematic breaking criterion which states $U/C \ge 1$. It is advantageous to use the breaking criterion of Barthelemy et al. [1] as it is straightforward to calculate in numerical models, and according to [1], is applicable to 2D and
3D, and has a narrow threshold. However, as stated by Banner and Peregrine [28] and Kurnia and van Groesen [20] since the wave is not of permanent form, the calculation of the wave crest velocity is non-trivial and must be calculated carefully. Barthelemy et al. [1] measured a forward and backward leaning cycle of a wave crest as it moves through a wave train, and it is at the point close to when the crest is moving at its slowest speed that the wave breaks. Therefore a method of calculating the crest velocity that can accurately capture this slow-down is essential to accurately predict breaking onset. However, calculation of the crest velocity is not necessarily straightforward when computing for multiple crests over a large domain, where the location of each crest must be determined.

The numerical model used in this study, HOS-NWT, is computationally efficient, open source code that solves for highly nonlinear wave propagation in a numerical wave tank using the High-Order Spectral (HOS) method. HOS-NWT and HOS-ocean, have demonstrated success in simulating highly accurate wave fields and wave generation, dispersion over flat and variable bathymetry, as well as nonlinear wave interaction and extreme sea states at low computational cost in part due to FFT in the HOS method (Ducrozet et al. [29], Bonnefoy et al. [30], Ducrozet et al. [31], Gouin et al. [32], Ducrozet et al. [33]).

Due to the assumptions of solving for potential flow, the HOS-NWT assumes a single-valued free surface and therefore cannot produce breaking waves. If a wave surpasses some limiting threshold that would cause a wave to break in nature, the HOS solver will instead calculate high-frequency oscillations that can lead to numerical instabilities, as discussed in Seiffert and Ducrozet [34], for example. The goal of implementing a wave-breaking mechanism into the model is to approximate the broken free surface as a single value. The first step to implementing a wave breaking mechanism into HOS-NWT is to accurately identify the time and location of a breaking wave.

To validate our model, a series of experiments is performed to generate breaking waves with different breaker types and strengths and at different ranges of primary frequencies. Much of the existing laboratory experiments on wave breaking is at a small scale, may not include the exact location and time of
breaking, or does not include breaking in irregular wave time series. Breaking waves for these experiments will be generated using energy focusing, modulational instability, steep irregular waves, and irregular waves with an energy focused wave superimposed. Particular attention is paid to accurately measuring the time and location of wave breaking.

## 2. Numerical model description

The High-Order Spectral method was originally developed to solve for the highly nonlinear wave propagation in an open ocean over flat bottom (Dommermuth and Yue [35], West et al. [36]). The open source solver HOS-ocean is an example of such model which has been widely validated in different configurations (Ducrozet et al. [29], Ducrozet et al. [33], Ducrozet et al. [37]). This model has been extended to more complex configurations, such as varying bathymetry (Gouin et al. [38] and Gouin et al. [32]), and the model discussed here, HOS-NWT, which is a numerical wave tank, (Bonnefoy et al. [30], and Ducrozet et al. [31]). Within HOS-NWT, the numerical domain has perfectly reflective side walls, a wavemaker to generate waves, and an absorbing zone to mimic the beach in a wave basin and to prevent reflections from the end wall.

Inside the domain the fluid is assumed to be incompressible and inviscid with irrotational flow, and therefore is governed by the Laplace equation. Following Zakharov [39], the free surface boundary conditions are written in terms of the surface quantities $\eta$ and the surface potential $\tilde{\phi}(x, t) = \phi(x, z = \eta, t)$, where $x$ is the horizontal $(x, y)$ vector. The kinematic and dynamic free surface boundary conditions are therefore written as

$$ \frac{\partial \eta}{\partial t} = (1 + |\nabla \eta|^2) \frac{\partial \phi}{\partial z} - \nabla \tilde{\phi} \cdot \nabla \eta \quad (1) $$

$$ \frac{\partial \tilde{\phi}}{\partial t} = -g\eta - \frac{1}{2} |\nabla \tilde{\phi}|^2 + \frac{1}{2}(1 + |\nabla \eta|^2) \left(\frac{\partial \phi}{\partial z}\right)^2 \quad (2) $$

on $z = \eta(x, t)$, where $\nabla$ is the horizontal gradient.

These nonlinear free surface boundary conditions enable the time-stepping of the two unknowns ($\eta$ and $\tilde{\phi}$) once the HOS method has evaluated the free surface
vertical velocity. The procedure relies on series expansion in wave steepness $\epsilon$ up to the so-called HOS order M, which is arbitrary.

The fixed boundary conditions use a no-flux condition, $\nabla\phi \cdot \mathbf{n} = 0$, where $\mathbf{n}$ is the vector normal to the boundary. These conditions define a spectral basis on which the velocity potential in the whole volume will be expanded. At the same time, surface quantities are also expressed on a spectral basis allowing the use of fast Fourier transforms (FFTs). The discretization is associated to the number of natural modes of the domain used for this expansion.

The specificity of HOS-NWT is in the presence of the following wavemaker boundary condition:

$$ \frac{\partial X}{\partial t} = \frac{\partial \phi}{\partial x} - (\nabla_v X) \cdot (\nabla_v \phi) \quad \text{on } x = X(y, z, t) \qquad (3) $$

where $\nabla_v = (\partial_y, \partial_z)$ is the vertical gradient. The motion of the wavemaker $X(y, z, t)$ is decomposed as $X(y, z, t) = f_v(z)\tilde{X}(y, t)$, where $f_v$ is a function describing the vertical geometry.

A specific procedure has been developed for an efficient and accurate description of the wavemaker, trying to preserve the pseudo-spectral formalism of the original HOS method. The main idea, following Agnon and Bingham [40], is to use the concept of additional potential. The potential $\phi$, solution of the whole problem (generation and propagation), is separated into two parts: i) one is the potential in the fixed geometry tank (no wavemaker) with its free surface and ii) one additional potential only accounting for the wavemaker. The latter is evaluated with perturbation expansions, defining the wavemaker order of nonlinearity.

In addition, an absorbing zone is modeled in the numerical model through a local modification of the free surface dynamic boundary condition. Further details on the model HOS-NWT can be found in Bonnefoy et al. [30] and Ducrozet et al. [31].
## 3. Experimental design

Experiments are conducted in the 2D wave flume at in the Hydrodynamics, Energetics & Atmospheric Environment Lab (LHEEA) at École Centrale de Nantes (ECN). The flume is 148 m long, 5 m wide and 3 m deep with a parabolic shaped absorbing beach approximately 8 m long.

Breaking waves are generated using energy focusing and modulational instability techniques, two widely used methods of generated controlled breaking waves in the laboratory. In addition, irregular waves with large steepness as well as irregular waves with energy focusing frequency and amplitude components superimposed will be generated to produce breaking waves in an irregular wave series. This will provide a wide range of different types of breakers and breaker strengths. In addition, these waves will be ran with several ranges of primary frequencies, generating breakers at different scales.

Measurements for surface elevation will be collected by a series of 14 resistance type wave gauges positioned in the region of breaking, as well as up wave and down wave of the breaking regions. Measurements are made at 128 Hz, and each wave was ran five times to assure repeatability of measurements. A repeatability analysis of surface elevation and frequency/amplitude spectra show wave probe measurements are well within 5% error. A butterworth filter is applied to the measurements during post-processing to filter out high frequency noise above 20 Hz.

Measurements of the wave breaking time and location are made by video recordings taken at the wave breaking location. A 4 meter measuring stick is affixed to the back wall of the flume to measure breaking location, and breaking time is measured by an LED light placed in the camera field of view that turns on when the wavemaker motion begins. Breaking onset is defined as the point at which the water surface begins to deform and become unstable. The location of breaking onset is measured as the horizontal distance from this point to the mean wavemaker position. One wave was recorded 5 times to ensure repeatability of breaking onset time and location measurements. Video analysis shows that the
210 breaking onset location and time measurements are also well within 5% error.

## 3.1. Breaking wave generation

211 Waves are generated with a flap-type wavemaker. The height above the flume bottom of the flap hinge is 0.48 m and the maximum stroke at the water depth 2.9 m is approximately ±0.3 m. A ramp of 4 s is applied to the beginning of the signal. A linear transfer function is used to convert the target free surface elevation to the corresponding wavemaker motion.

### 3.1.1. Energy focusing

218 Single plunging breaking waves are generated using a wave energy focusing technique described in Rapp and Melville [8], for example. Surface elevation at the wavemaker is described as

$$ \eta(0, t) = \sum_{n=1}^{N} a_n \cos(-k_n x_b - \omega_n(t - t_b)) \quad (4) $$

221 where $a_n$ is the amplitude of the nth wave component, $k_n$ is the wavenumber, $\omega_n = 2\pi f_n$ is the angular frequency, N is the number of frequency components, and $x_b$ and $t_b$ are the input breaking location and time. It is noted that due to nonlinear features associated with wave propagation, the input $x_b$ and $t_b$ are not exact, and the breaking wave tends to be within one period or wavelength of these values. The principle behind generating breaking waves using energy focusing is that a nonlinear superpositioning of waves of different amplitude and phase will converge at a determined location, causing a focusing of energy. Ideal generated waves using the energy focusing technique in deep water according to Sriram et al. [41] and Rapp and Melville [8] will have a broad spectrum as to produce a single breaking wave, and have N large enough to simulate a continuous spectrum.

### 3.1.2. Modulational instability

233 Multiple breaking waves are generated using modulational instability, a description of which can be found in Tulin and Waseda [9], for example. Surface elevation at the wavemaker is described by
$$ \eta(t) = a_0 \cos(\omega_0 t) + b \cos \left( \omega_1 t - \frac{\pi}{4} \right) + b \cos \left( \omega_2 t - \frac{\pi}{4} \right). \quad (5) $$

a₀ and ω₀ are the amplitude and angular frequency of the carrier wave, and b, ω₁ and ω₂ are the amplitude and frequencies of the sideband perturbations, where ω₁ = ω₀ - Δω/2, and ω₂ = ω₀ + Δω/2. The principle behind modulational instability is that nonlinear energy transfer from the center frequency to the sideband frequencies will cause waves to focus their energy and eventually break.

## 3.2. Tian, et al. (2012) scaled experiments

To have confidence in experimental measurements and to determine whether the wave breaking model will perform at a large scale, the experiments of Tian et al. [11] are reproduced at a 4:1 scale using Froude's similitude law. Tian's experiments are conducted in a wave flume 15 m long, 1.5 m wide and with a water depth of 0.62 m. Wave probes measured surface elevation at 14 locations along the flume at distances measured from the mean position of the wavemaker. Experimental setup parameters for Tian's experiments and the corresponding LHEEA experiments are given in Table 1. The water depth in the LHEEA flume remains constant at 2.9 m. However, since these are deep water waves, effects of not having an exact (4:1) scaling for water depth should be negligible.

In the Tian et al. [11] experiments, breaking waves are generated using a piston-type wave maker using both energy focusing and modulational instability. Waves generated using energy focusing have a center frequency of $f_c = 1.7$ Hz, a frequency bandwidth of $\Delta f = 1.4$ Hz and $N = 128$ frequency components. The breaking time and location are set as $t_b = 25.0$ s and $x_b = 7.0$ m. The frequency components are generated using these parameters, which are then summed according to Eqn. 3.1.2 to get the desired surface elevation. The initial input steepness $\epsilon_n = k_n a_n$ and corresponding parameters for wave generation in the LHEEA experiments are given in Table 2.

Table 3 gives the parameters of initial frequency $f_0 = \omega_0/2\pi$, ratio of frequency bandwidth to initial angular frequency ($\Delta\omega/\omega_0$) and initial wave steepness ($\epsilon = k_0 a_0$) used for waves generated using modulational instability, for both
<table><caption>Table 1: Experimental parameters from Tian et al. [11] and corresponding parameters for LHEEA experiments. *Length of flume corresponding to Tian experiments. **Water depth in the LHEEA flume is kept constant at 2.9 m.</caption><thead><tr><th>Parameter</th><th>Tian et al. [11] (1:1 scale)</th><th>LHEEA (4:1 scale)</th></tr></thead><tbody><tr><td>Flume length</td><td>15 m</td><td>60 m*</td></tr><tr><td>Water depth (h)</td><td>0.62 m</td><td>2.9 m**</td></tr><tr><td>Wave probes</td><td>[1.83, 2.53, 3.32, 4.04, 4.77<br>5.38, 5.99, 6.6, 7.22 7.89<br>8.57, 9.42, 10.23, 11.06] m</td><td>[7.32 10.12 13.28 16.16 19.08<br>21.52 23.96 26.4 28.88 31.56<br>34.28 37.68 40.92 44.24] m</td></tr></tbody></table>

<table><caption>Table 2: Input parameters for Tian et al. [11] and LHEEA experiments for waves generated using energy focusing. N = 128. Wave probe measurements will be collected at 128 Hz for 80 s.</caption><thead><tr><th rowspan="2">Wave</th><th colspan="5">Tian (1:1)</th><th colspan="5">LHEEA (4:1)</th></tr><tr><th>f<sub>c</sub> (Hz)</th><th>Δf (Hz)</th><th>ε<sub>n</sub></th><th>t<sub>b</sub> (s)</th><th>x<sub>b</sub> (m)</th><th>f<sub>c</sub> (Hz)</th><th>Δf (Hz)</th><th>ε<sub>n</sub></th><th>t<sub>b</sub> (s)</th><th>x<sub>b</sub> (m)</th></tr></thead><tbody><tr><td>EF1</td><td>1.7</td><td>1.4</td><td>0.0020</td><td>25.0</td><td>7.0</td><td>0.85</td><td>0.7</td><td>0.0020</td><td>50.0</td><td>28.0</td></tr><tr><td>EF2</td><td>1.7</td><td>1.4</td><td>0.0032</td><td>25.0</td><td>7.0</td><td>0.85</td><td>0.7</td><td>0.0032</td><td>50.0</td><td>28.0</td></tr><tr><td>EF3</td><td>1.7</td><td>1.4</td><td>0.0045</td><td>25.0</td><td>7.0</td><td>0.85</td><td>0.7</td><td>0.0045</td><td>50.0</td><td>28.0</td></tr><tr><td>EF4</td><td>1.7</td><td>1.4</td><td>0.0058</td><td>25.0</td><td>7.0</td><td>0.85</td><td>0.7</td><td>0.0058</td><td>50.0</td><td>28.0</td></tr></tbody></table>

the Tian and LHEEA experiments. The amplitude ratio for each of these waves is $b/a_0 = 0.5$. Energy focusing experiments are run for 80 s and modulational instability experiments are run for 100 s.

### 3.3. Large-scale wave breaking experiments

Breaking waves for a series of experiments will be generated using energy focusing, modulational instability, irregular waves with large wave steepness and irregular waves with added frequency components to generate a focused wave within the time series. A range of central or peak frequencies, frequency
<table><caption>Table 3: Input parameters for Tian et al. [11] and LHEEA experiments for waves generated using modulational instability. $b/a_0 = 0.5$. Wave probe measurements will be collected at 128 Hz for 100 s.</caption><thead><tr><th rowspan="2">Wave group</th><th colspan="3">Tian (1:1)</th><th colspan="3">LHEEA (4:1)</th></tr><tr><th>f<sub>0</sub> (Hz)</th><th>Δω/ω<sub>0</sub></th><th>ε = k<sub>0</sub>a<sub>0</sub></th><th>f<sub>0</sub> (Hz)</th><th>Δω/ω<sub>0</sub></th><th>ε = k<sub>0</sub>a<sub>0</sub></th></tr></thead><tbody><tr><td>BFI1710</td><td>1.7</td><td>0.166</td><td>0.10</td><td>0.85</td><td>0.166</td><td>0.10</td></tr><tr><td>BFI1712</td><td>1.7</td><td>0.200</td><td>0.12</td><td>0.85</td><td>0.200</td><td>0.12</td></tr><tr><td>BFI1716</td><td>1.7</td><td>0.266</td><td>0.16</td><td>0.85</td><td>0.266</td><td>0.16</td></tr><tr><td>BFI1720</td><td>1.7</td><td>0.330</td><td>0.20</td><td>0.85</td><td>0.330</td><td>0.20</td></tr></tbody></table>

bandwidths and initial wave steepness will be tested to obtain a range of breaker size, type and strength. Wave gauges for these experiments are positioned at 7.32, 14.5, 20.00, 21.20, 22.41, 23.60, 24.69, 26.05, 27.94, 30.06, 34.26, 37.88, 43.66 and 50.23 m from the mean position of the wavemaker.

Parameters for breaking waves generated using the energy focusing technique are given in Table 4. The number of frequency components is $N = 64$, and the input breaking time and location are $t_b = 50$ s and $x_b = 22$ m for EF6-54 and EF8-72, and $x_b = 24$ m for EF7-64. Surface elevation is recorded for 100 s, stopping well before the initial wave energy has reflected off the absorbing beach and returned to the last wave gauge (recalling that the wave flume is 148 m long.)

The wave parameters used for generating breaking waves using modulational instability are given in Table 5. $Δω/ω_0$ and $ε = k_0a_0$ are chosen following Fig. 17 in Tulin and Waseda [9] as the values which correspond to the range of frequencies where breaking occurs at a location close to the wavemaker. The parameter $b/a_0 = 0.5$ is selected so that waves are breaking in a region close to the wavemaker but not too close. Surface elevation measurements are recorded for 150 s.

Irregular waves are generated characterized by the JONSWAP spectrum with different peak periods ($T_p$) and initial steepness (given as the ratio of signif-
<table><caption>Table 4: Input parameters for waves generated using energy focusing are given along with the type of breaker and breaking time and location measured during experiments and calculated using HOS-NWT. The input breaking time and location are $t_b = 50$ s and $x_b = 22$ m for EF6-54 and EF8-72, and $x_b = 24$ m for EF7-64. * denotes an unbroken wave.</caption><thead><tr><th rowspan="2">Wave</th><th rowspan="2">f<sub>c</sub> (Hz)</th><th rowspan="2">Δf (Hz)</th><th rowspan="2">ε<sub>n</sub></th><th rowspan="2">Type</th><th colspan="2">Experiments</th><th colspan="2">HOS-NWT</th></tr><tr><th>x<sub>b</sub> (m)</th><th>t<sub>b</sub> (s)</th><th>x<sub>b</sub> (m)</th><th>t<sub>b</sub> (s)</th></tr></thead><tbody><tr><td rowspan="3">EF6-54</td><td rowspan="3">0.6</td><td rowspan="3">0.54</td><td>0.0047*</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>0.0048</td><td>spiller</td><td>23.61</td><td>50.38</td><td>22.75</td><td>50.03</td></tr><tr><td>0.0050</td><td>plunger</td><td>22.70</td><td>49.96</td><td>22.48</td><td>49.86</td></tr><tr><td rowspan="3">EF7-64</td><td rowspan="3">0.7</td><td rowspan="3">0.64</td><td>0.0051</td><td>plunger</td><td>22.53</td><td>49.9</td><td>22.38</td><td>49.8</td></tr><tr><td>0.0048</td><td>spiller</td><td>24.66</td><td>50.06</td><td>24.55</td><td>49.97</td></tr><tr><td>0.0054</td><td>plunger</td><td>21.83</td><td>47.26</td><td>21.68</td><td>47.18</td></tr><tr><td rowspan="3">EF8-72</td><td rowspan="3">0.8</td><td rowspan="3">0.72</td><td>0.0054</td><td>spiller</td><td>25.14</td><td>50.14</td><td>25.12</td><td>50.2</td></tr><tr><td>0.0060</td><td>plunger</td><td>21.74</td><td>47.14</td><td>21.57</td><td>47.04</td></tr><tr><td>0.0062</td><td>plunger</td><td>20.44</td><td>47.52</td><td>20.30</td><td>47.41</td></tr><tr><td></td><td></td><td></td><td>0.0070</td><td>plunger</td><td>20.92</td><td>47.68</td><td>20.67</td><td>47.55</td></tr><tr><td></td><td></td><td></td><td>0.0070</td><td>sm. spiller</td><td>24.46</td><td>50.70</td><td>-</td><td>-</td></tr></tbody></table>

icant wave height to peak wave length, $H_s/L_p$). For these waves the peakedness factor $\gamma$ remains constant at $\gamma = 3.3$. A second type of irregular wave will be generated with an energy focused wave superimposed on an irregular wave. This is accomplished by adding the frequency components for an energy focused wave to the components in a JONSWAP spectrum. This is similar to the strategy for developing an extreme wave found in Ducrozet et al. [29]. The parameters for these waves are given in Table 6. Surface elevation measurements for irregular wave experiments are recorded for 150 s.

## 4. Calculating wave-breaking onset

To calculate wave breaking onset, the criterion described by Barthelemy et al. [1] is used which is based on the ratio of local energy flux velocity to the
<table><caption>Table 5: Input parameters for waves generated using modulational instability. Type of breaking and breaking time and location measured during experiments and calculated using HOS-NWT is given for the first 10 recorded breaking waves.</caption><thead><tr><th rowspan="2">Wave</th><th rowspan="2">f<sub>0</sub> = ω<sub>0</sub>/2π (Hz)</th><th rowspan="2">Δω/ω<sub>0</sub></th><th rowspan="2">ε = k<sub>0</sub>a<sub>0</sub></th><th rowspan="2">Type</th><th colspan="2">Experiments</th><th colspan="2">HOS-NWT</th></tr><tr><th>x<sub>b</sub> (m)</th><th>t<sub>b</sub> (s)</th><th>x<sub>b</sub> (m)</th><th>t<sub>b</sub> (s)</th></tr></thead><tbody><tr><td rowspan="10">MI0719</td><td rowspan="10">0.7</td><td rowspan="10">0.317</td><td rowspan="10">0.19</td><td>spiller</td><td>33.17</td><td>32.66</td><td>31.66</td><td>32.09</td></tr><tr><td>plunger</td><td>32.95</td><td>41.14</td><td>32.85</td><td>41.13</td></tr><tr><td>spiller</td><td>32.14</td><td>49.44</td><td>31.77</td><td>49.25</td></tr><tr><td>spiller</td><td>31.44</td><td>57.70</td><td>31.03</td><td>57.53</td></tr><tr><td>spiller</td><td>34.74</td><td>60.42</td><td>34.09</td><td>60.19</td></tr><tr><td>spiller</td><td>31.48</td><td>66.50</td><td>29.90</td><td>65.61</td></tr><tr><td>spiller</td><td>32.56</td><td>68.18</td><td>33.02</td><td>68.32</td></tr><tr><td>spiller</td><td>32.13</td><td>76.60</td><td>32.16</td><td>76.56</td></tr><tr><td>spiller</td><td>31.38</td><td>84.90</td><td>31.01</td><td>84.64</td></tr><tr><td>spiller</td><td>34.50</td><td>87.60</td><td>34.35</td><td>87.45</td></tr><tr><td rowspan="10">MI0720</td><td rowspan="10">0.7</td><td rowspan="10">0.330</td><td rowspan="10">0.20</td><td>spiller</td><td>26.77</td><td>27.00</td><td>26.20</td><td>26.92</td></tr><tr><td>plunger</td><td>29.86</td><td>29.76</td><td>28.83</td><td>29.40</td></tr><tr><td>spiller</td><td>34.20</td><td>33.06</td><td>32.42</td><td>32.37</td></tr><tr><td>spiller</td><td>28.03</td><td>36.20</td><td>27.75</td><td>36.09</td></tr><tr><td>spiller</td><td>32.05</td><td>39.26</td><td>31.26</td><td>38.99</td></tr><tr><td>spiller</td><td>27.70</td><td>44.58</td><td>28.04</td><td>44.83</td></tr><tr><td>plunger</td><td>31.51</td><td>47.64</td><td>30.87</td><td>47.33</td></tr><tr><td>spiller</td><td>26.91</td><td>52.82</td><td>-</td><td>-</td></tr><tr><td>spiller</td><td>27.05</td><td>61.48</td><td>27.46</td><td>61.71</td></tr><tr><td>spiller</td><td>26.98</td><td>70.08</td><td>27.01</td><td>70.08</td></tr></tbody></table>
<table><caption>Table 6: Input parameters for irregular waves and irregular waves with a superimposed energy focusing wave are given along with the type of breaker and breaking time and location measured during experiments and calculated using HOS-NWT. Additional input parameters for the focusing component are N = 64, $x_b = 22$ m and $t_b = 50$ s. *IRR-4 did not contain any breaking waves.</caption><thead><tr><th rowspan="2">Wave</th><th colspan="2">Irregular</th><th colspan="3">Focusing</th><th rowspan="2">Type</th><th colspan="2">Experiments</th><th colspan="2">HOS-NWT</th></tr><tr><th>T<sub>p</sub></th><th>H<sub>s</sub>/L<sub>p</sub></th><th>f<sub>c</sub></th><th>δf</th><th>ε<sub>n</sub></th><th>x<sub>b</sub> (m)</th><th>t<sub>b</sub> (s)</th><th>x<sub>b</sub> (m)</th><th>t<sub>b</sub> (s)</th></tr></thead><tbody><tr><td>IRR-4*</td><td>2.1</td><td>0.030</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td rowspan="2">IRR-1</td><td rowspan="2">2.2</td><td rowspan="2">0.037</td><td rowspan="2">-</td><td rowspan="2">-</td><td rowspan="2">-</td><td>spiller</td><td>16.39</td><td>55.48</td><td>-</td><td>-</td></tr><tr><td>spiller</td><td>16.40</td><td>76.54</td><td>16.14</td><td>76.36</td></tr><tr><td rowspan="3">IRR-2</td><td rowspan="3">2.2</td><td rowspan="3">0.040</td><td rowspan="3">-</td><td rowspan="3">-</td><td rowspan="3">-</td><td>sm. spiller</td><td>16.59</td><td>34.36</td><td>-</td><td>-</td></tr><tr><td>sm. spiller</td><td>25.88</td><td>66.32</td><td>-</td><td>-</td></tr><tr><td>sm. spiller</td><td>16.12</td><td>98.86</td><td>16.41</td><td>98.98</td></tr><tr><td rowspan="5">IRR-5</td><td rowspan="5">2.1</td><td rowspan="5">0.035</td><td rowspan="5">-</td><td rowspan="5">-</td><td rowspan="5">-</td><td>sm. spiller</td><td>30.00</td><td>70.92</td><td>29.90</td><td>70.70</td></tr><tr><td>spiller</td><td>32.06</td><td>73.22</td><td>33.70</td><td>73.88</td></tr><tr><td>sm. spiller</td><td>20.02</td><td>89.54</td><td>19.71</td><td>89.35</td></tr><tr><td>sm. spiller</td><td>32.54</td><td>126.74</td><td>33.69</td><td>127.05</td></tr><tr><td>plunger</td><td>32.19</td><td>147.74</td><td>31.46</td><td>147.4</td></tr><tr><td>IRREF-1</td><td>2.2</td><td>0.02</td><td>0.7</td><td>0.64</td><td>0.005</td><td>plunger</td><td>22.27</td><td>49.72</td><td>22.76</td><td>49.93</td></tr><tr><td rowspan="2">IRREF-2</td><td rowspan="2">2.2</td><td rowspan="2">0.02</td><td rowspan="2">0.7</td><td rowspan="2">0.64</td><td rowspan="2">0.0055</td><td>plunger</td><td>20.05</td><td>47.24</td><td>20.33</td><td>47.42</td></tr><tr><td>spiller</td><td>23.32</td><td>50.14</td><td>22.97</td><td>49.98</td></tr><tr><td rowspan="7">IRREF-3</td><td rowspan="7">2.2</td><td rowspan="7">0.03</td><td rowspan="7">0.7</td><td rowspan="7">0.64</td><td rowspan="7">0.0055</td><td>plunger</td><td>22.40</td><td>147.74</td><td>19.77</td><td>145.20</td></tr><tr><td>spiller</td><td>20.10</td><td>47.26</td><td>19.85</td><td>47.09</td></tr><tr><td>spiller</td><td>22.90</td><td>49.82</td><td>-</td><td>-</td></tr><tr><td>spiller</td><td>23.17</td><td>49.94</td><td>-</td><td>-</td></tr><tr><td>spiller</td><td>29.18</td><td>52.02</td><td>29.31</td><td>52.01</td></tr><tr><td>plunger</td><td>21.09</td><td>145.74</td><td>21.29</td><td>145.80</td></tr><tr><td>spiller</td><td>23.66</td><td>148.3</td><td>24.34</td><td>148.61</td></tr></tbody></table>
local crest velocity. Recall, at the free surface the ratio of the local energy flux in the direction of wave propagation $F_x$ to local energy density $E$ reduces to the water particle velocity under the crest $U_x$. In two-dimensions, $U_x$ can be written in terms of the velocity potential at the free surface as,

$$U_x = \frac{\partial \phi}{\partial x}(z = \eta) \qquad (6)$$

and therefore, the breaking criterion can be written as follows:

$$B_x = \frac{F_x/E}{C_x} = \frac{U_x}{C_x} = \frac{\partial \phi}{\partial x} \frac{1}{C_x} > \text{threshold}, \qquad (7)$$

where $C_x$ is the crest velocity in the direction of wave propagation. A threshold of 0.85 is used to determine wavebreaking onset in calculations for this study, as suggested by the results of Barthelemy et al. [1] and Saket et al. [12].

To identify a random breaking wave at the instant of breaking onset, the ratio given by Eqn. (7) is computed at every point in the domain, and at every time step. Of course a crest cannot exist at every point in the domain at a given instant. Therefore at a given time step, to compute this ratio instantaneously over the whole domain, the instantaneous local phase velocity must be evaluated every where in the domain. At the location of a wave crest, this local phase velocity becomes an approximation of the crest velocity, and therefore this ratio approximates the breaking parameter of [1]. Furthermore, the breaking parameter is naturally at its maximum at the crest of a wave and nowhere else. This approximation is necessary due to the difficulties in calculating the crest velocity, which will be discussed below. We define the instantaneous local phase velocity as

$$C_x(x, t) = \sqrt{\frac{g \cdot \tanh(k(x, t) \cdot h)}{k(x, t)}} \qquad (8)$$

where $h$ is the water depth and the local wave number $k(x, t)$ is given by

$$k(x, t) = \frac{1}{\eta^2 + \mathbb{H}^2[\eta]} \left( \eta \frac{\partial \mathbb{H}[\eta]}{\partial x} - \mathbb{H}[\eta] \frac{\partial \eta}{\partial x} \right) \qquad (9)$$
where $\mathbb{H}[\eta(x, t)]$ is the partial Hilbert transform with respect to $x$. It is possible through some numerical error that a value of $k$ may be briefly negative, and as this would cause large errors in the calculation of the breaking parameter, any such negative values were filtered out. The calculation of crest velocity using the phase velocity approximation in a kinematic breaking condition has been used by Stansell and MacFarlane [13] and Kurnia and van Groesen [20], for example.

## 4.1. Calculation of wave crest velocity

It is of interest to assess whether the approximation of the crest velocity using a local phase velocity at the location of the crest is as accurate as using a more traditional crest tracking technique when calculating breaking waves. Two time series generated during laboratory experiments are chosen for investigation, EF6-54 with $\epsilon = 0.0048$, and IRREF-1, as they each produce only a single breaking wave. A numerical wave flume of 100 m is chosen for the calculations as it is long enough to cover each of the wave gauges in experiments sufficiently beyond the wave gauge furthest from the wavemaker, located at a distance of 50.25 m. Calculations are performed using 1025 modes, with an HOS order mHOS = 5 and wavemaker nonlinearity order 2, according to a sensitivity study which will be discussed in section 4.2. The time period that calculations are allowed to continue after wave breaking is short enough for these particular breaking waves, that instabilities due to the “broken” wave have yet to develop in the simulations. For larger breaking waves than the ones chosen for analysis in the figures, instabilities can develop quite quickly.

This study was conducted during post-processing using the known time and location of the breaking wave measured from experiments. Crest velocity and breaking parameter are calculated using three different methods. In the first method, crest velocity is calculated using a simple crest tracking method. The horizontal position of the wave crest $x_c(t)$ at time $t$ is measured at several time steps before breaking, and the velocity of the crest is taken as the distance the crest has traveled over time, or $C_{x1} = dx_c/dt$. This crest velocity $C_{x1}$ is used to calculate its corresponding breaking parameter, given as $B_{x1}$.
The following equation is given in Fedele [42], and Fedele et al. [43] as the exact crest speed in terms of the surface derivatives:

$$C_{x2} = -\frac{\eta_{,xt}}{\eta_{,xx}}, \qquad (10)$$

where $C_{x2}$ is the two-dimensional horizontal crest velocity in the x-direction, $\eta$ is the surface elevation, and the comma indicates partial derivatives, i.e. $\eta_{,x} = \partial\eta/\partial x$. Equation (10) derives from the condition that $\eta_{,x}(x_c(t), t) = 0$ at the location of the crest, $x_c(t)$ at time $t$. Crest velocity of the known breaking wave calculated using this equation ($C_{x2}$) is used to calculate its corresponding breaking parameter, given as $B_{x2}$.

To compare with the breaking parameter computed in numerical calculations, the maximum breaking parameter is calculated over several time steps before and after breaking using Eqns. (8) and (9) to calculate the phase velocity, and is given as $B_{x3}$. The maximum $B_{x3}$ naturally follows the wave crest and therefore, the local phase velocity used to calculate $B_{x3}$ is given as $C_{x3}$.

Figures 1(a) and 2(a) show calculations for each of the crest velocities described above, for the simulated waves EF6-54 with $\epsilon = 0.0048$, and IRREF-1, respectively. These time series are plotted for a short time window just before and after wave breaking, where $t_{br}$ is breaking onset time measured during experiments. Each of these methods predicts a slowdown of crest velocity near the time of breaking, although $C_{x3}$ appears to over estimate this slowdown in the case of the energy focused wave. Overall, $C_{x3}$ is able to give a good approximation of the crest velocity when compared with $C_{x1}$ and $C_{x2}$.

Some large spikes are observed when calculating $C_{x1}$ and $C_{x2}$. These large spikes are the result of relatively large shifts in the location of the peak of the wave of interest. The wave crest is located by finding the position of the maximum surface elevation at a given time step. However, the position of the maximum surface elevation can rapidly change as wave energy focuses during a breaking event, which can be observed in Fig. 3. This figure plots surface elevation for several time steps before breaking and includes vertical lines indicating where the maximum surface elevation is located at that time step. The rapid
<figure><img src="image_2.png" /><figcaption>Figure 1: (a) Crest velocity is calculated using crest tracking ($C_{x1}$), Eqn.(10) ($C_{x2}$) and approximated using Hilbert transform ($C_{x3}$) for the wave EF6-54 and $\epsilon = 0.0048$. (b) The ratio of $C_{x3}$ to reference velocity is calculated for four reference velocities, $C_0$, $C_{02}$, $C_{03}$, and $C_{04}$. (c) Breaking onset parameter $B_x$ is calculated corresponding to the crest velocities in (a) along with normalized crest velocity $C_{x3}/C_0$, and normalized water particle velocity at the wave crest $U_{x,cr}$.</figcaption></figure>

change in position of the maximum surface elevation causes these large fluctuations in crest velocity. Furthermore, when calculating $C_{x1}$ and $C_{x2}$, significant filtering had to be applied to the data to remove high frequency oscillations due to the numerical evaluation of second order derivatives. This illustrates the level of difficulty involved when using a crest tracking technique to locate a breaking wave. Using crest tracking to identify multiple random breaking waves over a large domain, and in a time domain solver, without knowing ahead of time where they will break, would require significant computational effort and may be impractical. Calculation of the breaking parameter instantaneously and
<figure><img src="image_3.png" /><figcaption>Figure 2: (a) Crest velocity is calculated using crest tracking ($C_{x1}$), Eqn.(10) ($C_{x2}$) and approximated using Hilbert transform ($C_{x3}$) for the wave IRREF-1. (b) The ratio of $C_{x3}$ to reference velocity is calculated for four reference velocities, $C_0$, $C_{02}$, $C_{03}$, and $C_{04}$. (c) Breaking onset parameter $B_x$ is calculated corresponding to the crest velocities in (a) along with normalized crest velocity $C_{x3}/C_0$, and normalized water particle velocity at the wave crest $U_{x,cr}$.</figcaption></figure>

everywhere in the domain means we do not need to identify the peak of every potential breaking wave in the domain and track its location over time, eliminating these computational challenges. Rather, we are able to calculate the exact time and location of a breaking wave at the moment the breaking parameter surpasses the given threshold.

It has been noted by Banner et al. [25], [42] and Barthelemy et al. [1], for example, that the crest slows down to approximately 0.8 times a reference linear wave velocity when the wave is about to break. Here we investigate the relationship between the crest velocity approximation $C_{x3}$ and four different reference
<figure><img src="image_4.png" /><figcaption>Figure 3: Surface elevation is plotted for several time steps up until breaking. The vertical lines represent the location of the maximum surface elevation at that time step.</figcaption></figure>

velocities. The first reference velocity, $C_0$ is calculated by first decomposing the series into its modal spectrum at each time step, and then taking the average value of $k$ that corresponds to the peak energy over time, and then using this value of $k$ to calculate $C_0$ using the linear dispersion relation. $C_{02}$ is calculated by determining the peak frequency estimated from the time series taken at the location of the numerical wavemaker and dispersion relation to obtain $C_{02}$. Using the same time series as in $C_{02}$, $C_{03}$ is calculated using a spectrally weighted wave frequency (as in Banner et al. [25]) where

$$ \bar{f} = \frac{\int fS(f)df}{\int S(f)df} \qquad (11) $$

and $C_{03} = g/2\pi\bar{f}$. Finally, $C_{04} = g/2\pi f_c$, where $f_c$ is simply the input central frequency. The ratio of $C_{x3}$ to each of these reference velocities plotted with
normalized breaking time is given in Figs. 1(b) and 2(b). These figures show ratios of approximately 0.5 to 0.7, smaller than the approximation of 0.8 given by [25], [42] and [1]. For the case of the energy focused wave, it appears an underestimation of the crest velocity approximation may partially contribute to this difference. However, for the irregular wave, all three methods of calculating crest velocity ($C_{x1}$, $C_{x2}$, and $C_{x3}$) are nearly identical, so there must be some other reason the ratios are not closer to 0.8.

In Figs. 1(c) and 2(c), the breaking onset thresholds $B_{x1}$, $B_{x2}$, and $B_{x3}$ are presented along with a solid horizontal line indicating the breaking onset threshold value of 0.85. Also included in these figures are the normalized velocity for $C_{x3}$, and maximum horizontal water particle velocity under the crest $U_{x,cr}$ normalized with the value $U_{x,cr}$ at wave breaking onset.

For the energy focused wave EF6-54 with $\epsilon = 0.0048$, given in Fig. 1(c), it is evident that each of the breaking parameters are successful at predicting wave breaking when compared with experiments, although $B_{x3}$ calculates breaking slightly early. (Recall, $t_{br}$ in the figure is breaking onset for the single breaking wave measured during experiments.) In fact, there is excellent agreement between each of the three breaking parameters for the large non-breaking waves before and after the breaking wave as well.

For the irregular wave IRREF-1, given in Fig. 2(c), there is excellent agreement between each of the breaking parameters and each are successful at predicting breaking onset when compared with experiments. It is noted some irregularities in $C_{x3}$ observed approximately 1.5 wave periods after breaking in Fig. 2(a) do not adversely affect the calculation of the breaking parameter in Fig. 2(c).

Although agreement is good between each of the breaking parameters for both the energy focused wave and the irregular wave in Figs. 1(c) and 2(c), it is worth noting that it required a significant amount of smoothing of $C_{x1}$, $C_{x2}$ and $B_{x1}$, $B_{x2}$ in order to achieve this agreement.

This investigation suggests calculation of wave breaking onset by use of the current breaking parameter which uses an approximation of crest velocity is
successful when compared with crest tracking. Additionally, these figures suggest that it is the combination of a decrease in crest velocity and an increase of water particle velocity under a steepening wave that leads to wave breaking.

## 4.2. Sensitivity study

A sensitivity study is conducted using the breaking wave series EF2 and BFI1716 from Tian et al. [11], to determine the best number of modes to use for the calculation of surface elevation and breaking onset as well as to compare the large scale LHEEA laboratory measurements with measurements made by Tian et al. [11]. For this study, calculations are stopped when the breaking onset threshold is surpassed, indicating wave breaking. Figures 4 and 5 show surface elevation measured during the LHEEA and Tian experiments and calculated using HOS-NWT using 513, 1025 and 2049 modes in the x-direction. The first 7 wave gauges are presented, which are the gauges located upwave to the breaking wave. At 513 modes, the maximum wave number solved is $k_{max} = 16.12 \text{ m}^{-1}$,

<figure><img src="image_5.png" /><figcaption>Figure 4: Surface elevation is computed using HOS-NWT for n1 = 513, 1025 and 2049 modes and compared with measurements made during the LHEEA and Tian experiments for the wave EF2.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Figure 5: Surface elevation is computed using HOS-NWT for n1 = 513, 1025 and 2049 modes and compared with measurements made during the LHEEA and Tian experiments for the wave BFI1716.</figcaption></figure>

at 1025 modes $k_{max}$ = 32.20 m⁻¹, and at 2049 modes $k_{max}$ = 64.37 m⁻¹. This means for the wave EF2 with a center frequency of $f_c$ = 0.85 Hz equating to a center wave number of $k_c$ = 2.9 m⁻¹, calculations made using 513 modes will compute up to the fifth harmonic, calculations made using 1025 modes will calculate the tenth harmonic and calculations made using 2049 will calculate up to the twenty-first harmonic of the center wave number, for example.

The difference observed between each of these calculations for surface elevation is very small. The greatest difference is in calculating breaking onset, where at 513 modes, wave breaking is not calculated until much later when compared with calculations made using 1025 and 2049 modes. A comparison of wave energy is given in Fig. 6, which gives the long time integration of the surface variance ($<\eta^2> = \int_0^{t_{br}} \eta(x, t)^2 dt$) for EF2 up to the point of breaking. The steep drop off in energy is not due to any calculated energy dissipation, but rather due to the time series being halted at wave breaking, so there is very little energy at the seventh wave gauge, for example. According to this figure,
<figure><img src="image_7.png" /><figcaption>Figure 6: The long time integration of the surface variance ($<\eta^2> = \int_0^\tau \eta(x, t)^2 dt$) is computed using HOS-NWT for n1 = 513, 1025 and 2049 modes and compared with measurements made during the LHEEA and Tian experiments for the wave EF2.</figcaption></figure>

it is evident there is convergence at 1025 modes. It is noted that although it appears the lower number of modes, 513, agrees better with experimental measurements, this is somewhat misleading. There is very little difference between calculations for surface elevation for each of the number of modes. The differences in calculations and measurements observed in Fig. 6 are largely due to very small differences in the smaller surrounding waves, that when added up over the entire domain, appear to be more significant than they are.

An interesting observation that can be made from Figs. 4 and 5 is that numerical calculations are nearly identical to surface elevation measured during the LHEEA experiments, whereas there are some small phase differences between calculations and Tian's experiments. This is because waves generated during Tian's experiments are not generated exactly according to Eqn. 8. Rather, some small adjustments of wave amplitude were made during experiments to develop breaking waves without also adjusting wave number *k*, and therefore,
the initial steepness parameter $\epsilon_n$ is not exact, and waves generated in the Tian experiments are not exactly generated according to the theory presented in Sec. 3.1. The LHEEA experiments are generated exactly according to theory, and therefore can be reproduced numerically using the theory from Sec. 3.1, instead of using wave gauge data as input, for example.

Due to this analysis, the number of modes chosen for calculations is 1025, as it is sufficient to capture wave breaking and surface elevation. Note that this choice of discretization is far larger than the one necessary to obtain an accurate solution of wave propagation in the case of non-breaking waves (Ducrozet et al. [31]).

A spatial interpolation is applied to the surface elevation and velocity before breaking onset is calculated at 5 times the original spatial discretization using a cubic spline. This allows the breaking onset to be calculated at a higher spatial resolution, thus increasing the accuracy. The value of 5 was chosen according to a sensitivity study showing convergence at around this value. The wavemaker non-linearity order is chosen as 2 and the HOS order is chosen as 5 after additional sensitivity studies, including calculations made using HOS order of 8, show the convergence of the solution at these selected values. Most of the non-linearity is taken into account with mHOS = 5, and it is the typical choice for practical applications of large-scale ocean sea states as it provides the best compromise in accuracy and CPU time balance.

# 5. Results and Discussion

Figure 7 shows still images of video captured for the wave EF6-54 and $\epsilon_n$ = 0.0051 at $t$ = 49.80, 49.90, 50.16, 50.24, 50.32, 50.46, 50.60 and 51.44 s. Wave breaking onset time is calculated by HOS-NWT at 49.80 s and estimated as 49.90 s during experiments. At these instances, the vertical white line indicates the measurement for the location of the wave crest. The location of the horizontal ruler is known with respect to the mean position of the wavemaker, and therefore we can extract the position of the wave crest at times 49.80 s and 49.90 s. At $t$
513 = 50.16 s, the wave has just began to overturn and is collapsing at t = 50.24 s
514 and 50.32 s. Evidence of splash up and whitewater appears at t = 50.46 s and
515 50.60 s, and all of the wave-breaking energy loss appears to be complete at t =
516 51.44 s.

517 Figure 8 plots the location of wave breaking onset ($x_b$) calculated by HOS-NWT versus the location measured during experiments, and Fig. 9 plots the time of wavebreaking onset ($t_b$) calculated by HOS-NWT and measured during the LHEEA experiments. Measurements and calculations for $x_b$ and $t_b$ are also given in Tables 4, 5 and 6 for breaking waves generated using energy focusing, modulational instability and irregular waves, respectively. As there were over 50 breaking waves total in the modulational instability waves, only the first 10 breaking waves from each case are listed in Table 5. Dashed lines present under the HOS-NWT column indicate there was no breaking wave calculated close to that time and location. For many of the wave cases, there were additional breaking waves in both the experiments and calculations that are not presented. Only waves where measurements were made for breaking onset and location during experiments are presented along with the corresponding calculations. The type of breaker is also presented in Tables 4 - 6, and is determined during the experiments based on visual observation.

532 Breaking time and location are predicted accurately by HOS-NWT when compared with experimental measurements. Breaking onset location and time are measured during experiments by visual observation of a video recording, therefore the determination of breaking onset is somewhat subjective, providing a good explanation for differences between measurements and calculations. However, the agreement is still excellent for each of the wave breaker types and for the different types of wave breaking generation. This study confirms the wave breaking parameter $B_x$ provides an accurate means of predicting breaking onset regardless of breaker type, strength, or scale.

541 It is interesting to note, based on visual observation and calculations, wave-breaking onset occurs well before the appearance of white capping, on the order of half a wave period. This is evident by the breaking wave presented in Fig. 7
<figure><img src="image_8.png" /><figcaption>Figure 7: Still images of video captured for the wave EF6-54 and $\epsilon_n = 0.0051$ at $t = 49.80, 49.90, 50.16, 50.24, 50.32, 50.46, 50.60$ and $51.44$ s. Wave breaking onset is estimated at $t_{br} = 49.90$ s and at this time, a measurement of the location of the crest with respect to the measuring stick is made, indicated by the vertical white line, and is used to calculate the distance to the mean position of the wavemaker.</figcaption></figure>
544 where breaking onset occurs at $t = 49.90$ s but white-capping is not observed until after $t = 50.46$ s.

<figure><img src="image_9.png" /><figcaption>Figure 8: Location of wavebreaking onset ($x_b$) measured during experiments is compared with calculations using HOS-NWT.</figcaption></figure>

545

546

## 6. Conclusions

547

548

549

550

551

Breaking waves are studied experimentally and numerically using the non-linear, potential flow solver HOS-NWT. Single and multiple breaking waves are generated by energy focusing and modulational instability techniques, as well as by steep irregular wave propagation and irregular waves with an energy focusing wave superimposed.

552

553

554

555

556

557

558

559

A wave breaking onset criteria is implemented in HOS-NWT which calculates wave breaking when the ratio of water particle velocity at the surface to wave crest velocity exceeds a value of 0.85. To calculate the onset of a breaking wave at some random location in the domain, the breaking parameter is calculated at every point in the domain and at every time step. To achieve this, the instantaneous local phase velocity is calculated at every point in the domain and at every time step by use of the partial Hilbert transform and dispersion relation. A comparison of the local phase velocity with more traditional crest
<figure><img src="image_10.png" /><figcaption>Figure 9: Time of wavebreaking onset ($t_b$) measured during experiments is compared with calculations using HOS-NWT.</figcaption></figure>

tracking methods shows that at the peak of a wave about to break, the local phase velocity provides a good approximation of the crest velocity, and the prediction of wave breaking using this method achieves excellent results.

The wave-breaking onset criteria implemented in the model is successful at predicting breaking events as demonstrated by the excellent agreement between calculations and experimental measurements for wave-breaking onset time and location. The success of the model when predicting breaking waves for a wide range of waves tested indicates the wave-breaking onset criteria is robust for different types and strengths of breakers and at different scales, and could potentially be applied in other nonlinear potential flow solvers.

Breaking waves are highly nonlinear and rapid processes, making the prediction of onset challenging. Difference between the time and location of breaking events in calculations and measurements can be explained by the nature of observing breaking time and location during experiments and the somewhat subjective nature of the definition of wave-breaking onset when making visual observations. Now that a robust method of predicting breaking onset in HOS-NWT is validated, the next step is to determine a means to calculate wave-
583 Acknowledgements

<footnote n="577">577 breaking associated energy losses. It is left for future work to determine success
578 of the breaking onset parameter for varying depth as well as in a multidirec-
579 tional wave field. By calculating wave-breaking onset and energy dissipation in
580 HOS-NWT, not only will the solver be able to accurately calculate more ener-
581 getic sea states, but the potential for numerical instabilities caused by breaking
582 waves that can lead to failure of the model is greatly reduced.</footnote>

<footnote n="584">584 The authors are grateful for the experimental data provided by Dr. Zhigang
585 (Kevin) Tian.</footnote>
[1] X. Barthelemy, M. L. Banner, W. L. Peirson, F. Fedele, M. Al-lis, F. Dias, On a unified breaking onset threshold for gravity waves in deep and intermediate depth water, submitted. URL: https://arxiv.org/pdf/1508.06002.pdf.

[2] A. V. Babanin, D. Chalikov, I. R. Young, I. Savelyev, Numerical and laboratory investigation of breaking of steep two-dimensional waves in deep water, Journal of Fluid Mechanics 644 (2010) 433.

[3] D. Chalikov, A. V. Babanin, Simulation of Wave Breaking in One-Dimensional Spectral Environment, Journal of Physical Oceanography (2012).

[4] J. Grue, A. Jensen, Orbital velocity and breaking in steep random gravity waves, Journal of Geophysical Research: Oceans 117 (2012) 1–16.

[5] M. Perlin, W. Choi, Z. Tian, Breaking Waves in Deep and Intermediate Waters, Annual Review of Fluid Mechanics 45 (2013) 115–145.

[6] A. V. Babanin, T. Waseda, T. Kinoshita, a. Toffoli, Wave Breaking in Directional Fields, Journal of Physical Oceanography 41 (2011) 145–156.

[7] A. Babanin, D. Chalikov, I. Young, I. Savelyev, Predicting the breaking onset of surface water waves, Geophysical Research Letters 34 (2007).

[8] R. J. Rapp, W. K. Melville, Laboratory measurements of deep-water breaking waves 331 (1990) 735–800.

[9] M. P. Tulin, T. Waseda, Laboratory observations of wave group evolution, including breaking effects, Journal of Fluid Mechanics 378 (1999) 197–232.

[10] Z. Tian, M. Perlin, W. Choi, Evaluation of a deep-water wave breaking criterion, Physics of Fluids 20 (2008) 066604.
611

612

613

[11] Z. Tian, M. Perlin, W. Choi, An eddy viscosity model for two-dimensional breaking waves and its validation with laboratory experiments, Physics of Fluids 24 (2012).

614

615

616

[12] A. Saket, W. Peirson, M. Banner, X. Barthelemy, M. Allis, On the threshold for wave breaking of two-dimensional deep water wave groups in the absence and presence of wind., Journal of Fluid Mechanics 811 (2017) 642–658.

617

618

619

[13] P. Stansell, C. MacFarlane, An Experimental Investigation of Breaking Criteria Based on Wave Phase Speeds, Journal Of Physical Oceanography 32 (2002) 1269–1283.

620

621

[14] C. H. Wu, H. M. Nepf, Breaking criteria and energy losses for three-dimensional wave breaking, Journal of Geophysical Research 107 (2002).

622

623

[15] M. P. Tulin, M. Landrini, Breaking waves in the ocean and around ships, in: Twenty-Third Symposium on Naval Hydrodynamics, 2001, pp. 713–745.

624

625

[16] W. K. Melville, P. Matusov, Distribution of breaking waves at the ocean surface., Nature 417 (2002) 58–63.

626

627

628

[17] M. L. Banner, W. L. Peirson, Wave breaking onset and strength for two-dimensional deep-water wave groups, Journal of Fluid Mechanics 585 (2007) 93–115.

629

630

631

[18] Z. Tian, M. Perlin, W. Choi, Energy dissipation in two-dimensional unsteady plunging breakers and an eddy viscosity model, Journal of Fluid Mechanics 655 (2010) 217–257.

632

633

[19] L. Shemer, On kinematics of very steep waves, Natural Hazards and Earth System Science 13 (2013) 2101–2107.

634

635

[20] R. Kurnia, E. van Groesen, High order Hamiltonian water wave models with wave-breaking mechanism, Coastal Engineering 93 (2014) 55–70.
636

637

638

[21] J. Song, M. L. Banner, On Determining the Onset and Strength of Breaking for Deep Water Waves. Part I: Unforced Irrotational Wave Groups, Journal of Physical Oceanography 32 (2002) 2541–2558.

639

640

641

[22] M. Derakhti, J. Kirby, Breaking-onset, energy and momentum flux in unsteady focused wave packets, Journal of Fluid Mechanics 790 (2016) 553–581.

642

643

644

[23] M. Derakhti, J. T. Kirby, F. Shi, G. Ma, Wave breaking in the surf zone and deep-water in a non-hydrostatic RANS model. part 1: Organized wave motions, Ocean Modelling (2016).

645

646

647

[24] M. Derakhti, J. T. Kirby, F. Shi, G. Ma, Wave breaking in the surf zone and deep-water in a non-hydrostatic RANS model. part 2: Turbulence and mean circulation, Ocean Modelling (2016).

648

649

650

[25] M. L. Banner, X. Barthelemy, F. Fedele, M. Allis, A. Benetazzo, F. Dias, W. L. Peirson, Linking Reduced Breaking Crest Speeds to Unsteady Non-linear Water Wave Group Behavior, Physical Review Letters 112 (2014).

651

652

653

[26] S. T. Grilli, J. Skourup, I. A. Svendsen, An efficient boundary element method for nonlinear water waves, Engineering Analysis with Boundary Elements 6 (1989) 97 – 107.

654

655

656

[27] S. T. Grilli, P. Guyenne, F. Dias, A fully non-linear model for three-dimensional overturning waves over an arbitrary bottom, Int. J. Numer. Meth. Fluids 35 (2001) 829–867.

657

658

[28] M. Banner, D. H. Peregrine, Wave breaking in deep water, Annual Review of Fluid Mechanics 25 (1993) 373–397.

659

660

661

[29] G. Ducrozet, F. Bonnefoy, D. Le Touzé, P. Ferrant, 3-D HOS simulations of extreme waves in open seas, Natural Hazards and Earth System Science 7 (2007) 109–122.
[30] F. Bonnefoy, G. Ducrozet, D. Le Touzé, P. Ferrant, Advances in Numerical Simulation of Nonlinear Water Waves, volume 11 of *Advances in Coastal and Ocean Engineering*, World Scientific, 2009, pp. 129–164.

[31] G. Ducrozet, F. Bonnefoy, D. Le Touzé, P. Ferrant, A modified High-Order Spectral method for wavemaker modeling in a numerical wave tank, European Journal of Mechanics, B/Fluids 34 (2012) 19–34.

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

[33] G. Ducrozet, F. Bonnefoy, D. Le Touzé, P. Ferrant, HOS-ocean: Open-source solver for nonlinear waves in open ocean based on high-order spectral method, Computer Physics Communications (2016).

[34] B. R. Seiffert, G. Ducrozet, Deep water wave-breaking in a high-order spectral model, in: 31st International Workshop on Water Waves and Floating Bodies, Plymouth, MI, 2016.

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

[36] B. J. West, K. A. Brueckner, R. S. Janda, D. M. Milder, R. L. Milton, A new numerical method for surface hydrodynamics, J. Geophys. Res 92 (1987) 803–11.

[37] G. Ducrozet, F. Bonnefoy, P. Ferrant, On the equivalence of unidirectional rogue waves detected in periodic simulations and reproduced in numerical wave tanks, Ocean Engineering 117 (2016) 346–358.

[38] M. Gouin, G. Ducrozet, P. Ferrant, Development and validation of a non-linear spectral model for water waves over variable depth, European Journal of Mechanics-B/Fluids 57 (2016) 115–128.
[39] V. Zakharov, Stability of periodic waves of finite amplitude on the surface of a deep fluid, J. App. Mech. Tech. Phys. 9 (1968) 86-94.

[40] Y. Agnon, H. B. Bingham, A non-periodic spectral method with application to nonlinear water waves, European Journal of Mechanics-B/Fluids 18 (1999) 527-534.

[41] V. Sriram, T. Schlurmann, S. Schimmels, Focused wave evolution using linear and second order wavemaker theory, Applied Ocean Research 53 (2015) 279-296.

[42] F. Fedele, Geometric phases of water waves, EPL (Europhysics Letters) 107 (2014) 69001.

[43] F. Fedele, C. Chandre, M. Farazmand, Kinematics of fluid particles on the sea surface: Hamiltonian theory., Journal of Fluid Mechanics 801 (2016) 260-288.