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

HAL
open science

# Breaking waves and their applications to the offshore wind sector

Jean-François Filipot

► To cite this version:

Jean-François Filipot. Breaking waves and their applications to the offshore wind sector: Habilitation to drive research (Habilitation à Diriger des Recherches, HDR). Engineering Sciences [physics]. Université de Bretagne Occidentale, 2022. tel-04184079

HAL Id: tel-04184079

https://hal.science/tel-04184079

Submitted on 21 Aug 2023

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.

Public Domain
Mémoire de recherche présenté par

Jean-François Filipot

Pour l'obtention du diplôme
d'Habilitation à Diriger des Recherches
de l'Université de Bretagne Occidentale

# Les déferlantes et leurs applications au domaine
de l'éolien en mer

## Breaking waves and their applications to the
offshore wind sector

Soutenu le 13 Mai 2022
devant le jury composé de

### Rapporteurs

<table><tr><td>Bruno Castelle</td><td>Professeur, Université de Bordeaux</td></tr><tr><td>Xiaoli Guo Larsen</td><td>Scientifique Senior, DTU Wind Energy</td></tr><tr><td>Lili Kimmoun</td><td>Maîtresse de conférence HDR, Centrale Marseille</td></tr></table>

### Examinateurs/trices

<table><tr><td>France Floc'h</td><td>Maîtresse, Université de Bretagne Occidentale</td></tr><tr><td>Marc Prevosto</td><td>Chercheur, indépendant</td></tr></table>

### Directrice de recherche

<table><tr><td>Marie-Noëlle Bouin</td><td>Chercheuse (HDR), Météo France</td></tr></table>
Acknowledgements

I am very grateful to many people, since I have inspired by different mentors along
my carrier. First I would like to thank X. Carton for the confidence he provided
me, showing that I could access any of the top colleges in France. Of course, I
am indebted to F. Ardhuin who taught me almost everything I know about wave
physics, which is probably less that 10% of what he actually knows. Pr. Cheung,
my advisor at U. Hawaii, taught me how to produce efficient and good research
material, defining exactly what was the expected results from the beginning.
B. Chapron, beyond his brilliant and often crazy ideas, taught me that it was
important to keep having fun while working, a philosophy that I attemp to apply
on a daily basis. I would like to also thank M. Prevosto, who taught me wisdom
and pragmatism along our different common projects. My discussions with my
friend M. Delpey have given me faith that it was worth continuing research even
when we are overwhelmed by extra-research duties. In this regards, I would like
to thank all my FEM colleagues who have always encouraged me to pursue my
research, although they have had to bear some my other duties to do so. I would
also like to thank M.N. Bouin who accepted to act as my research director and
for her guidance in the writing up of the manuscript. My colleagues F. Leckler
and R. Duarte, and different PhD students (A. Varing, P. Platzer, F. Hulin) and
postdocs (A. Ruju, P. Guimaraes, C. Stringari, P. Renaud, M. Battle) are major
contributors to my research activities. I am also indebted to all the entities from
the private and public sectors which invested in the research activities we are
carrying at FEM about the interactions between breaking waves and offshore
wind turbines.

Enfin, je voudrais aussi dire á ma femme Emmanuelle et mes enfants Côme
et Gauthier que si mon regard est souvent perdu dans le(s) vague(s) mon coeur
leur sera toujours réservé.
# Contents

<table><tr><td>1</td><td>Introduction</td><td>11</td></tr><tr><td>2</td><td>The breaking onset criteria</td><td>15</td></tr><tr><td>2.1</td><td>Introduction</td><td>15</td></tr><tr><td>2.2</td><td>The geometric breaking criteria</td><td>15</td></tr><tr><td>2.2.1</td><td>Wave steepness approaches</td><td>15</td></tr><tr><td>2.2.2</td><td>Wave front local steepness approaches</td><td>16</td></tr><tr><td>2.3</td><td>The dynamic breaking criteria</td><td>17</td></tr><tr><td>2.4</td><td>The kinematic breaking criterion</td><td>18</td></tr><tr><td>2.5</td><td>Discussions and conclusions</td><td>19</td></tr><tr><td colspan="2">Appendices</td><td>21</td></tr><tr><td>2.A</td><td>Paper from Varing & al. (2021a)</td><td>21</td></tr><tr><td>3</td><td>Breaking statistics models and observations</td><td>37</td></tr><tr><td>3.1</td><td>Introduction</td><td>37</td></tr><tr><td>3.2</td><td>Breaking wave probability parameterizations</td><td>37</td></tr><tr><td>3.2.1</td><td>Breaking probability parameterizations based on the Gaussian field theory</td><td>37</td></tr><tr><td>3.2.2</td><td>Observation-based breaking probability parameterizations</td><td>38</td></tr><tr><td>3.3</td><td>Parameterization for the spectral density of breaking crests</td><td>40</td></tr><tr><td>3.4</td><td>A new wave breaking parameterization</td><td>41</td></tr><tr><td>3.4.1</td><td>Main underlying assumptions</td><td>41</td></tr><tr><td>3.4.2</td><td>Introduction of a linear-equivalent breaking onset criterion</td><td>43</td></tr><tr><td>3.5</td><td>Observations of breaking waves statistics</td><td>43</td></tr><tr><td>3.6</td><td>Discussions & Conclusions</td><td>44</td></tr><tr><td colspan="2">Appendices</td><td>47</td></tr><tr><td>3.A</td><td>Paper from Stringari & al. (2021b)</td><td>47</td></tr><tr><td>3.B</td><td>Paper from Stringari & al. (2021a)</td><td>67</td></tr><tr><td>4</td><td>Spectral models and coastal observations</td><td>81</td></tr><tr><td>4.1</td><td>Wave breaking source terms</td><td>81</td></tr></table>
<table>
  <tbody>
    <tr>
      <td>4.1.1</td>
      <td>The spectral wave models</td>
      <td>81</td>
    </tr>
    <tr>
      <td>4.2</td>
      <td>Developments in WW3 &amp; SWAN</td>
      <td>83</td>
    </tr>
    <tr>
      <td>4.2.1</td>
      <td>Development of the WW3 model</td>
      <td>83</td>
    </tr>
    <tr>
      <td>4.2.2</td>
      <td>Evaluation of WW3 source terms accuracy in extreme conditions</td>
      <td>84</td>
    </tr>
    <tr>
      <td>4.2.3</td>
      <td>Development of the SWAN model</td>
      <td>85</td>
    </tr>
    <tr>
      <td>4.3</td>
      <td>Model validation in coastal areas</td>
      <td>87</td>
    </tr>
    <tr>
      <td>4.4</td>
      <td>Discussions &amp; Conclusion</td>
      <td>88</td>
    </tr>
    <tr>
      <td colspan="2"><strong>Appendices</strong></td>
      <td><strong>89</strong></td>
    </tr>
    <tr>
      <td>4.A</td>
      <td>Paper from Ruju &amp; al. (2020)</td>
      <td>89</td>
    </tr>
    <tr>
      <td>4.B</td>
      <td>Paper from Filipot &amp; Cheung (2012)</td>
      <td>109</td>
    </tr>
    <tr>
      <td>4.C</td>
      <td>Paper from Varing &amp; al. (2021b)</td>
      <td>123</td>
    </tr>
    <tr>
      <td><strong>5</strong></td>
      <td><strong>Breaking waves effects on offshore structures</strong></td>
      <td><strong>141</strong></td>
    </tr>
    <tr>
      <td>5.1</td>
      <td>Introduction</td>
      <td>141</td>
    </tr>
    <tr>
      <td>5.2</td>
      <td>Past slamming loads field experiments</td>
      <td>143</td>
    </tr>
    <tr>
      <td>5.3</td>
      <td>Slamming loads on offshore structures</td>
      <td>143</td>
    </tr>
    <tr>
      <td>5.4</td>
      <td>Database of the geometry and kinematics of breaking waves</td>
      <td>146</td>
    </tr>
    <tr>
      <td>5.5</td>
      <td>Discussions &amp; Conclusion</td>
      <td>148</td>
    </tr>
    <tr>
      <td colspan="2"><strong>Appendices</strong></td>
      <td><strong>149</strong></td>
    </tr>
    <tr>
      <td>5.A</td>
      <td>Paper from Filipot &amp; al. (2019)</td>
      <td>149</td>
    </tr>
    <tr>
      <td>5.B</td>
      <td>Paper from Guimaraes &amp; al. (2020)</td>
      <td>171</td>
    </tr>
    <tr>
      <td><strong>6</strong></td>
      <td><strong>Perspectives</strong></td>
      <td><strong>185</strong></td>
    </tr>
    <tr>
      <td>6.1</td>
      <td>Modeling waves in extreme conditions</td>
      <td>185</td>
    </tr>
    <tr>
      <td>6.1.1</td>
      <td>Deep and intermediate water wave breaking dissipation parameterization</td>
      <td>185</td>
    </tr>
    <tr>
      <td>6.1.2</td>
      <td>Shallow water wave breaking dissipation parameterization</td>
      <td>186</td>
    </tr>
    <tr>
      <td>6.1.3</td>
      <td>Spatial variability</td>
      <td>187</td>
    </tr>
    <tr>
      <td>6.2</td>
      <td>Estimation of slamming loads on offshore wind turbines</td>
      <td>189</td>
    </tr>
    <tr>
      <td>6.2.1</td>
      <td>Determining the severity of breaking waves</td>
      <td>189</td>
    </tr>
    <tr>
      <td>6.2.2</td>
      <td>Linear-Equivalent breaking criterion and severity</td>
      <td>190</td>
    </tr>
    <tr>
      <td>6.3</td>
      <td>Application to wind characterization</td>
      <td>191</td>
    </tr>
    <tr>
      <td>6.3.1</td>
      <td>Wave breaking and wind retrieval from space</td>
      <td>192</td>
    </tr>
    <tr>
      <td>6.3.2</td>
      <td>Air-sea coupled modeling</td>
      <td>195</td>
    </tr>
    <tr>
      <td>6.3.2.1</td>
      <td>Effect of breaking waves on the momentum flux at the air-sea interface through air flow separation</td>
      <td>195</td>
    </tr>
    <tr>
      <td>6.3.2.2</td>
      <td>Effect of breaking waves on the momentum and heat flux at the air-sea interface through sea spray generation</td>
      <td>195</td>
    </tr>
    <tr>
      <td>6.4</td>
      <td>Development of real scale laboratories</td>
      <td>196</td>
    </tr>
  </tbody>
</table>
<table><tr><td>CONTENTS</td><td>5</td></tr><tr><td>6.4.1 Perpetuation of La Jument laboratory</td><td>197</td></tr><tr><td>6.4.2 Implementation of a monitoring system on a real scale offshore wind turbine</td><td>197</td></tr><tr><td>7 Conclusion</td><td>199</td></tr><tr><td>A Curriculum Vitae</td><td>215</td></tr></table>
# Glossary

**BOSZ** Boussinesq model for Ocean and Surf Zone. 18

**CAPEX** Capital expenditures. 191

**CARAVELE** FEM research project, addressing the characterization of the wind resource for the offshore wind sector, with a focus on extreme conditions
More references here: https://www.france-energies-marines.org/en/projects/caravele/. 193, 195

**CASSIOWPE** FEM research project, addressing the air-sea coupled modeling for the characterization of the wind resource in the Gulf of Lion, France.
More references here: https://www.france-energies-marines.org/en/projects/cassiowpe/. 185, 196

**Cerema** Centre d'études et d'expertise sur les risques, l'environnement, la mo-bilité et l'aménagement. 144

**DiMe** FEM research project dealing with the characterization of breaking waves with applications to the design of Ocean Renewable Energy converters.
More references here: https://www.france-energies-marines.org/en/projects/dime/. 12, 18, 37, 41, 84, 185, 189, 193, 199

**DIMPACT** FEM research project, addressing slamming loads of floating off-shore wind turbines. More references here: https://www.france-energies-marines.org/en/projects/dimpact/. 12, 43, 143, 148, 186, 190, 191, 197, 199

**FEM** France Energies Marines. 12, 84, 144, 185, 188, 190, 193, 196, 199

**HZG** Helmholtz Zentrum Geesthacht. 144, 145

**Ifremer** Institut Francais de Recherche pour l'Exploitation de la Mer. 144, 190

**IRPHE** Institut de Recherche sur les Phénomènes Hors Equilibre. 187

**LCoE** Levelized Cost of Energy. 11, 12, 191
**LOPS** Laboratoire d'Océanographie Physique et Spatiale. 196

**MARLIT** Project dealing with the means of assessing and mitigating local coastal risks due to storm surges.. 187, 188

**NRCS** Normalized Radar Cross Section. 192

**RSMAS** Rosenstiel School of Marine and Atmospheric Science, U. Miami, Florida, USA. 193

**SAR** Synthetic Aperture Radar. 87, 187, 189

**Shom** Service Hydrographique et Océanographique de la Marine. 144

**SUSTAIN** SURGE-STRUCTURE-ATMOSPHERE INTERACTION FACILITY, U. Miami, Florida, USA. 193

**SWAN** Simulating WAves Nearshore. 4, 82, 85-87

**ULS** Ultimate Limit State. 37, 185, 189, 191

**WW3** WAVEWATCHIII. 83, 84, 185, 187
Abstract

This document presents research efforts dedicated to improve the characterization
of breaking waves and of their multiple effects on offshore wind turbines.

The manuscript focuses first on the definition of the breaking onset criterion
that is important for the design of offshore wind turbines. It further presents
recent advances in the parameterization of the wave breaking statistics, needed
to identify the design sea states, when slamming loads are considered. This is
followed by a description of efforts to improve the spectral wave models, used to
obtain the wave climatology at offshore wind farms location.

Observational works to capture area of wave energy focusing (potentially fea-
turing breaking waves) are then introduced. The document further presents ex-
perimental efforts to observe interactions between large amplitude breaking waves
and offshore structures.

Finally, different perspectives targeting the improvement of the above-cited
research axes are given, and extension to other applications are discussed. These
applications cover the effects of breaking waves on the characterization of the wind
properties with spaceborne measurements and ocean-wave-atmosphere coupled
numerical models.
# Chapter 1

# Introduction

This manuscript presents different contributions to the wave breaking research field with its application to different offshore wind sector challenges. The corresponding research work primarily intends to reduce the uncertainties in the design of the turbines but carries further implications in the assessment of the wind resource.

The offshore wind industry is growing quickly, stimulated by the need for greener energy production and the continuous reduction of its levelized cost of energy (LCOE). As an illustration, for the first time in 2010, the global capacity of offshore wind exceeded 1 gigawatt (GW) while in 2018, 4.3 GW of new offshore wind capacity was installed. The projections given by International Energy Agency Wind foreshadow an increase by at least 15-fold by 2040 at global scale. This growth scenario implies a series of improvements to make the offshore wind energy production even more reliable and profitable.

Offshore wind turbines are deployed in areas of high winds that often generate severe sea states. These sea states carry breaking waves of different scales that have a wide spectrum of implications (sometime surprising) on the offshore wind sector.

To further understand their impact, let us enumerate some of the breaking waves key characteristics :

1. Breaking waves may behave as walls of water (the horizontal fluid velocity at crest equaling the crest phase speed)

2. They are the main sink of wave energy, through generation of turbulence and work against the bubbles buoyancy

3. They induce currents that may dominate the nearshore hydrodynamics

4. Their sharp crests induced the so-called "air-flow" separation that modifies the wind drag over the sea surface
5. They constitute a source of sea spray that contributes to heat and momentum fluxes at the air-sea interface

6. When observed from space, their white foam carries a distinctive signature that alters electromagnetic signal backscatter.

From these statements, it follows that:

a. breaking waves may induce large impulsive loads on offshore wind turbines, also known as slamming loads, that can be the design loads (from 1.)

b. the breaking-induced dissipation needs to be correctly parameterized in numerical models to simulate the wave conditions at an offshore wind farm (from 2.)

c. breaking waves should be properly captured in the nearshore area to understand and model the sediment transport at bottom-fixed offshore wind farm and the morphodynamics of the power cable landing area (from 3.)

d. breaking waves effect should be taken into account in the modeling of the atmospheric and wave conditions especially in storm/hurricane conditions where air-sea interactions are a dominant process (from 4. and 5.)

e. breaking waves should be accounted for wind retrieval from satellite-borne measurements (from 6.)

Hence, increasing the knowledge about breaking waves may contribute to refine load estimations on the turbines, reduce uncertainties on extreme wind and wave conditions and to quantify the sediment transport at bottom-fixed offshore wind farms (and at power cables landing zones). In fine, this should lead to more reliable offshore wind turbines, a lower LCoE and possibly an increased acceptability of the offshore wind activities. These are the main sources of motivation of the present work.

Despite contributions from some of the greatest scientific minds, such as the 2021 Physics Nobel Prize Klaus Hasselmann, wave breaking remains a poorly understood process and all items listed above still represent open challenges with a great potential for improvements.

This document recalls research activities conducted, in particular, in the frame of the DiMe and DIMPACT projects in strong collaboration with several PhD students, postdoctoral candidates and permanent scientists from FEM (France Energies Marines) and other institutes. It addresses a series of challenges in the field of the definition, statistics, observations and modeling of breaking waves, namely:

* the definition of the breaking onset criterion (Chapter 2)
* the derivation of new wave breaking statistics parameterizations and related observational efforts (Chapter 3)

* the development of spectral wave models and observational means for their validation (Chapter 4)

* the observation of breaking waves interactions with offshore structures (Chapter 5)

Perspectives are finally given to present future research works either related to the above-cited topics or new challenges addressing the impact of breaking waves on wind retrieval from spaceborne instruments and on air-sea interactions.
# Chapter 2

# The breaking onset criteria

## 2.1 Introduction

Being able to determine when and where breaking onset occurs is critical for the Ultimate Limit State (ULS) design of offshore wind turbines to correctly simulate slamming loads on the foundations with numerical wave models. It further directly impacts the definition of appropriate wave breaking statistics parameterizations, as discussed in Chapter 3.

Breaking onset parameters are generally classed into three different classes: the geometric, dynamic and kinematic criteria. The present chapter will recall the main research efforts that have been dedicated to each of these classes with a focus on the work of Filipot et al. [2010] and Varing et al. [2021b] for the geometric and kinematic criterion respectively.

## 2.2 The geometric breaking criteria

### 2.2.1 Wave steepness approaches

Miche [1944] demonstrated that the limit form of regular waves over flat bottom is reached when:

$$ \frac{kH_{max}}{0.88 \tanh(kd)} = 1, \qquad (2.1) $$

with $k$, $H_{max}$ and $d$ the wavenumber, limit wave height and water depth, respectively. On one hand, this result extends the limit wave steepness $ak = 0.443$, with $a$ the wave amplitude, derived by Stokes [1880] and in the other hand, it reduces to $H/d = \gamma$ in shallow water. The so-called breaker index, $\gamma$, is equals to 0.88 under Miche's waves assumption but generally depends on the beach slope [Allsop et al., 1999], the offshore wave steepness [Battjes and Stive, 1985], the non-dimensional water depth [Ruessink et al., 2003] among other parameters.
<figure><img src="image_2.png" /></figure>

Figure 2.1: From Filipot et al. [2010]. Evolution of Miche’s limit criterion as a function of $X = \tanh(kd)$. The black squares show the analytical solution for a regular nonlinear wave at breaking onset over a flat bottom, the red crosses give the value of the same parameter but estimated using a linear wave of same potential energy, the blue squares correspond to a polynomial fit derived by Filipot et al. [2010] and the black crosses are a fit proposed by Ruessink et al. [2003] for shallow water applications.

Filipot et al. [2010] revisited Miche [1944] work to propose a breaking onset criteria valid from deep to shallow water and suitable for the development of wave breaking dissipation source terms in spectral models. Indeed, in such numerical wave models, the wave nonlinearities that modify the wave shape prior to wave breaking cannot be captured. Filipot et al. [2010] proposed a way around by computing a linear-equivalent breaking criteria, assessing the properties of a linear wave with the same potential energy as Miche’s limit wave (see figure 2.1).

### 2.2.2 Wave front local steepness approaches

Because of the observed variability of the Miche parameter at breaking onset for realistic irregular 2D or 3D waves, [Tian et al., 2008, 2012], other scientists explored the wave front local steepness as another potential breaking onset criterion [Bonmarin, 1989; Rapp and Melville, 1990]. Still, even with this approach, the quest for an universal steepness-based breaking parameter appears hopeless: for instance Bonmarin [1989] found that the front face steepness is highly dependent on the type of breaking waves, varying from 0.38 to 0.61 for 2D spilling and plunging breakers, respectively. Three-dimensional effects further disqualify front face steepness breaking onset criteria as reported by She et al. [1994]: 3D waves are able to reach significantly higher front face steepness that equivalent 2D waves. Finally, the way breaking waves are generated (bottom effects, Ben-
jamin Feir instability or wave packet focusing) and the wind effect may also affect the front face steepness at breaking onset as shown by the results of *Phillips and Banner* [1974]; *Touboul et al.* [2006]; *Kharif et al.* [2008]; *Diorio et al.* [2009].

## 2.3 The dynamic breaking criteria

There are two main families of dynamic breaking criteria : those based on acceleration considerations and those relying on energetic considerations. The first family arose from the pioneering work of *Longuet-Higgins* [1963] who demonstrated that a limit Stokes (or Miche) regular deep water wave has a particle acceleration at the crest of 0.5 g. This finding was later exploited by *Snyder and Kennedy* [1983]; *Kennedy and Snyder* [1983]; *Snyder et al.* [1983] to provide a deep water wave breaking probability parameterization. Still, *Ochi and Tsai* [1983] and *Ramberg and Griffin* [1987] reported that for real irregular wave, the value of the crest acceleration at breaking onset was closer to 0.4 while *Dawson et al.* [1993] found an optimal value of 0.335. As for the limit wave steepness, such scatter in the breaking threshold questions the generalization of acceleration-based breaking onset parameters.

The second family of dynamic breaking criteria investigates the energy flux within wave packets, as a potential indicator for the occurrence of breaking. *Song and Banner* [2002] reported from numerical investigations of 2D focused wave packets that the mean energy convergence rate at the center of wave groups determines the evolution towards breaking. These findings were supported by studies from *Banner and Peirson* [2007a]; *Tian et al.* [2008]. Later, the same concept of energy flux was applied to the crest region by *Perlin et al.* [2013]; *Derakhti et al.* [2016]; *Saket et al.* [2017, 2018], providing a more robust and tractable framework.

Surfing on the same approach, the work by *Barthelemy et al.* [2018] has recently become a reference in the wave breaking community. Arguing, again that breaking onset occurs due to an excess of energy flux in steep wave crests, they introduce the following dynamical breaking onset parameter:

$$ \mathbf{B} = \mathbf{F}(E \| c \|) \qquad (2.2) $$

With E the energy density, c the wave phase speed and F the local energy flux. When evaluated at the sea surface, B reduces to $B_x = u_c/c$, with $u_c$ the horizontal fluid velocity at the crest. $B_x$ is hence equivalent to a kinematic criterion, as discussed in the next section. Investigating nonlinear numerical simulations of focused wave packets in intermediate and deep water, *Barthelemy et al.* [2018] main finding is that a no-return point, (after which a wave will inevitably and shortly evolve toward breaking) is reached for $u_c/c \simeq 0.85$. The physical reason why this is observed remains unclear, but since then, this criterion has been
applied to capture breaking onset in nonlinear wave models [e.g. Seiffert and Ducrozet, 2018; Simon et al., 2019].

## 2.4 The kinematic breaking criterion

The kinematic breaking parameter, $u_c/c$, is currently gaining a renewed interest recently with the results from Barthelemy et al. [2018] presented above. Still, the value of the threshold (0.85) identified by these authors, that departs from one, continues to trigger the interest of the wave breaking community.

Experimental investigations also suggest that $u_c/c$ at breaking onset shows an appreciable and unexplained scatter. Using Particle Image Velocimetry (PIV) measurements Perlin et al. [1996] found $u_c/c = 0.74$ just prior to breaking for deep water plunging breakers while Chang and Liu [1998] reported $u_c/c = 0.87$ studying shallow water monochromatic breaking waves. Qiao and Duncan [2001], reported $u_c/c$ in the range [0.75, 0.95] for spilling breakers, sensed with PIV. Similar orders of magnitude were further reported by Melville and Matusov [2002]; Banner and Peirson [2007b] and Tian et al. [2010].

This apparent observational scatter of $u_c/c$ at breaking onset, combined with the powerful but intriguing findings from Barthelemy et al. [2018] motivated, during the DiMe project, a research effort detailed in Varing et al. [2021b] (see Appendix 2.A).

The initial goal of the latter authors was to verify the validity of Barthelemy et al. [2018]'s work for shallow water breaking waves, with the final objective to implement the corresponding criterion in the BOSZ model (Boussinesq model for Ocean and Surf Zone) [Roeber et al., 2010]. This study yielded an additional and original result: the definition of a new kinematic breaking onset criterion that offers a framework to reconcile theory and lab observations. Varing et al. [2021b] exploited a database of both breaking and nearly breaking waves generated by the Fully Nonlinear Potential Flow Model of Grilli et al. [1989] and Grilli and Subramanya [1996]. This kind of approach provides reliable information on the wave geometry and fluid velocity within the crests until breaking onset. For practical reasons, the study focuses on solitary and pseudo-regular waves propagating in shallow water but the results should carry potential for a generalization to any kind of waves. For the wave families studied, we showed that using $||u_m||$, the magnitude of the maximum surface fluid velocity in the front face region, instead of $u_c$, reduces the scatter of the breaking kinematic criterion. Intuitively this finding translates the fact that the fluid instability and associated jet that initiates breaking may occur slightly ahead of the tip of the crest. Using $||u_m||/c$, instead of $u_c/c$, further yields values closer to unity (the "theoretical" value).

Varing et al. [2021b]'s results further supported the validity for shallow water waves of Barthelemy et al. [2018]'s no-return point ($u_c/c = 0.85$), as also demonstrated by Derakhti et al. [2020]. Finally, Varing et al. [2021b]'s findings
may explain why the observational studies investigating the kinematic criterion by tracking $u_c$ instead of $\|\mathbf{u}_m\|$, reported values less than unity.

## 2.5 Discussions and conclusions

This chapter gives an overview of the research dedicated to the identification of the breaking onset criterion. The kinematic breaking onset criterion that appears to be the most universal and accurate one is finally arising after decades of efforts dedicated to other criterion and especially of geometric type, obviously because the quantities of interest were more easily accessible from observations and models. With the continuous development of measurement methods (especially PIV) and numerical models, the kinematic criterion is becoming the reference criterion. Indeed, the development of fully nonlinear models, capturing the entire flow dynamics in the crests renders this type of breaking criterion more and more applicable.

Being able to capture breaking onset precisely is a key information to estimate the parameters of interest for slamming loads on offshore wind turbines (see Chapter 5) but also to ignite wave breaking dissipation at the right instant in the numerical wave simulations. It is finally critical to assess the breaking waves statistics as presented in the next chapter.
# Appendix

## 2.A Paper from Varing & al. (2021a)
<figure><img src="image_3.png" /></figure>

# A new definition of the kinematic breaking onset criterion validated with solitary and quasi-regular waves in shallow water

Audrey Varing<sup>a</sup>, Jean-Francois Filipot<sup>a, *</sup>, Stephan Grilli<sup>b</sup>, Rui Duarte<sup>a</sup>, Volker Roeber<sup>c</sup>, Marissa Yates<sup>d, e</sup>

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

## ARTICLE INFO

**Keywords:**

Wave breaking onset
Depth-induced breaking
Kinematic wave breaking criterion
Nearshore waves
Fully nonlinear potential flow model

## ABSTRACT

A large body of work has been devoted to the accurate detection and simulation of wave breaking in coastal areas. It is a key process for a wide range of engineering activities and environmental issues. This has motivated the development of a variety of breaking onset criteria, such as kinematic criteria based on a maximum value (usually unity) of the ratio $u_c/c$, of the horizontal particle velocity at the wave crest $u_c$ to its phase velocity $c$, both taken in the direction of wave propagation. Here, we investigate numerically the validity of this criterion in capturing breaking onset for solitary and quasi-regular two-dimensional shallow water waves using the Fully Nonlinear Potential Flow (FNPF) model by Grilli and Subramanya (1996). With this model, the propagation up to overturning of solitary waves over plane slopes, and solitary and quasi-regular waves over a submerged bar, both initially specified as numerically exact FNPF waves, is simulated. In all cases, waves break as spilling or plunging breakers, initiated by the formation of a breaker jet near the wave crest. Results show that the location of the maximum fluid velocity $||u||_m$ on the free surface closely coincides with the location where the overturning jet is initiated. Based on this, a new breaking onset criterion is proposed as $||u||_m/c \simeq 1$, which is shown to be more universal for accurately detecting wave breaking initiation than existing criteria based on the crest velocity.

## 1. Introduction

Wave breaking is a critical nearshore process for a wide range of engineering activities and environmental issues. Breaking waves cause the largest hydrodynamic loads on marine and coastal structures. In addition, energy dissipation resulting from wave breaking is also the primary control of the wave height cross-shore evolution, the wave setup magnitude, and the generation of infragravity waves and wave-induced currents that drive sediment erosion and transport (Svendsen, 1984; Longuet-Higgins and Stewart, 1964; Longuet-Higgins, 1970; MacMahan et al., 2006; Battjes et al., 2004; Symonds et al., 1982). Accordingly, a significant amount of work has been devoted to the accurate simulation or parameterization of wave breaking in nearshore numerical wave propagation models. Wave models not based directly on the primitive Navier-Stokes equations (e.g. Derakhti et al., 2018; Derakhti and Kirby, 2016; Roeber, 2010), require wave breaking to be

parameterized. Breaking must first be detected using an appropriate breaking onset criterion, and then an appropriate amount of energy dissipation must be specified over the grid area deemed to be breaking. The latter has been achieved based on different methods such as empirical eddy viscosity, sponge layer or shock capturing algorithm (e.g. Kennedy et al., 2000; Shi et al., 2012). Over the past half century, many breaking onset criteria have been proposed for a variety of practical applications based on theoretical studies, numerical simulations, laboratory experiments, or field observations (e.g Kazolea et al., 2014; Tian et al., 2010). Such criteria are typically classified into three categories: geometric, kinematic, and dynamic.

### 1.1. Geometric breaking criteria

By studying regular waves propagating over a flat bottom, Miche (1944) showed that wave breaking is initiated when $kH/\tanh(kh) =$

<footnote><sup>a</sup> France Energies Marines, Bâtiment Cap Océan, 525 avenue Alexis de Rochon, 29280, Plouzané, France</footnote>

<footnote><sup>b</sup> Department of Ocean Engineering, University of Rhode Island, Narragansett, RI, USA</footnote>

<footnote><sup>c</sup> Université de Pau et des Pays de l'Adour, E2S UPPA chair HPC-Waves, Allée du Parc Montaury, 64600, Anglet, France</footnote>

<footnote><sup>d</sup> Saint-Venant Hydraulics Laboratory, Université Paris-Est (joint research unit EDF R&D, Cerema, ENPC), 6 quai Watier, BP 49, 78401, Chatou, France</footnote>

<footnote><sup>e</sup> Cerema, Technical Department of Water, Sea and Rivers, 29280, Plouzané, France</footnote>

<footnote>* Corresponding author.
E-mail address: jean.francois.filipot@france-energies-marines.org (J.-F. Filipot).</footnote>
0.88, with $k = 2\pi/\lambda$ the wavenumber, $H$ the wave height, $\lambda$ the wavelength, and $h$ the water depth. In shallow water, Miche's criterion transforms into the breaker index criterion, $\gamma = H/h = \gamma_m$, where $\gamma_m$ is a function of the beach slope and the incident wave steepness. Widely used in the literature (e.g. Battjes and Janssen, 1978; Thornton and Guza, 1983), this approach has proven to be adequate for initiating wave breaking dissipation in phase-averaged spectral models (e.g. Filipot and Arduhuin, 1978-2012). In phase-resolving wave models, a maximum slope criterion for the wave front face has also been used to detect the onset of breaking (e.g. Grilli et al., 2019; Guignard and Grilli, 2001; Roeber, 2010; Tissier et al., 2012), but the maximum slope value used varies widely between different studies (e.g. Cienfuegos et al., 2010; Madsen et al., 1997; Schäffer et al., 1993; Sørensen et al., 1998; Tissier et al., 2012), which poses the question of the generality and applicability of such breaking onset criteria.

## 1.2. Kinematic breaking criteria

Kinematic breaking criteria are typically based on a threshold value of the ratio between the horizontal fluid velocity at the wave crest (defined as the maximum free surface elevation) $u_c$, and the local wave phase speed $c$, both taken in the direction of wave propagation (e.g. Kurnia and Van Groesen, 2014; Stansell and MacFarlane, 2002). Wave breaking is assumed to occur when $u_c$ exceeds $c$, i.e. $u/c \ge 1$, corresponding to crest overturning and the local formation of a breaker jet (see e.g. Wei et al. (1995) for an illustration of this for solitary waves overturning on plane slopes). It is important to note that, in numerical models, difficulties arise to compute accurately $u_c$ and $c$ for rapidly evolving breaking waves, and sometimes $c$ is approximated by the shallow water linear wave phase velocity $\sqrt{gh}$. Experimental studies performed in deep water (Stansell and MacFarlane, 2002; Wu and Nepf, 2002; Banner and Peirson, 2007; Tian et al., 2010) and shallow water (Itay and Liberzon, 2017) have found different thresholds for $u_c/c$, depending on the experimental configurations, ranging from 0.7 to 1.05.

## 1.3. Dynamic breaking criterion

Dynamic breaking criteria have been proposed in the last two decades, but mostly for deep water waves. Empirical findings based on observations of modulated deep water wave groups suggest that wave breaking occurs when the local wave energy flux within a group exceeds a given threshold (Song and Banner, 2002). These findings are supported by the experimental results of Tian et al. (2008, 2010) and Banner and Peirson (2007). Barthelemy et al. (2018) recently proposed, based on a numerical analysis of three-dimensional directionally focused wave packets in deep and intermediate water, a wave breaking onset criterion $B_x$ based on the ratio of the local energy flux to the energy density, normalized by the local crest speed magnitude. At the free surface, this criterion reduces to $B_x = u_c/c$ (Barthelemy et al., 2018), showed that if $B_x$ exceeded the threshold 0.85 – 0.86, wave breaking inevitably occurs within a fraction of a wave period. Thus, the threshold that they propose is not strictly speaking a breaking onset criterion, but is a precursor to a breaking onset criterion. Although their work is still ongoing, Grilli et al. (2018) recently reported that the same $B_x = 0.85 \pm 0.02$ criterion was able to predict similarly the eventual onset of wave breaking in shallow water conditions. The recent work of Derakhti et al. (2020) confirms the validity of Barthelemy et al. (2018)'s breaking onset criterion for shallow water waves.

In the present work, based on numerical simulations with the 2D Fully Nonlinear Potential Flow (FNPF) model of Grilli and Subramanya (1996), the validity of kinematic breaking onset criteria in shallow water is explored, and a new criterion is proposed that, for the waves studied here, improves breaking onset detection in comparison to the conventional kinematic criterion $u_c/c$.

The paper is organized as follows: in Section 2, the new breaking

criterion is defined, and in Section 3 the model used in the study is introduced briefly. Comparisons between $u_c/c$ and the new breaking criterion for different breaking wave conditions in shallow water are presented in Sections 4 to 6. Finally, conclusions are presented in Section 7.

# 2. Definition of a new kinematic breaking criterion

Stansell and MacFarlane (2002) studied the conventional $u_c/c$ kinematic criterion using a variety of methods to compute the wave phase speed. However, they concluded that all of the methods they used allowed confirming the validity of this criterion. Here, the objective is instead to focus on the other parameter in the criterion, the surface fluid velocity. Numerical results presented hereafter will show, for a large range of incident wave conditions and sloping bottom geometries, that when breaking is initiated by wave overturning (i.e. for spilling and plunging breakers), the maximum fluid velocity at the surface $||u||_m$, does not necessarily occur at the wave crest, but rather ahead of the crest on the front face of the wave. Breaking onset thus refers to the instant in time when part of the wave front face becomes vertical. It is thus possible to track more accurately the location where the instability leading to the overturning of the wave crest develops (see Fig. 3). This is consistent with the findings of (Perlin et al., 1996) and (Qiao and Duncan, 2001) who found, based on PIV measurements of deep water spilling and plunging breakers, that the largest flow velocity occurs along the uppermost part of the wave front, but not necessarily at the crest itself. This observation led to investigate the skill of a new kinematic criterion at capturing wave breaking onset in shallow water, defined as:

$$ \frac{||u||_m}{c} \sim 1 \qquad (1) $$

In the following, the performance of the new criterion, defined in Equation (1), is compared to that of $u_c/c \simeq 1$ for capturing accurately wave breaking onset.

# 3. The FNPF model and some definitions

## 3.1. The FNPF model

Detailed characteristics of shoaling waves near and at the breaking point will be simulated with the FNPF model of Grilli et al. (1989) and Grilli and Subramanya (1996) (hereafter referred to as 2D-NWT, i.e. 2D numerical wave tank), such as phase speed and particle velocities, and the performance of breaking onset criteria will be evaluated.

In the 2D-NWT, mass conservation is solved with a BEM using higher-order elements (at least cubic) ensuring local continuity of the first-derivatives of the geometry and field variables along the free surface. The nonlinear kinematic and dynamic free surface boundary conditions are integrated in time using an explicit second-order Eulerian-Lagrangian scheme (this formalism was first proposed by Longuet-Higgins and Cokelet (1976). An adaptive time step, based on a mesh Courant number, allows achieving accurate results when the distance between free surface nodes, here identical to Lagrangian particles, decreases. This typically occurs near and at the wave crests and subsequently within the tip of breaker jets. A summary of the model equations is presented in Appendix A, and more details about the mathematical model and numerical methods can be found in Grilli et al. (1989) and Grilli and Subramanya (1996).

Models based on FNPF theory have been shown to predict accurately the wave shape geometry of wave shoaling up to and into the early stages of wave overturning, before touchdown of the breaker jet on the free surface (before blow-up of the model). Grilli et al. (1994a, 1994b, 1997) provide extensive validation of the wave geometry with the 2D-NWT for the simulation of shoaling and breaking of solitary waves over slopes and submerged breakwaters by comparing to laboratory
data. Grilli and Horrillo (1997), Grilli (1998) and Grilli et al. (2019) provide similar results for regular waves. Wei et al. (1995) has also used this 2D-NWT as a standard of accuracy to validate their fully-nonlinear Boussinesq-type wave model.

Since the wave geometry is accurately simulated by the model at and during breaking, this indirectly validates the underlying wave kinematics at and near the surface. Therefore, the 2D-NWT model can be used to gain physical insight into complex nonlinear wave phenomena such as breaking onset based on kinematic parameters. In the present computations, the accuracy of the numerical results will be assessed by verifying the global conservation of wave volume and total energy.

## 3.2. Computation of $u_c/c$ and $||u||_m/c$ in the 2D-NWT

Computing the breaking onset criteria, $u_c/c$ and $||u||_m/c$, at each time step of the 2D-NWT simulations requires calculating the horizontal fluid velocity at the wave crest, the maximum fluid velocity on the surface, and the wave phase speed.

According to potential flow theory, the fluid velocity is given by $\mathbf{u} = \nabla\varphi$ (with $\varphi$ the velocity potential), or in a local orthogonal coordinate system, tangential and normal to the free surface, as $\mathbf{u} = (u_s, u_n)$, with the normal velocity defined as $u_n = \partial\varphi/\partial n$ and the tangential velocity as $u_s = \partial\varphi/\partial s$. Since $(\varphi, \partial\varphi/\partial n)$ are the working variables of the 2D-NWT, they are computed at free surface nodes at each time step $t$ of the simulation (using a cubic interpolation in between nodes). The tangential derivatives are then computed using a 5th-order sliding polynomial interpolation (Grilli and Subramanya, 1996). The horizontal velocity along the free surface is found by projecting the local velocity components as $u = u_s \cos\beta + u_n \sin\beta$, with $\beta$ the angle between the horizontal axis and the tangent to the free surface. At the wave crest, the latter is zero and thus $u_c = u_s$. Similarly, the velocity magnitude is computed as $||\mathbf{u}|| = \sqrt{u_s^2 + u_n^2}$, and its maximum $\mathbf{u}_m$ is calculated as the upper bound of the local maximum within each boundary element close to the crest.

The phase speed $c$ is known to be difficult to calculate accurately or to measure at arbitrary points of a non-permanent wave form (e.g. Wei et al., 1995). Stansell and MacFarlane (2002) compared three different definitions of phase speed from: (i) linear wave theory, (ii) a Hilbert transform, and (iii) the speed of a wave crest. Here, particularly since solitary waves are primarily considered, the speed of the wave crest (i.e. a local maximum in surface elevation) is used, which can be computed with a simple wave crest tracking method. Thus, in each simulation, the location of the maximum surface elevation, $x_c(t)$ is first calculated at each time $t$, using cubic interpolation in the BEM. Once model simulations are completed, a 4th-order polynomial is fit to the part of this time series including a short time before and after breaking onset, $x_{cfit}$, and the crest velocity (phase speed) is calculated analytically as $c = dx_{cfit}/dt$ from this polynomial fit.

## 3.3. Breaking onset definition

The time of breaking onset is denoted by $t_b$ and defined as the time when the front face of the wave becomes vertical within a boundary element (i.e. there is a vertical tangent with $\beta = \pi/2$). Although the cubic discretization will play a small role, for a fine spatial discretization, this will also be approximately the time a free surface node overtakes another one in the wave propagation direction. Once this occurs, the free surface slope exceeds the vertical (i.e. $\beta > \pi/2$), and the wave overturns, passing a point of no return and starting to break. The comparison of the two kinematic breaking criteria, $u_c/c$ and $\mathbf{u}_m/c$, will be done in the following applications based on this definition of breaking onset. It is interesting to note here that this particular definition of wave breaking onset may sometimes but not always correspond to other definitions, such as the wave reaching its maximum elevation.

# 4. Solitary wave propagation over a plane beach

The performance of the kinematic breaking criteria, $u_c/c$ and $\mathbf{u}_m/c$ is first investigated with the 2D-NWT for a series of simulations of solitary wave shoaling and overturning over slopes (e.g. Fig. 1). Although solitary waves are idealized limiting form (i.e. with an infinite wavelength and no trough), long waves that do not strictly occur in nature, many observations have shown that long and nonlinear shallow water swells can be approximated by a succession of solitary waves as they approach the shore since they have increasingly steep and narrow crests separated by long shallow troughs (Peregrine, 1983). A practical advantage of using solitary waves in FPNF models is that the computations are interrupted when a breaker jet impacts the free surface, which only occurs once per simulation with this type of wave.

Here, fully nonlinear solitary waves are specified as the initial condition in the constant depth ($h_0$) region of the computational domain (Fig. 1) based on the numerically exact geometry and kinematics (potential $\varphi$, normal velocity $\partial\varphi/\partial n$, and elevation $\eta$) computed with the method of Tanaka (1986). This method provides the wave properties to an accuracy of at least 9 significant digits.

In the following, prime variables denote dimensionless variables where length is normalized by $h_0$ and time by $\sqrt{h_0/g}$.

## 4.1. Model parameters

The computational domain used to simulate solitary wave shoaling and overturning is shown in Fig. 1. The length of the domain $L'$ is function of the plane slope $s = 1:100$ to $1:15$ and the initial equivalent length of the solitary wave $\lambda'$ (see Table 1). Solitary waves of model input initial height $H_0' = 0.2$ to $0.7$ were simulated with each solitary wave generated such that its minimum surface elevation is a given fraction $e_z = 0.001$ of its maximum elevation $H_0'$, which yields the value of the equivalent length $\lambda'$ (note that $\lambda'$ is also a weak function of the free surface discretization and domain geometry, so $\lambda'$ varies slightly for the same incident waves propagating over different slopes). Hence, the computational domain is longer for smaller incident waves (smaller solitary waves are longer) and milder slopes. Initial wave heights and slopes are combined to create 22 test cases of different types of wave breaking. Both the physical and numerical parameters for each of these test cases are listed in Table 1. Following Grilli et al. (1997), the corresponding surf-similarity parameter, $S_0 = 1.521s/\sqrt{H_0'}$ of these solitary waves ranges from 0.018 to 0.185, indicating spilling ($S_0 < 0.025$) or plunging ($0.025 < S_0 < 0.30$) wave breaking.

The BEM computational domain boundary is discretized with $N$ nodes, including $N_f$ nodes on the free surface. To increase the accuracy of the simulations, the latter is a large fraction of $N$ (typically over 80%, see below). Because the domain length and wavelength vary, the number of nodes per wave length (nd. per $\lambda'$) is different in each simulation. So-called cubic “Mid-Interval-Interpolation” (MII) boundary elements

<figure><img src="image_5.png" /><figcaption>Fig. 1. Definition sketch of numerical simulations with the 2D-NWT for solitary wave propagation over a plane slope.</figcaption></figure>
<table><caption>**Table 1**

Physical parameters and numerical results for computations of solitary wave shoaling over a plane beach.</caption><thead><tr><th>Simulations</th><th>s</th><th>H<sub>0</sub>'</th><th>S<sub>0</sub></th><th>λ'</th><th>nd. per λ'</th><th>L'</th><th>u<sub>c</sub>/c</th><th>||u||<sub>m</sub>/c</th><th>H<sub>b</sub>/h<sub>b</sub></th></tr></thead><tbody><tr><td>S100_H70</td><td>1:100</td><td>0.70</td><td>0.018</td><td>19.28</td><td>202</td><td>80.0</td><td>1.008</td><td>1.079</td><td>1.022</td></tr><tr><td>S100_H60*</td><td>1:100</td><td>0.60</td><td>0.020</td><td>20.24</td><td>212</td><td>80.0</td><td>0.995</td><td>1.074</td><td>1.036</td></tr><tr><td>S100_H50</td><td>1:100</td><td>0.50</td><td>0.022</td><td>21.52</td><td>184</td><td>98.0</td><td>0.994</td><td>1.064</td><td>1.043</td></tr><tr><td>S100_H40</td><td>1:100</td><td>0.40</td><td>0.024</td><td>23.16</td><td>198</td><td>98.0</td><td>0.996</td><td>1.071</td><td>1.048</td></tr><tr><td>S100_H30</td><td>1:100</td><td>0.30</td><td>0.028</td><td>25.8</td><td>188</td><td>115.0</td><td>1.016</td><td>1.085</td><td>1.049</td></tr><tr><td>S100_H20</td><td>1:100</td><td>0.20</td><td>0.034</td><td>30.47</td><td>222</td><td>115.0</td><td>1.006</td><td>1.028</td><td>1.058</td></tr><tr><td>S050_H60</td><td>1:50</td><td>0.60</td><td>0.039</td><td>20.17</td><td>260</td><td>80.0</td><td>0.992</td><td>1.084</td><td>1.181</td></tr><tr><td>S050_H40</td><td>1:50</td><td>0.40</td><td>0.048</td><td>30.32</td><td>348</td><td>80.0</td><td>0.988</td><td>1.087</td><td>1.259</td></tr><tr><td>S050_H20</td><td>1:50</td><td>0.20</td><td>0.068</td><td>30.36</td><td>318</td><td>80.0</td><td>0.995</td><td>1.072</td><td>1.257</td></tr><tr><td>S035_H60</td><td>1:35</td><td>0.60</td><td>0.056</td><td>20.16</td><td>328</td><td>51.5</td><td>0.978</td><td>1.081</td><td>1.314</td></tr><tr><td>S035_H40</td><td>1:35</td><td>0.40</td><td>0.069</td><td>23.17</td><td>292</td><td>66.5</td><td>0.966</td><td>1.079</td><td>1.373</td></tr><tr><td>S035_H30</td><td>1:35</td><td>0.30</td><td>0.079</td><td>25.87</td><td>326</td><td>66.5</td><td>0.968</td><td>1.085</td><td>1.407</td></tr><tr><td>S035_H20</td><td>1:35</td><td>0.20</td><td>0.097</td><td>30.31</td><td>382</td><td>66.5</td><td>0.977</td><td>1.101</td><td>1.462</td></tr><tr><td>S020_H70</td><td>1:20</td><td>0.70</td><td>0.091</td><td>19.23</td><td>424</td><td>38.0</td><td>0.914</td><td>1.026</td><td>1.649</td></tr><tr><td>S020_H60</td><td>1:20</td><td>0.60</td><td>0.098</td><td>20.22</td><td>446</td><td>38.0</td><td>0.922</td><td>1.046</td><td>1.724</td></tr><tr><td>S020_H50</td><td>1:20</td><td>0.50</td><td>0.108</td><td>51.41</td><td>390</td><td>46.0</td><td>0.931</td><td>1.061</td><td>1.792</td></tr><tr><td>S020_H40</td><td>1:20</td><td>0.40</td><td>0.120</td><td>23.16</td><td>422</td><td>46.0</td><td>0.934</td><td>1.073</td><td>1.870</td></tr><tr><td>S020_H20</td><td>1:20</td><td>0.20</td><td>0.170</td><td>30.23</td><td>478</td><td>53.0</td><td>0.892</td><td>1.046</td><td>2.180</td></tr><tr><td>S015_H60</td><td>1:15</td><td>0.60</td><td>0.131</td><td>20.18</td><td>376</td><td>44.985</td><td>0.869</td><td>1.003</td><td>2.176</td></tr><tr><td>S015_H50</td><td>1:15</td><td>0.50</td><td>0.143</td><td>21.47</td><td>400</td><td>44.985</td><td>0.868</td><td>1.012</td><td>2.296</td></tr><tr><td>S015_H40*</td><td>1:15</td><td>0.40</td><td>0.160</td><td>23.19</td><td>432</td><td>44.985</td><td>0.866</td><td>1.018</td><td>2.441</td></tr><tr><td>S015_H30</td><td>1:15</td><td>0.30</td><td>0.185</td><td>25.77</td><td>480</td><td>44.985</td><td>0.845</td><td>1.007</td><td>2.669</td></tr></tbody></table>

are used to interpolate the solution in between nodes on the free surface, and 3-node quadratic isoparametric elements on other parts of the boundary (Grilli and Subramanya, 1996). The bottom discretization is typically much coarser than the free surface discretization, but to

increase accuracy in shallower water, the discretization is non-uniform over the slope with the distance between nodes decreasing gradually with decreasing depth. Consistent with the Eulerian-Lagrangian formalism, the BEM nodes on the free surface gradually converge near

<figure><img src="image_6.png" /><figcaption>**Fig. 2.** Absolute maximum numerical errors with respect to initial values computed with the method of Tanaka (1986), as a function of $S_0$ for the wave volume $||ε_v||$ and total energy $||ε_e||$.</figcaption></figure>

<footnote>*Simulations appearing in Fig. 3.</footnote>
the crest region during the shoaling phase preceding breaking onset, which also increases the accuracy of the results in shallower water.

A convergence study as a function of the discretization, detailed in Appendix B, was first performed to establish the minimum discretization required to ensure both accurate and converged results, particularly for the parameters of interest (flow velocity and phase speed) at breaking onset. Shoaling of a solitary wave with $H'_{0} = 0.6$ was simulated over 1:15 and 1:100 slopes, and maximum relative errors of the solitary wave volume and total energy (as compared to initial values computed with the method of Tanaka (1986)) were computed, together with the values of $u_c/c$, $u_m/c$, and $c$ at breaking onset. Nine different discretizations were tested for each case and it was found that using $N = 1038$ and $N_f = 838$ ensured that both errors remained smaller than 0.01% while the values of the breaking onset parameters had clearly converged to within less than 0.2% at least three configurations coarser than the selected one ($N = 1038$ and $N_f = 838$). These values of $N$ and $N_f$ were thus used in the 22 test cases studied here, with the model parameters and results shown in Table 1. [Note that in this table $H_b$ and $h_b$ denote the breaking wave height (crest elevation for a solitary wave) and water depth, respectively, and $\gamma = H_b/h_b$ is the breaker index, at $t_b$.]

To assess further the accuracy of the numerical results for these cases, the maximum relative errors of wave volume and energy conservation were also calculated up to breaking onset. Fig. 2 shows the maximum absolute errors of wave energy $||\epsilon_e||$ and volume $||\epsilon_v||$, calculated for wave propagation from $t = 0$ to $t_b$. At most, the maximum error for energy reaches 0.2% and for volume 0.015%, which is sufficient to ensure reliable estimates of the breaking onset parameters ($u_c$, $u_m$, $c$) based on the convergence study in Appendix B. [Note that this convergence study also showed that the breaking onset parameters converged for coarser discretizations (and hence for larger numerical errors) than achieved with this discretization.]

## 4.2. Comparison of $u_c/c$ and $||u||_m/c$ breaking onset criteria

Some of the solitary wave cases simulated in this study (parameters listed in Table 1) were similar to the cases modeled and analyzed by Grilli et al. (1997). Additional simulations were performed for different slopes (e.g. 1:20 and 1:50) and incident wave heights. Fig. 3 shows the free surface elevation simulated with the 2D-NWT at four different time steps (with the second being $t_b$) for a mild spilling breaker (top panel $H'_{0} = 0.6$ on a 1:100 slope) and for a more intense plunging breaker (bottom panel, $H'_{0} = 0.4$ on a 1:15 slope). The symbols on each free

<figure><img src="image_7.png" /><figcaption>Fig. 3. Solitary wave propagation simulated with the 2D-NWT up to overturning of the wave crest for: (top panel, spilling breaker) a wave with an initial height $H'_{0} = 0.6$ over a 1:100 slope, and (bottom panel, plunging breaker) a wave with an initial height $H'_{0} = 0.4$ over a 1:15 slope. Curves a,b,c,d are results obtained for increasing times (discussed in the text). In each case, curve b corresponds to the breaking onset time $t_b$, at which the front face of the wave becomes vertical. The locations of the crest (maximum elevation) and maximum velocity $||u||_m$ are marked by circle and square symbols, respectively.</figcaption></figure>

surface profile mark the location of the maximum elevation (i.e. wave crest) and that of the maximum velocity $u_m$. They are in general different, although they are much closer for the spilling breaker than for the plunging breaker. At $t_b$, $||u||_m$ occurs approximately at the location where the wave front face is vertical and a breaker jet will soon be emerging, leading to wave overturning.

Fig. 4 shows the breaking onset parameter values $u_c/c$ and $||u||_m/c$ computed for each case in Table 1, either as a function of time $t' - t'_b$, or at $t' = t'_b$ as a function of $S_0$. Both $u_c/c$ and $||u||_m/c$ increase continuously with time and, because the maximum fluid velocity is larger than the horizontal velocity at the crest, $||u||_m/c$ reaches larger values at $t'_b$ than $u_c/c$. Fig. 4 (b) and (d) show that at breaking onset, $u_c/c \in [0.845-1.015]$ while $u_m/c \in [1.003-1.101]$. More specifically, $u_c/c$ value at breaking onset decreases with an increasing $S_0$. For only spilling breakers and small plunging breakers (i.e. $S_0 < 0.8$), the 95% confidence interval of $u_c/c$ is $u_c/c = 0.992 \pm 0.028$ at breaking onset, but this value becomes much smaller for intense plunging breakers, $u_c/c = 0.902 \pm 0.075$. Considering all test cases, $u_c/c = 0.951 \pm 0.104$. In contrast, the corresponding 95% confidence interval of $u_m/c$ values at breaking onset is narrower, with $||u||_m/c \simeq 1.058 \pm 0.057$ for all the cases in Table 1, which provides a good prediction of breaking onset even for strong plunging waves.

This figure highlights two significant findings. First, consistent with the work of (Derakhti et al., 2020), it shows that the results of Barthelemy et al. (2018) can be extended to shallow water: the parameter $B_x = u_c/c \approx 0.85$ is thus an accurate universal precursor to breaking initiation. Second, $||u||_m/c$ also agrees with this finding and is even closer to $||u||_m/c \approx 1$ at breaking initiation, which likely results from $||u||_m$ capturing better the location where overturning starts as compared to $u_c$.

Finally, the differences between $||u||_m$ and the maximum horizontal particle velocity $u_m$ on the free surface (not detailed in the paper) showed that $u_m$ and $||u||_m$ are almost identical until breaking, as the vertical component of the free surface velocity $w$ is very small in the crest area. Then $||u||_m$ becomes larger than $u_m$ a short time before breaking, particularly for plunging breakers, and the difference between these velocities continues to increase beyond breaking onset.

# 5. Solitary wave propagation over a submerged bar

In this application, the variation of $u_c/c$ and $u_m/c$ is investigated for solitary waves shoaling over a submerged bar with a 1:20 front slope, a crest at $0.75h_0$, and a 1:10 back slope. The bar crest length is $5h_0$ (see Fig. 5). Beji and Battjes (1994) first proposed this bar geometry and used it to perform both laboratory experiments and numerical simulations of regular wave shoaling and energy transfers between higher-order harmonics. Grilli et al. (2019) validated the 2D-NWT with measurements from several of the test cases in these experiments. The present study uses the same bar geometry but instead evaluates the properties of breaking solitary waves propagating over the bar.

## 5.1. Model parameters

Table 2 provides the physical and model parameters for the simulated test cases, as well as the numerical results of the estimated breaking onset parameter values. Fig. 5 shows the computational domain used in the 2D-NWT. Here, to reduce the computational effort, the domain is truncated, with a vertical wall specified at the toe of the bar back slope ($x' = 57.5$). The computations were stopped before any reflections from the far end of the 2D-NWT could propagate back into the area of interest over the bar. The model input wave height in the different test cases varied from $H'_{0} = 0.1106$ to 0.3, with corresponding incident equivalent wavelengths varying from $\lambda' = 20.04$ to 29.78. The BEM discretization for this case was selected following the same approach as for the previous application with a plane slope. The number of nodes per wavelength varies from 292 to 434 on the free surface.
<figure><img src="image_8.png" /><figcaption>Fig. 4. Computations of solitary wave shoaling and overturning over a plane slope. Breaking onset criteria computed: (a,c) as a function of $t' - t'_b$ and (b,d) at $t'_b$ as a function of the surf-similarity parameter $S_0$. Each breaking wave is identified by a color and label identical to those listed in Table 1, which shows the related parameters values and results. Labels indicate both slope and wave height with, for example, S020_H70 referring to a wave with $H'_0 = 0.70$ propagating over a 1:20 slope. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)</figcaption></figure>

In this application, both breaking and nearly breaking solitary waves are simulated. For the latter, the waves shoal and nearly reach breaking onset over the bar crest, but do not break and continue propagating beyond the bar. With this bar geometry, numerical experiments showed that incident solitary waves with $H'_0 < 0.1175$ did not break (e.g. Fig. 6, top panel), whereas those with larger incident wave heights did break (e. g. Fig. 6, bottom panel).

## 5.2. Comparison of $u_c/c$ and $||u||_m/c$ breaking onset criteria

Fig. 7 compares the breaking onset parameter values $u_c/c$ and $u_m/c$ computed for each case in Table 2, either as a function of time $t' - t'_b$, or at $t' = t'_b$ as a function of $H'_0$. For non-breaking cases, the reference time used (instead of $t_b$) corresponds to the time the crest reaches its maximum elevation. Fig. 7 (a) shows that the conventional breaking onset criterion $u_c/c \simeq 1$ does not allow distinguishing between breaking and non-breaking waves. Further, consistent with the findings of
<figure><img src="image_9.png" /><figcaption>Fig. 5. Computational domain used in 2D-NWT simulations of solitary wave propagating over a submerged bar.</figcaption></figure>

Barthelemy et al. (2018) and Derakhti et al. (2020), Fig. 7 (b) shows that only waves for which $u_c/c$ becomes larger than 0.85 during shoaling on the bar will evolve towards breaking (see also Table 2). In Fig. 7 (c) and (d), this appears to apply to $u_m/c$ as well, which remains below 0.85 for non-breaking waves. This is likely because the location of $u_m$ is very close to the crest for the considered waves, hence $||u||_m \approx u_c$ for waves that do not evolve towards breaking. In contrast, for breaking wave cases, the 95% confidence intervals are computed as $u_c/c \approx 0.906 \pm 0.044$ and $u_m/c \approx 1.032 \pm 0.044$ at $t_b$. Therefore, while $u_c/c \approx 0.85$ allows distinguishing breaking from non-breaking cases, this difference appears to be more pronounced with $u_m/c \approx 1$.

# 6. Quasi-regular wave propagation over a submerged bar

Finally, the shoaling of quasi-regular waves over the same submerged bar as in the previous application is simulated to assess whether the same conclusions apply regarding the breaking onset criteria. Running simulations with quasi-regular (or even irregular) waves in the 2D-NWT is possible, but becomes more computationally expensive.

## 6.1. Model parameters

The physical and model parameters for the simulated test cases are shown in Table 3. Fig. 8 shows the computational domain used in the 2D-NWT. As compared to the previous computational domain used to propagate and shoal solitary waves, the length of the constant depth region offshore of the bar is reduced such that the toe of the back slope of the bar is now at $x' = 42.5$, and an absorbing beach (AB) of length $l'_{AB} = 20$ is specified to minimize reflections from the far-end of the 2D-NWT (see Grilli and Horrillo (1997) for more details). Numerically exact FNPF regular waves are generated on the offshore boundary, as “zero-mass-flux” stream function waves (see Grilli and Horrillo (1997) for more details). The model input wave period is $T' = 4$ in all cases, and the model input wave height is varied from $H'_0 = 0.06$ to 0.2 with an incident wavelength of approximately $\lambda' = 12$. The BEM discretization was selected following the same approach as in previous applications. The number of nodes per wavelength is approximately 160 on the free surface. Since the 2D-NWT simulation is terminated once a breaker jet impacts the free surface, a fully stationary regular wave train can not be

<caption>**Table 2**</caption>

Physical parameters and numerical results for computations of solitary waves propagating over a submerged bar. The domain length is $L' = 57.5$ in all cases.

<table><thead><tr><th>Simulations</th><th>H'<sub>0</sub></th><th>&lambda;'</th><th>nd. per &lambda;'</th><th>u<sub>c</sub>/c</th><th>||u||<sub>m</sub>/c</th><th>Breaking</th><th>H<sub>b</sub>/h<sub>b</sub></th></tr></thead><tbody><tr><td>HBS1160</td><td>0.1106</td><td>29.78</td><td>434</td><td>0.776</td><td>0.782</td><td>No</td><td>0.52</td></tr><tr><td>HBS1165*</td><td>0.1165</td><td>29.78</td><td>434</td><td>0.793</td><td>0.801</td><td>No</td><td>0.52</td></tr><tr><td>HBS1170</td><td>0.1170</td><td>29.64</td><td>432</td><td>0.826</td><td>0.834</td><td>No</td><td>0.53</td></tr><tr><td>HBS1175</td><td>0.1175</td><td>29.64</td><td>432</td><td>0.874</td><td>1.016</td><td>Yes</td><td>0.51</td></tr><tr><td>HBS1500</td><td>0.15</td><td>26.62</td><td>388</td><td>1.080</td><td>1.006</td><td>Yes</td><td>1.17</td></tr><tr><td>HBS2000*</td><td>0.2</td><td>23.60</td><td>344</td><td>0.926</td><td>1.050</td><td>Yes</td><td>1.34</td></tr><tr><td>HBS2500</td><td>0.25</td><td>21.55</td><td>314</td><td>0.885</td><td>1.025</td><td>Yes</td><td>1.53</td></tr><tr><td>HBS3000</td><td>0.3</td><td>20.04</td><td>292</td><td>0.923</td><td>1.067</td><td>Yes</td><td>1.74</td></tr></tbody></table>

<figure><img src="image_10.png" /><figcaption>Fig. 6. Simulated free surface elevation of solitary waves propagating over a submerged bar: (top panel) $H'_0 = 0.1165$ (non-breaking case) and (bottom panel) $H'_0 = 0.2$ (breaking case). Curves a,b,c represent different time steps. In the top panel, curve b corresponds to the time the crest reaches its maximum elevation, while in the bottom panel curve b corresponds to $t = t_b$, i.e. breaking onset.</figcaption></figure>

achieved. In the following, we will therefore only study the first breaking wave in the wave train, which is thus quasi-regular and not purely regular for the reasons explained above.

Both breaking and nearly breaking waves are investigated. For the non-breaking cases, only the first wave of the quasi-regular wave train is studied for consistency with the breaking cases. Numerical experiments showed that, for this bar geometry, incident regular parameters with $H'_0 < 0.0708$ ($H'_0$ is the input parameter of the model) did not break (e.g. Fig. 9 top panel), whereas those with a larger incident wave heights did break (e.g. Fig. 9 bottom panel).

Fig. 10: Breaking onset criteria computed for the first wave in the wave train of regular waves (period $T' = 4$) propagating over a submerged bar: (a,c) as a function of $t' - t'_b$ and (b,d) at $t'_b$ as a function of $H'_0$. Each case is identified by a color and a label identical to those listed in Table 3, which shows the related parameters values and results. Note, for the non-breaking cases, the reference time used (instead of $t_b$) corresponds to the time the crest reaches its maximum elevation.

## 6.2. Comparison of $u_c/c$ and $||u||_m/c$ breaking onset criteria

Similar to solitary waves, Fig. 10 (a) and (b) show that the conventional breaking onset criterion $u_c/c \approx 1$ does not allow distinguishing between breaking and non-breaking waves. Consistent with the findings of Barthelemy et al. (2018) and Derakhti et al. (2020), Fig. 10 (b) and (d) show that only waves for which $u_c/c$ or $||u||_m/c$ becomes larger than 0.85 during shoaling over the bar will evolve towards breaking (see also Table 3). For breaking wave cases, the 95% confidence intervals are wider than for solitary waves, $u_c/c \approx 0.926 \pm 0.173$ and $||u||_m/c \approx 1.027 \pm 0.09$ at $t_b$. $||u||_m/c$ is closer to a constant value than $u_c/c$.

# 7. Discussion and conclusions

In this study, numerical simulations were performed with a 2D-NWT to assess the ability of two kinematic breaking criteria $u_c/c$ and $u_m/c$ to predict correctly wave breaking onset for solitary and quasi-regular waves shoaling and breaking as spilling or plunging breakers over plane slopes or a mildly sloping bar.

The two breaking criteria only differ in the definition of the fluid velocity considered at the free surface: the first criteria $u_c/c$ uses the horizontal component of the fluid velocity at the wave crest, while the second criteria $u_m/c$ uses the maximum fluid velocity on the front face of the wave. These two velocities are usually quite similar up to close to breaking onset, but may then differ significantly. Here, the results show that the maximum velocity at the free surface occurs closer to the

<footnote>*Simulations appearing in Fig. 6.</footnote>
<figure><img src="image_11.png" /><figcaption>**Fig. 7.** Breaking onset criteria computed for the solitary waves propagating over a submerged bar: (a,c) as a function of $t' - t'_b$ and (b,d) at $t'_b$ as a function of $H'_0$. Each wave is identified by a color and a label identical to those listed in Table 2, which shows the related parameters values and results. Note, for the non-breaking cases, the reference time used (instead of $t_b$) corresponds to the time the crest reaches its maximum elevation. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>Fig. 8. Computational domain used in 2D-NWT simulations of quasi-regular wave propagating over a submerged bar. An absorbing beach of length $l'_{AB}$ is specified at the far-end of the domain, starting at $x' = 42.5$.</figcaption></figure>

<caption>Table 3</caption>

Physical parameters and numerical results for computations of regular waves propagating over a submerged bar (Fig. 5). The wave period is $T' = T/\sqrt{g/h_0} = 4$ and the domain length is $L' = 62.4$ for all cases.

<table><thead><tr><th>Simulations</th><th>H'<sub>0</sub></th><th>λ'</th><th>nd. per λ</th><th>u<sub>c</sub>/c</th><th>||u||<sub>m</sub>/c</th><th>Breaking</th><th>H<sub>b</sub>/h<sub>b</sub></th></tr></thead><tbody><tr><td>HBR0600</td><td>0.060</td><td>12.00</td><td>161</td><td>0.460</td><td>0.460</td><td>No</td><td>0.41</td></tr><tr><td>HBR0650*</td><td>0.065</td><td>12.00</td><td>161</td><td>0.552</td><td>0.553</td><td>No</td><td>0.48</td></tr><tr><td>HBR0700</td><td>0.070</td><td>12.00</td><td>161</td><td>0.751</td><td>0.763</td><td>No</td><td>0.58</td></tr><tr><td>HBR0708</td><td>0.0708</td><td>12.00</td><td>161</td><td>1.092</td><td>1.100</td><td>Yes</td><td>0.64</td></tr><tr><td>HBR0900</td><td>0.090</td><td>12.02</td><td>161</td><td>0.942</td><td>1.062</td><td>Yes</td><td>0.69</td></tr><tr><td>HBR1200*</td><td>0.12</td><td>12.05</td><td>162</td><td>0.861</td><td>0.999</td><td>Yes</td><td>0.81</td></tr><tr><td>HBR1500</td><td>0.15</td><td>12.10</td><td>162</td><td>0.863</td><td>0.984</td><td>Yes</td><td>0.84</td></tr><tr><td>HBR2000</td><td>0.2</td><td>12.18</td><td>163</td><td>0.873</td><td>0.989</td><td>Yes</td><td>0.84</td></tr></tbody></table>

location of the initiation of the overturning jet than the wave crest velocity (i.e. at maximum elevation). This finding, which is supported by some experimental studies (Perlin et al., 1996; Qiao and Duncan, 2001), motivated the present investigation of the performance of the $u_m/c$ criterion in capturing the onset of wave breaking.

In the configurations considered in this paper, including the propagation of solitary waves over a plane beach, and solitary and quasi-regular waves over a submerged bar, $||u||_m/c$ was found to predict the breaking onset with a higher accuracy than $u_c/c$, which varied over a wider range. Considering all the results summarized in Tables 1–3, the Root-Mean-Square (RMS) errors have been computed and they are minimized considering the empirical thresholds 0.95 and 1.05 for $u_c/c$ and $||u||_m/c$ respectively. The RMS errors are then 6.3 and 3.2%, with respect to $u_c/c = 0.95$ and $||u||_m/c = 1.05$, respectively. In addition to this improvement in terms of RMS errors, the results here suggest that the initiation of overturning for spilling or plunging breakers is captured more accurately when considering the location along the free surface of the maximum fluid velocity.

These findings may help explain the scatter observed in breaking thresholds reported in the literature for kinematic breaking criteria based on the horizontal fluid velocity at the crest (e.g. Banner and Peirson, 2007; Saket et al., 2018; Stansell and MacFarlane, 2002; Tian et al., 2010; Wu and Nepf, 2002). However, the results presented here also showed that using $u_c$ in the criterion is acceptable when considering the onset of spilling and small plunging breakers. For these cases, $u_c/c$ is indeed close to 1 since the breaker jet forms close to the crest location. In contrast, for strong plunging breakers, $u_c/c < 1$ at breaking onset, since $u_c$ underestimates the velocity representative of the breaker jet. In such cases, the plunging jet forms farther from the crest location. Therefore, the new kinematic criterion based on the maximum velocity observed along the free surface may be more universal, allowing improvements in the prediction of breaking onset for both spilling and plunging wave breaking. In all of the cases of breaking and non-breaking waves presented here, the results are consistent with the conclusions of Barthelemy et al. (2018) that waves with $B_x = u_c/c > 0.855 \pm 0.005$ will inevitably evolve towards breaking. This also appears to apply to $||u||_m/c$.

<figure><img src="image_13.png" /><figcaption>Fig. 9. Simulated free surface elevations for quasi-regular waves of period $T' = 4$ shoaling over a submerged bar: (top panel) $H'_{0} = 0.065$ (non-breaking case) and (bottom panel) $H'_{0} = 0.12$ (breaking case). Curves a,b,c represent different time steps. In the top panel, curve b corresponds to the time the crest reaches its maximum elevation, while in the bottom panel curve b corresponds to $t = t_b$, i. e. the time of breaking onset.</figcaption></figure>

This work demonstrates that $||u||_m/c \approx 1$ can be used as a robust and accurate breaking onset criterion for identifying the initiation of breaking for both solitary and quasi-regular breaking waves in shallow water conditions.

The implementation of this new breaking onset criterion thus requires calculating $||u||_m$ and $||c||$ throughout the numerical domain. As previously mentioned, Stansell and MacFarlane (2002) confirmed the validity of a variety of different methods for calculating the phase velocity $c$. However, the calculation of the maximum fluid velocity at the free surface $||u||_m$ deserves more attention, and the implemented method depends on the type of wave propagation model considered. For example, this term can be calculated from the velocity potential at the free surface in FNFP models using a variety of different approaches (e.g. High Order Spectral (HOS) (Gouin et al., 2016; Seiffert et al., 2017), finite difference (Bingham and Zhang, 2007), finite difference-spectral methods (Raoult et al., 2016), or coupled local modes (Belibassakis and Athanassoulis, 2011), among others) to solve the Zakharov equations (Zakharov, 1968) at the free surface. Seiffert et al. (2017) already successfully implemented Barthelemy et al. (2018)'s criterion in a FNFP-HOS model, which supports the practical applicability of our new definition of the kinematic criterion to phase-resolving wave models, provided there is an additional step in the computations, to identify $u_m$ instead of $u_c$. In the case of commonly used Boussinesq-type models (e.g. Kennedy et al., 2000b; Roeber et al., 2010; Shi et al., 2011; Kazolea et al., 2014, among others), the horizontal fluid velocity is solved at an arbitrary depth (close to mid-depth in general), and thus an additional method is required to estimate the velocity at the free surface, such as extrapolating the mid-depth velocity to the free surface using a quadratic velocity profile. This extrapolation can be challenging, since the necessary terms to compute the velocity at a level other than that around mid-depth involve second-order derivatives. Especially over an irregular bathymetry with an irregular free surface in the 2D horizontal plane, the numerical solution can be quite spiky, which can lead to unreliable values of the recomputed velocity at the free surface. Similar to the $B_x$ breaking criterion, the applicability of our new breaking criterion is thus directly related to the accuracy of this parametric vertical reconstruction, but in our case with an additional step to identify the location of the maximum velocity at each computational time step. The application of these types of breaking onset criteria in weakly nonlinear models should be further investigated. Using the newly proposed criterion instead of $B_x$ may still allow achieving a higher nonlinearity before dissipation is triggered in fully nonlinear wave models.

Also note that the validity of $||u||_m/c$ in deep water conditions and for irregular and 3D waves has not been assessed yet. Grilli et al. (1997)

<footnote>*Simulations appearing in Fig. 9.</footnote>
<figure><img src="image_14.png" /><figcaption>**Fig. 10.** Compares the breaking onset parameter values $u_c/c$ and $||u||_m/c$ computed for each case in Table 3, either as a function of time $t' - t_b$, or at $t' = t_b$ as a function of $H_0$. As in the previous section, for non-breaking cases, the reference time used (instead of $t_b$) corresponds to the time the crest reaches its maximum elevation. The evolution of the wave breaking onset criteria is recorded for the first wave of the wave train in each case.</figcaption></figure>

have shown that there exists self-similarity of breaking geometries and properties, indicating a “local loss of memory” of the flow for the phenomenon that has caused breaking, suggesting that whether the pre-breaking waves are solitary, quasi-regular, or even irregular should not significantly affect the wave crest flow velocities at breaking onset. Still, rigorous validation of our new definition of the kinematic criterion in other situations would require further research.

## Acknowledgments

This work has been done in the framework of the DiMe and CAR-AVELE projects, which benefit from French government support managed by the Agence Nationale de la Recherche (ANR) under the program "Investissement d'Avenir" ANR-10-IEED-0006-14 and ANR-10-IEED-0006-26. The numerical simulations were performed on the High Performance Computing (HPC) cluster Datarmor of the "Pôle de Calcul Intensif pour la Mer" (PCIM) (http://www.ifremer.fr/pcim). Volker Roeber acknowledges financial support from the Isite program Energy &
## Appendix A - Summary of the numerical model governing equations and numerical implementation

The velocity potential $\varphi(x, t)$ is used to represent inviscid irrotational 2D flows in the vertical plane $(x,z)$, and the velocity is defined by $\mathbf{u} = \nabla\varphi = (u, w)$. The continuity equation in the fluid domain $\Omega(t)$ with boundary $\Gamma(t)$ is the Laplace equation for the potential:

$$ \nabla^2 \varphi = 0 \quad \text{in } \Omega(t) \qquad (A1) $$

The 2D-NWT model simulates 2D free-surface flows of an ideal fluid. Under such conditions, Green's second identity (with free space Green's function $G(x, x_l) = -(1/2\pi)\log|x - x_l|$) makes it possible to transform the continuity equation for the velocity potential into a boundary-integral equation (BIE).

$$ \alpha(\mathbf{x}_l)\varphi(\mathbf{x}_l) = \int_{\Gamma_l(\mathbf{x})} \left[ \frac{\partial\varphi}{\partial n}(\mathbf{x})G(\mathbf{x}, \mathbf{x}_l) - \varphi(\mathbf{x})\frac{\partial G(\mathbf{x}, \mathbf{x}_l)}{\partial n} \right] d\Gamma_l(\mathbf{x}) \qquad (A2) $$

in which $\mathbf{x} = (x, z)$ and $\mathbf{x}_l = (x_l, z_l)$ are position vectors for points on the boundary, $\mathbf{n}$ is the unit outward normal vector, and $\alpha(\mathbf{x}_l)$ is a geometric coefficient. Equation (A2) is solved by a BEM using a set of collocation nodes on the boundary and high-order elements to interpolate in between the collocation nodes. Integrals in A2 are evaluated numerically, and the resulting algebraic system of equations is assembled and solved for the equivalent discretized problem. On the free surface $\Gamma_f(t)$, $\varphi$ satisfies the nonlinear kinematic and dynamic boundary conditions,

$$ \frac{Dr}{Dt} = \mathbf{u} = \nabla\varphi \quad \text{on } \Gamma_f(t) \qquad (A3) $$

$$ \frac{D\varphi}{Dt} = -gz + \frac{1}{2}\nabla\varphi \cdot \nabla\varphi - \frac{p_a}{\rho} \quad \text{on } \Gamma_f(t) \qquad (A4) $$

respectively, with $\mathbf{r}$ the position vector of a free surface fluid particle, $g$ the acceleration of gravity, $z$ the vertical coordinate (positive upwards, with $z=0$ at the undisturbed free surface), $p_a$ the atmospheric pressure, $\rho$ the fluid density, and the material derivative being defined as,

$$ \frac{D}{Dt} \equiv \frac{\partial}{\partial t} + \mathbf{u} \cdot \nabla \qquad (A5) $$

Along the stationary bottom $\Gamma_B$ and the other fixed boundary $\Gamma_{r2}$, a no-flow condition is prescribed as,

$$ \frac{\partial \varphi}{\partial n} = 0 \quad \text{on } \Gamma_B \text{ and } \Gamma_{r2} \qquad (A6) $$

The time updating is performed by integrating the fully nonlinear free-surface boundary conditions A3 and A4 using a Lagrangian Taylor series expansions of the free surface position $\mathbf{r}$ and potential $\varphi$.

## Appendix B. - Convergence of the model results as a function of the discretization

Here, a convergence analysis is presented for two solitary waves of incident height $H'_{0} = 0.6$ shoaling over different plane slopes (see Fig. 1). The configuration $S100\_H60$, with $s = 1:100$, leads to spilling wave breaking, while $S015\_H60$, with $s = 1:15$, leads to plunging wave breaking. These two configurations have different spatial resolutions since their domain lengths $L'$ also differ (see Table 1). The sensitivity of the breaking onset parameters is evaluated for a variety of spatial resolutions in order to verify that the simulations are both accurate in terms of mass and energy conservation and have adequately converged.

The reference discretization of the BEM grid is the highest resolution case, which has $N = 1038$ nodes distributed along the free surface ($N_f = 838$ nodes) and the bottom, left, and right boundaries of the computational domain. Here, different spatial resolutions are tested, defined as a percentage of nodes with respect to the reference discretization (which thus has 100% of the nodes).

<caption>Table B1</caption>

Convergence of the position of wave breaking onset $x_c$ as a function of the computational domain discretization for the $S100\_H60$ configuration.

<table><thead><tr><th>Simulations</th><th>N<sub>f</sub></th><th>nd. per λ</th><th>x<sub>c</sub> at t<sub>0</sub></th></tr></thead><tbody><tr><td>100% of nodes</td><td>838</td><td>212</td><td>44.47</td></tr><tr><td>90% of nodes</td><td>755</td><td>191</td><td>44.46</td></tr><tr><td>80% of nodes</td><td>671</td><td>170</td><td>44.47</td></tr><tr><td>70% of nodes</td><td>587</td><td>149</td><td>44.40</td></tr><tr><td>60% of nodes</td><td>503</td><td>127</td><td>44.39</td></tr><tr><td>50% of nodes</td><td>419</td><td>106</td><td>44.39</td></tr><tr><td>40% of nodes</td><td>335</td><td>85</td><td>44.21</td></tr><tr><td>30% of nodes</td><td>251</td><td>64</td><td>43.95</td></tr><tr><td>20% of nodes</td><td>167</td><td>42</td><td>41.72</td></tr></tbody></table>
<caption>Table B2</caption>

Convergence of the position of wave breaking onset $x_c$ as a function of the computational domain discretization for the S015_H60 configuration.

<table><thead><tr><th>Simulations</th><th>N<sub>f</sub></th><th>nd. per λ</th><th>x<sub>c</sub> at t<sub>b</sub></th></tr></thead><tbody><tr><td>100%of nodes</td><td>838</td><td>376</td><td>40.24</td></tr><tr><td>90%of nodes</td><td>755</td><td>339</td><td>40.24</td></tr><tr><td>80%of nodes</td><td>671</td><td>301</td><td>40.25</td></tr><tr><td>70%of nodes</td><td>587</td><td>263</td><td>40.24</td></tr><tr><td>60%of nodes</td><td>503</td><td>226</td><td>40.23</td></tr><tr><td>50%of nodes</td><td>419</td><td>188</td><td>40.21</td></tr><tr><td>40%of nodes</td><td>335</td><td>150</td><td>40.17</td></tr><tr><td>30%of nodes</td><td>251</td><td>116</td><td>39.91</td></tr><tr><td>20%of nodes</td><td>167</td><td>77</td><td>37.69</td></tr></tbody></table>

<figure><img src="image_15.png" /><figcaption>Fig. B1. Convergence analysis of the flow velocities $u_c$ (left) and $u_m$ (center) and of the wave phase velocity $c$ (right) at $t_b$ as a function of the BEM grid discretization for the two beach slope configurations.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Fig. B2. Convergence analysis of the breaking onset criteria $u_c/c$ (left) and $u_m/c$ (right) at $t_b$ as a function of the BEM grid discretization for the two beach slope configurations.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>**Fig. B3.** Convergence analysis of the maximum errors of mass (||ε<sub>v</sub>||, left) and energy (||ε<sub>e</sub>||, right) of the solitary wave at t<sub>b</sub> as a function of the BEM grid discretization for the two beach slope configurations.</figcaption></figure>

Figs. B1, B2 and B3 and Tables B1 and B2 show the simulation results for both beach slope configurations and all of the grid discretizations for the flow velocities u<sub>c</sub> and u<sub>m</sub>, the wave phase velocity c, the breaking onset criteria u<sub>c</sub>/c and ||u||<sub>m</sub>/c, and the maximum relative errors (in absolute value) of the wave volume and energy, ||ε<sub>v</sub>|| and ||ε<sub>e</sub>||, with respect to the initial values computed with the method of Tanaka (1986).

For the S015_H60 configuration, reference 100% node discretization has a very fine resolution, so the parameters in Figures B1 and B2 start converging for simulation with only a 60% node discretization. In contrast, the parameters for the S100_H60 configuration only start converging for the simulation with a 80% node discretization. Importantly, Figure B3 shows that both relative numerical errors on mass and energy of the incident wave decrease significantly to acceptable levels for the converged simulations (to less than or about 0.01%).

For the two configurations tested here, the convergence analysis results show that the reference discretization provides both converged and accurate results at t<sub>b</sub>. This discretization is the basis for all the simulations considered in this paper.

## References

* Banner, M.L., Peirson, W.L., 2007. Wave breaking onset and strength for two-dimensional deep-water wave groups. J. Fluid Mech. 585, 93–115.

* Barthelemy, X., Banner, M., Peirson, W., Fedele, F., Allis, M., Dias, F., 2018. On a unified breaking onset threshold for gravity waves in deep and intermediate depth water. J. Fluid Mech. 841, 463–488.

* Battjes, J.A., Janssen, J.P.F.M., 1978. Energy Loss and Set-Up Due to Breaking of Random Waves. In: Proceedings of the 16th International Conference on Coastal Engineering. ASCE, pp. 569–587.

* Beji, S., Battjes, J.A., 1994. Numerical simulation of nonlinear wave propagation over a bar. Coast Eng. 23, 1–16.

* Belibassakis, K.A., Athanassoulis, G.A., 2011. A coupled-mode system with application to nonlinear water waves propagating in finite water depth and in variable bathymetry regions. Coast. Eng. 58, 337–350.

* Bingham, H.B., Zhang, H., 2007. On the accuracy of finite-difference solutions for nonlinear water waves. J. Eng. Math. 58, 211–228.

* Cienfuegos, R., Barthélemy, E., Bonneton, P., 2010. Wave-breaking model for boussinesq-type equations including roller effects in the mass conservation equation. J. Waterw. Port, Coast. Ocean Eng. 136 (1), 10–26.

* Derakhti, M., Kirby, J.T., 2016. Breaking-onset, energy and momentum flux in unsteady focused wave packets. J. Fluid Mech. 790, 553–581.

* Battjes, J. A., Bakkenes, H. J., Janssen, T. T., & van Dongeren, A. R. (2004). Shoaling of subharmonic gravity waves. *Journal of Geophysical Research: Oceans*, **109**(C2).

* Derakhti, M., Banner, M.L., Kirby, J.T., 2018. Predicting the breaking strength of gravity water waves in deep and intermediate depth. J. Fluid Mech. 848, R2.

* Derakhti, M., Kirby, J.T., Banner, M.L., Grilli, S.T., Thomson, J., 2020. A unified breaking onset criterion for surface gravity water waves in arbitrary depth. Journal of Geophysical Research: Oceans e2019JC015886.

* Filipot, J.-F., Arduhuin, F., 1978–2012. A unified spectral parameterization for wave breaking: from the deep ocean to the surf zone. J. Geophys. Res.: Oceans (C11), 117.

* Gouin, M., Ducrozet, G., Ferrant, P., 2016. Development and validation of a non-linear spectral model for water waves over variable depth. Eur. J. Mech. B Fluid 57, 115–128.

* Grilli, S.T., 1998. Depth inversion in shallow water based on nonlinear properties of shoaling periodic waves. Coast Eng. 35, 185–209.

* Grilli, S.T., Horillo, J., 1997. Numerical generation and absorption of fully nonlinear periodic waves. J. Eng. Mech. 123 (10), 1060–1069.

* Grilli, S.T., Subramanya, R., 1996. Numerical modeling of wave breaking induced by fixed or moving boundaries. Comput. Mech. 17 (6), 374–391.

* Grilli, S.T., Skourup, J., Svendsen, I.A., 1989. An efficient boundary element method for nonlinear water waves. Eng. Anal. Bound. Elem. 6 (2), 97–107.

* Grilli, S., Subramanya, R., Svendsen, I., Veeramony, J., 1994a. Shoaling of solitary waves on plane beaches. J. Waterw. Port, Coast. Ocean Eng. 120 (6), 609–628.

* Grilli, S.T., Losada, M.A., Martin, F., 1994b. Characteristics of solitary wave breaking induced by breakwaters. J. Waterw. Port, Coast. Ocean Eng. 120 (1), 74–92.

* Grilli, S., Svendsen, I., Subramanya, R., 1997. Breaking criterion and characteristics for solitary waves on slopes. J. Waterw. Port, Coast. Ocean Eng. 123 (3), 102–112.

* Grilli, S., Derakhti, M., Kirby, J., Banner, M., 2018. A Unified Formulation for Predicting the Breaking Onset of Gravity Water Waves from Deep to Shallow Water: Validation Cases Using a Fully Nonlinear Potential Flow Model. In: Fall Meeting Abstract. AGU. OS31D-1799.

* Grilli, S.T., Horillo, J., Guignard, S., 2019. Fully nonlinear potential flow simulations of wave shoaling over slopes: spilling breaker model and integral wave properties. Water Waves 1–35.

* Guignard, S., Grilli, S.T., et al., 2001. Modeling of Wave Shoaling in a 2D-NWT Using a Spilling Breaker Model. In: The Eleventh International Offshore and Polar Engineering Conference (ISOPE01, Stavanger, Norway, June 2001), vol. 3. International Society of Offshore and Polar Engineers, pp. 116–123.

* Itay, U., Liberzon, D., 2017. Lagrangian kinematic criterion for the breaking of shoaling waves. J. Phys. Oceanogr. 47 (4), 827–833.

* Kazolea, M., Delis, A.I., Synolakis, C.E., 2014. Numerical treatment of wave breaking on unstructured finite volume approximations for extended boussinesq-type equations. J. Comput. Phys. 271, 281–305.

* Kennedy, A.B., Chen, Q., Kirby, J.T., Dalrymple, R.A., 2000a. Boussinesq modeling of wave transformation, breaking, and runup. i: 1d. J. Waterw. Port, Coast. Ocean Eng. 126 (1), 39–47.
Kennedy, A.B., Chen, Q., Kirby, J.T., Dalrymple, R.A., 2000b. Boussinesq modeling of wave transformation, breaking, and runup. I: 1D. J. Waterway, Port, Coast. and Ocean Eng. 126 (1), 39–47.

Kurnia, R., Van Groesen, E., 2014. High order Hamiltonian water wave models with wave-breaking mechanism. Coast. Eng. 93, 55–70.

Longuet-Higgins, M.S., 1970. Longshore currents generated by obliquely incident sea waves, 1. J. Geophys. Res. 75, 6778–6789.

Longuet-Higgins, M.S., Cokelet, E.D., 1976. The deformation of steep surface waves on water. Part I. a numerical method of computation. Proc. Roy. Soc. Lond. A 350, 1–26.

Longuet-Higgins, M.S., Stewart, R.W., 1964. Radiation stress in water waves, a physical discussion with applications. Deep Sea Res. 11, 529–563.

MacMahan, J.H., Thornton, E.B., Reniers, A.J., 2006. Rip current review. Coast Eng. 53, 191–208.

Madsen, P.A., Sørensen, O., Schäffer, H., 1997. Surf zone dynamics simulated by a boussinesq type model. part i. model description and cross-shore motion of regular waves. Coast. Eng. 32 (4), 255–287.

Miche, A., 1944. Mouvements ondulatoires de la mer en profondeur croissante ou décroissante. forme limite de la houle lors de son déferlement. application aux digues maritimes. Troisième partie. Forme et propriétés des houles limites lors du déferlement. Croissance des vitesses vers la rive. Annales des Ponts et Chaussées Tome 114, 369–406.

Peregrine, D.H., 1983. Breaking waves on beaches. Annu. Rev. Fluid Mech. 15 (1), 149–178.

Perlin, M., He, J., Bernal, L.P., 1996. An experimental study of deep water plunging breakers. Phys. Fluid. 8 (9), 2365–2374.

Qiao, H., Duncan, J.H., 2001. Gentle spilling breakers: crest flow-field evolution. J. Fluid Mech. 439, 57–85.

Raoult, C., Benoit, M., Yates, M., 2016. Validation of a fully nonlinear and dispersive wave model with laboratory non-breaking experiments. Coast Eng. 114, 194–207.

Roeber, V., 2010. Boussinesq-type Model for Nearshore Wave Processes in Fringing Reef Environment, Ph.D. thesis. University of Hawaii at Manoa, Honolulu.

Roeber, V., Cheung, K.F., Kobayashi, M.H., 2010. Shock-capturing boussinesq-type model for nearshore wave processes. Coast. Eng. 57 (4), 407–423.

Saket, A., Peirson, W.L., Banner, M.L., Allis, M.J., 2018. On the influence of wave breaking on the height limits of two-dimensional wave groups propagating in uniform intermediate depth water. Coast. Eng. 133, 159–165.

Schäffer, H.A., Madsen, P.A., Deigaard, R., 1993. A boussinesq model for waves breaking in shallow water. Coast. Eng. 20 (3–4), 185–202.

Seiffert, B.R., Ducrozet, G., Bonnefoy, F., 2017. Simulation of breaking waves using the high-order spectral method with laboratory experiments: wave-breaking onset. Ocean Model. 119, 94–104.

Shi, F., Kirby, J.T., Tehranirad, B., Harris, J.C., Grilli, S., 2011. Funwave-tvd, fully nonlinear boussinesq wave model with tvd solver, documentation and user's manual. Center Appl. Coastal Res., Univ. Delaware, Newark, DE, USA, Res. Rep. 1. No CACR-11-04 2.

Shi, F., Kirby, J.T., Harris, J.C., Geiman, J.D., Grilli, S.T., 2012. A high-order adaptive time-stepping TVD solver for Boussinesq modeling of breaking waves and coastal inundation. Ocean Model. 43, 36–51.

Song, J.-B., Banner, M.L., 2002. On determining the onset and strength of breaking for deep water waves. Part I: unforced irrotational wave groups. J. Phys. Oceanogr. 32, 2541–2558. http://ams.allenpress.com/archive/2541-2558/32/9/pdf/i1520-0485-32-9-2541.pdf.

Sørensen, O.R., Madsen, P.A., Schaffer, H.A., 1998. Nearshore Wave Dynamics Simulated by Boussinesq Type Models. In: Coastal Engineering, vol. 1998, pp. 272–285.

Stansell, P., MacFarlane, C., 2002. Experimental investigation of wave breaking criteria based on wave phase speeds. J. Phys. Oceanogr. 32, 1269–1283. http://ams.allenpress.com/archive/1520-0485/32/5/pdf/i1520-0485-32-5-1269.pdf.

Svendsen, I.A., 1984. Wave heights and set-up in a surf zone. Coast Eng. 8, 303–329.

Symonds, G., Huntley, D.A., Bowen, A.J., 1982. Two-dimensional surf beat: long wave generation by a time-varying breakpoint. J. Geophys. Res.: Oceans 87 (C1), 492–498.

Tanaka, M., 1986. The stability of solitary waves. Phys. Fluid. 29 (3), 650–655.

Thornton, E.B., Guza, R.T., 1983. Transformation of wave height distribution. J. Geophys. Res. 88 (C10), 5925–5938.

Tian, Z., Perlin, M., Choi, W., 2008. Evaluation of a deep-water wave breaking criterion. Phys. Fluids 20 (6), 066604.

Tian, Z., Perlin, M., Choi, W., 2010. Energy dissipation in two-dimensional unsteady plunging breakers and an eddy viscosity model. J. Fluid Mech. 655, 217–257.

Tissier, M., Bonneton, P., Marche, F., Chazel, F., Lannes, D., 2012. A new approach to handle wave breaking in fully non-linear boussinesq models. Coast. Eng. 67, 54–66.

Wei, G., Kirby, J.T., Grilli, S.T., Subramanya, R., et al., 1995. A fully nonlinear boussinesq model for surface waves. part 1. highly nonlinear unsteady waves. J. Fluid Mech. 294 (7), 71–92.

Wu, C.H., Nepf, H.M., 2002. Breaking criteria and energy losses for three-dimensional wave breaking. J. Geophys. Res. 107 (C10), 3177. https://doi.org/10.1029/2001JC001077.

Zakharov, V.E., 1968. Stability of periodic waves of finite amplitude on the surface of a deep fluid. J. Appl. Mech. Tech. Phys. 2, 190–194.
# Chapter 3

# Breaking statistics
parameterizations and observations

## 3.1 Introduction

Being able to quantify the statistics of occurrence of (large) breaking waves in storm seas is a challenge but is required to identify the most dangerous sea states to consider in the ULS design of offshore wind turbines. Breaking statistics are further involved in the parameterization of the energy dissipation in the spectral numerical wave models applied to generate the hindcasts, from which the design conditions are assessed.

In this chapter we shall examine different wave breaking statistics parame-
terizations proposed in the literature and introduce a new one, developed in the
DiMe project. The first section covers the parameterizations for the breaking
probabilities, the second one addresses the models for the averaged spectral crest
length densities from *Phillips* [1985] and the last one introduces the new model
for the probability of occurrence of dominant breaking waves by *Stringari et al.*
[2021b] (DiMe project).

## 3.2 Breaking wave probability parameterizations

3.2.1 Breaking probability parameterizations based on the
Gaussian field theory

Breaking probability is traditionally obtained by computing the ratio of the num-
ber of breaking waves to the total number of waves passing a fixed point over a
statistically representative period of time. It may be a bulk or a spectral quantity
if the analysis attempts to distinguish between the different wave scales that can
be present in a given wave field.
For deep water applications, a first family of wave breaking probability pa-
rameterizations relates to the exploitation of the Gaussian field theory [Rice,
1944; Longuet-Higgins, 1957] and the corresponding wave height-period theoret-
ical joint distributions [e.g. Cavanie et al., 1976].

Ochi and Tsai [1983] developed a spectral breaking probability by integrating, at a given wave period, the joint distribution of Cavanie et al. [1976] for elevations above an arbitrary wave height threshold corresponding to a modified Stokes-like limiting steepness. Srokosz [1986] followed a quite similar path, but relied on the mono-variate distribution of local maxima elevations derived by Cartwright and Longuet-Higgins [1956], combined with a breaking criterion based on the maximal downward acceleration applied to regular waves. They modified Cartwright and Longuet-Higgins [1956]'s distribution to consider only the local maxima with a downward acceleration higher than $\alpha g$. While $\alpha = 0.5$ for regular waves [Longuet-Higgins, 1963], Srokosz [1986] proposed $\alpha = 0.4$ arguing that this value leads to a perfect match with the parameterizations of Ochi and Tsai [1983]. Later, Papadimitrakis [2005] proposed an approach sharing many similarities with the one of Ochi and Tsai [1983] but with a breaking criterion that was wave scale-dependent and a composite function of the local wave steepness, wind forcing and long wave modulation.

In shallow water, one of the most significant effort in that field was brought
by Thornton and Guza [1983] who developed a breaking probability model rely-
ing on the shallow water version of Miche parameter. Their approach inspired
Chawla and Kirby [2002] who studied breaking statistics in presence of current
gradient and later Filipot et al. [2010] who developed a formulation, designed to
capture breaking probabilities of different wave scales and applicable from the
deep ocean to the surf-zone. For each wave scale, they built a breaking wave
height distribution from which the breaking probability can be readily obtained.
Their model was calibrated using deep, intermediate and shallow water data (see
figure 3.2.1 for the deep water validation) and implemented in the WaveWatchIII
model (Chapter 4).

3.2.2 Observation-based breaking probability parameteri-
zations

While the previous work relied on appealing theoretical frameworks, only a few
of them were comforted by validation with breaking waves statistics obtained
in the field. This was probably part of the motivation that triggered a series
of studies and papers by M. Banner, A. Babanin and their colleagues. The
first of this series, Banner et al. [2000], examined a large dataset of dominant
breaking waves, obtained in lakes, closed seas and open ocean conditions. The
main findings of this study was to highlight the correlation of the dominant wave
<figure><img src="image_18.png" /><figcaption>Figure 3.2.1: Comparison between observations presented in *Banner et al.* [2000] (black symbols) and breaking models ($Q_{\text{mod}}$) by *Banner et al.* [2000] (solid black line) and *Filipot et al.* [2010] (magenta triangle for the final version). The variable $\epsilon$ on the x-axis is the dominant wave steepness, given by equation 3.1. Reproduced from *Filipot et al.* [2010].</figcaption></figure>
breaking statistics with the dominant wave steepness:

$$ \epsilon = k_p H_p / 2 \qquad (3.1) $$

where $k_p$ is the dominant wave number and $H_p$ is the dominant wave height obtained by considering the energy density in the interval $[0.7, 1.3]f_p$, with $f_p$ the peak frequency. They further reported a threshold behavior: dominant breaking waves where only observed in wave fields having $\epsilon$ greater than 0.055 and proposed a quadratic dependence of the dominant wave breaking rate as a function of $\epsilon$. This study, was adapted to intermediate water depth by Babanin et al. [2001] who introduced in particular an empirical water depth dependence in Banner et al. [2000]'s formulation. Shortly after, Banner et al. [2002] extended their 2000's effort to wave scales between 1 and 2.48 times the peak frequency, $f_p$. To achieve this, they based their approach on the spectral saturation, noted $\sigma(k) = k^4 E(k)$. Using an innovative methodology they were able to compute the breaking statistics by an iterative decomposition of the sea surface into a range of wave scales (the so-called Riding Wave Removal Method, detailed in Schulz [2009]). They reported again a threshold behavior in terms of $\sigma(k)$, but with a dependence in the wave scales i.e., the shorter the waves, the lower the threshold. Banner et al. [2002] related this dependence to a possible effect of the wave directionality on the breaking statistics and scaled the spectral saturation by the directional spreading, evaluated at each wave scale, which resulted in a common breaking saturation threshold for all wave scales. Since this work, different studies dealing with the modulation of short waves by long waves suggest that the modulation would be another candidate to explain the observed decrease in the spectral saturation breaking threshold for shorter wave scales [see for instance Veras Guimaraes, 2018; Peureux et al., 2021].

## 3.3 Parameterization for the spectral density of breaking crests

Measuring the spectral distribution of the breaking probability causes difficulties illustrated by figure 3.3.1 and well described by Phillips [1985]: "There is clearly some association of the breaking events with waves of different scales, but it is difficult to make the association in an unambiguous way if we consider only the surface configuration at one given instant. A breaking crest may indeed be a local maximum in the instantaneous surface configuration but there is no guarantee that a local wavelength of the breaking wave can be defined clearly."

To circumvent this problem Phillips [1985] proposed the following solution: "It seems more satisfactory to use the velocity $c$ of the breaking front as a measure of the scale of the breaking, since this is a well-defined quantity that might (conceptually at any rate) be measured from ciné images of the sea surface."
These statements led Phillips to introduce a new spectral wave breaking statistics: the distribution $\Lambda(c)$ such as $\Lambda(c)dc$ "represents the average total length per unit surface area of breaking fronts that have velocities in the range $[c, c + dc]$."

This distribution $\Lambda(c)$ can be derived from video images from the ocean surface. Although the practical detection of this quantity comes with challenges, [see for instance *Banner et al.*, 2014], it now represents the basis of the modern wave breaking observations. *Melville and Matusov* [2002] proposed an entirely empirical parameterization as a function of wind speed based on observations of breakers obtained by aerial images. *Makin and Kudryavtsev* [2002b] introduced a model for $\Lambda(c)$ combining Gaussian field theory and empirical findings from *Banner et al.* [2000]. The most widely used parameterization to date, is certainly that of *Sutherland and Melville* [2013]. Their approach relies on a dimensional analysis of the problem to identify the physical variables involved in the breaking process. They concluded that the governing parameters are the ballistic speed of the falling crest, the wave age, the average wave steepness and the non-dimensional fetch.

## 3.4 Introduction of a new wave breaking parameterization based on the Gaussian field theory

### 3.4.1 Main underlying assumptions

Analysing the past parameterizations of breaking statistics, it appears that none of them accounts explicitly and jointly for :

1. The stochastic nature of the waves, that must play a decisive role in the breaking process, especially in deep water, through group modulation effect.

2. The most widely accepted breaking onset criterion, i.e. the kinematic one detailed in the previous chapter.

The DiMe project provided an alternative approach that allows to incorporate the two above-cited missing ingredients into a new wave breaking statistics parameterization presented in *Stringari et al.* [2021b] (see Appendix 3.A).

The new wave breaking parameterization relies on the Gaussian theory-derived joint distribution of the phase velocity $c$ and horizontal fluid velocity $u$, both quantities being estimated at local maxima (or local crests) in the space domain. This formulation $p(u, c)$, derived by M. Prevosto (Ifremer), gives the breaking probability $Q$, by the following integration:

$$Q = \int_{u>Ac} \int_{0}^{\infty} p(u, c)dcdu \qquad (3.2)$$
<figure><img src="image_19.png" /><figcaption>Figure 3.3.1: Method employed by Filipot et al. [2010] to compute the breaking probability of different wave scales from observations collected in lake George, Australia by Babanin et al. [2001]. The first panel shows the total wave elevation in blue, and the black spikes indicate detection of breaking wave-induced bubbles. The four other sub-panels display the 4 different wave scales obtained by a spectral filtering around the frequencies 0.55, 1, 1.86, 3.45 $f_p$ (from top to bottom) and a correction for second-order nonlinear contributions using Forristall [2000]'s method. The red bars indicate the level of nonlinearities while the yellow bars mark the detection of a breaking wave. The main limitations of this approach are somewhat arbitrary definition of the wave scales and the nonlinear correction that was not accurate enough (leading to contamination of the shorter wave scales by peak waves harmonic contributions).</figcaption></figure>
Here *A* is the breaking threshold above which Gaussian (linear) waves are supposed to break, as introduced in Chapter 2. *Stringari et al.* [2021b] presents a validation of this model, focusing on the dominant waves and introducing a breaking criterion adapted to Gaussian (or linear) waves.

### 3.4.2 Introduction of a linear-equivalent breaking onset criterion

The model proposed by *Stringari et al.* [2021b] has only one free parameter, *A* the kinematic breaking threshold for linear waves. We introduced a theoretical threshold based on an equivalence between a nonlinear regular wave at breaking onset and a linear wave of same mechanical energy. Using the results from *Cokelet* [1977], we found *A* = 0.3852 as a linear-equivalent breaking criterion for regular waves. However, the comparison of our model to different observational datasets, suggested that this value of *A* leads to a significant underestimation of the dominant wave breaking statistics. The reasons why it occurs are still unclear, although the effects of waves irregularity, wind forcing or long wave modulations are credible candidates.

Further investigations led us to propose an empirical tuning of the *A* parameter, yielding *A* = 0.24 through a comparison with $\Lambda(c)$ observations obtained in young seas states by *Schwendeman et al.* [2014].

Further efforts are also currently conducted in the DIMPACT project by F. Hulin (PhD student) to explore the existence of a linear-equivalent breaking onset criterion using nonlinear wave simulations (see Chapter 6).

One of the main conclusion of *Stringari et al.* [2021b] paper is that there is a lack of field observations, especially in storm conditions for validating the wave breaking statistics parameterizations. This statement is one of the main motivation of the current field experiment from La Jument lighthouse (see Chapter 6). Another challenge resides in the analysis of sea surfaces video observations to extract meaningful breaking waves statistics. The next section introduces a new methodology, based on machine learning, to detect breaking crests in video records and quantify the $\Lambda(c)$ distribution.

## 3.5 Observations of breaking waves statistics

Capturing breaking statistics in the field has been achieved by visual inspection of video records [*Banner et al., 2000*; *Babanin et al., 2001*; *Banner et al., 2002*], analysis of the sound of breaking induced-bubbles *Manasseh et al.* [2006] and more recently by automatic analysis of video-records to obtain the $\Lambda(c)$ distribution [*Phillips et al., 2001*; *Mironov and Dulov, 2007*; *Sutherland and Melville, 2013*].

One challenge to extract wave breaking statistics from images of the sea surface is to isolate the patches of active foam (corresponding to wave crests actively
44CHAPTER 3. BREAKING STATISTICS MODELS AND OBSERVATIONS

breaking) from the passive or fossil foam left in the wake of breaking waves. Dif-
ferent methods have been implemented to extract the active foam patches and
derive the Λ(c) distributions. These methods [e.g. Mironov and Dulov, 2007;
Kleiss and Melville, 2011] however appear to be very sensitive to sun glitter,
light exposure and camera settings [Schwendeman and Thomson, 2015] making
their generalization to other field observations delicate. Furthermore, Banner
et al. [2014] showed that the estimation of Λ(c) was strongly affected by the
method applied to compute the speed and length of the crests flagged as actively
breaking.

The paper from Stringari et al. [2021a] (see Appendix 3.B) presents the appli-
cation of Artificial Intelligence (AI) methods to extract breaking waves from video
images. It relies on an heavy visual inspection of thousand of images to produce
a relevant dataset to train deep learning neural networks. The machine learning
algorithm shows an improvement of about 50% in the capture of breaking crests
from video records, compared to conventional (non AI-based) approaches. The
other force of this work is to provide the community with an open tool whose
ability to detect breaking waves will increase through the growth of the learning
database, to be expanded by the scientific community. As will be discussed in
Chapter 6, this effort will be used to feed improvements in wave breaking statis-
tics models (discussed in this chapter), spectral wave models dissipation source
terms (Chapter 4) and air-sea coupled models (Chapter 6).

3.6 Discussions & Conclusions

Parameterizing wave breaking statistics is a hard task because of the nonlinear
and stochastic nature of the physical process at hand. The new wave breaking
probability parameterization from Stringari et al. [2021b] carries an interesting
potential but will require more intense validation with field and wave tank ob-
servations. The main limitation is the linear-equivalent breaking threshold that
needs further adjustment. The parallel that is drawn between linear and non
linear waves, underlying our approach, also deserves to be investigated in more
details. Ongoing work to progress on these open questions are further introduced
in Chapter 6.

The lack of more complete observational databases further slows down the
scientific progresses. This lack of reliable information is particularly obvious
in storm conditions for obvious reasons related to instrument survivability at
sea. This motivates the perpetuation of La Jument lighthouse observatory (see
Chapter 6).

While machine learning can bring decisive information to estimate breaking
statistics from video observations [Stringari et al., 2021a], it is also a possible tool
to model breaking statistics. Depending on the performance of our "physical"
wave breaking statistics model, that needs to be further evaluated with additional
observations, deep learning algorithms constitute alternative candidates to help capturing the physical processes not incorporated in our model. This will however, again require large datasets of observations, covering different environmental conditions and a wide spectrum of environmental parameters.
# Appendix

## 3.A Paper from Stringari & al. (2021b)
JGR Oceans

RESEARCH ARTICLE

10.1029/2020JC016943

**Key Points:**

* A new probabilistic wave breaking model based on Gaussian field theory is presented for dominant, wind-sea waves

* Wave breaking probabilities are modeled from the joint probability density between wave phase speed and particle orbital velocity

* The proposed model performs well when compared to six other historical models using three field data sets

**Correspondence to:**

C. E. Stringari and J.-F. Filipot,
Caio.Stringari@france-energies-marines.org;
Jean.Francois.Filipot@france-energies-marines.org

**Citation:**

Stringari, C. E., Prevosto, M., Filipot, J.-F., Leckler, F., & Guimarães, P. V. (2021). A new probabilistic wave breaking model for dominant wind-sea waves based on the Gaussian field theory. *Journal of Geophysical Research: Oceans*, 126, e2020JC016943. https://doi.org/10.1029/2020JC016943

Received 29 OCT 2020
Accepted 11 MAR 2021

# A New Probabilistic Wave Breaking Model for Dominant Wind-Sea Waves Based on the Gaussian Field Theory

C. E. Stringari¹ ID, M. Prevosto², J.-F. Filipot¹, F. Leckler¹, and P. V. Guimarães¹,³

¹France Energies Marines, Plouzané, France, ²Institut Français de Recherche pour l'Exploitation de la Mer, Plouzané, France, ³PPGOceano, Federal University of Santa Catarina, Florianópolis, Brazil

**Abstract** This study presents a novel method for obtaining the probability of wave breaking ($P_b$) of deep water, dominant wind-sea waves (i.e., waves made of the energy within ±30% of the peak wave frequency) derived from Gaussian wave field theory. For a given input wave spectrum, we demonstrate how it is possible to derive a joint probability density function between wave phase speed (c) and horizontal orbital velocity at the wave crest (u) from which a model for $P_b$ can be obtained. A nonlinear kinematic wave breaking criterion consistent with the Gaussian framework is further proposed. Our model would allow, therefore, the application of the classical wave breaking criterion (i.e., wave breaking occurs if $u/c > 1$) in spectral wave models which, to the authors' knowledge, has not been done to date. Our results show that the proposed theoretical model has errors in the same order of magnitude as six other historical models when assessed using three field data sets. With optimization of the proposed model's single free parameter, it can become the best performing model for specific data sets. Although our results are promising, additional, more complete wave breaking data sets collected in the field are needed to comprehensively assess the present model, especially in regards to the dependence on phenomena such as direct wind forcing, long wave modulation, and wave directionality.

**Plain Language Summary** Waves will break if the speed of the water particles on the wave crest is greater than the speed of the wave itself, causing the wave crest to overtake the front part of the wave, leading to wave breaking. Precisely simulating real ocean waves requires, therefore, a particle-by-particle description of the water motion, which is too expensive for current computers to handle in real-world applications. Instead, wave models describe waves by means of their statistical properties, that is, wave properties are averaged over a large number of waves. Here, we present a mathematical formulation that allows us to calculate the combined probability between the speed of particles on the wave crest and the wave speed, based only on statistical properties. From these combined probabilities, we model the probability of wave breaking. Our results indicate that our model performed relatively well when compared to six other models using three historical data sets. Because of a lack of observed data to assess our model, we recommend that future research should focus on collecting more wave breaking data measured in the field. Future advances along this line of research could lead, for example, to improvements on operational weather forecast models.

# 1. Introduction

A robust description of wave breaking is a crucial aspect of wave modeling. It is via wave breaking that most of the wave energy is dissipated and a precise formulation of this phenomenon is required to obtain reliable models. Despite its importance, energy dissipation due to wave breaking is still modeled as a semi-empirical process due to the difficulty of representing physically derived wave breaking criteria on phase-averaged wave models (Ardag & Resio, 2020; Ardhuin et al., 2010; Banner et al., 2000, 2002, 2014; Battjes & Janssen, 1978; Filipot & Ardhuin, 2012; Filipot et al., 2010; Thornton & Guza, 1983; Zieger et al., 2015). The available probabilistic (i.e., parametric or empirical) formulations included in these models have been derived from limited data sets and without rigorous theoretical frameworks and, therefore, currently lack a solid physical background. While the current operational (spectral) models are capable of reproducing field observations of integrated spectral parameters (e.g., significant wave height, peak wave period, and peak wave direction) with good accuracy, it remains unclear if their wave breaking parameterizations are entirely reliable. This knowledge gap partly occurs because limited research has focused on wave breaking statistics derived from field data, especially when it comes to wave breaking observations distributed as
a function of wave scales (e.g., wave frequency or wave phase speed). The research developed here has, therefore, important implications for air-sea flux parameterizations (Kudryavtsev et al., 2014), safety at sea (Kjeldsen et al., 1980), and design of offshore structures (Filipot et al., 2019), all of which directly rely on the properties of breaking waves.

Historically, parametric wave breaking models have been constructed from two different approaches: the first approach considers wave statistics (wave steepness, most frequently) derived from a wave-by-wave analysis of the surface elevation time series collected at a single point location where wave breaking occurrences are synchronously identified (using video data, most frequently). The wave breaking probability (i.e., the ratio between the total number of breaking waves over the total number of waves during a given period of time) can then be expressed as a bulk quantity (Alsina & Baldock, 2007; Chawla & Kirby, 2002; Janssen & Battjes, 2007; Thornton & Guza, 1983) or can be distributed over wave frequency (f), wavenumber (k), or wave speed (c) ranges, referred as to “wave scales” by the wave modeling community (Banner et al., 2002; Eldeberky & Battjes, 1996; Filipot et al., 2010).

The second approach follows from Phillips (1985), who defined the distribution $\Lambda(c)dc$ as the “average total length per unit surface area of breaking fronts that have velocities in the range $c$ to $c + dc$.” This approach, therefore, relates to the analysis of sea surface images in which individual wave breaking patches are tracked in space and time. The main motivation for introducing this new concept was clearly stated in Phillips (1985): “There is clearly some association of the breaking events with waves of different scales, but it is difficult to make the association in an unambiguous way if we consider only the surface configuration at one given instant. A breaking crest may indeed be a local maximum in the instantaneous surface configuration but there is no guarantee that a local wavelength of the breaking wave can be defined clearly. It seems more satisfactory to use the velocity $c$ of the breaking front as a measure of the scale of the breaking.” This quotation clearly identifies the limitations of directly relying on the analysis of single point elevation time series. Different parameterizations have been proposed to quantity $\Lambda(c)dc$ from theoretical (Phillips, 1985) or empirical considerations (Melville & Matusov, 2002; Romero, 2019; Sutherland & Melville, 2013). However, Phillips' (1985) framework remains controversial, particularly regarding its practical application, given that different interpretations of his concepts can generate differences of several orders of magnitude in the calculations of $\Lambda(c)dc$ and its moments (Banner et al., 2014). For a detailed review of commonly used parametric wave breaking models, please refer to Appendix A.

Interestingly, while the ratio between the horizontal orbital velocity at the crest ($u$) to wave phase speed ($c$) appears the most reliable parameter to determine wave breaking occurrence (Barthelemy et al., 2018; Derakhti et al., 2020; Saket et al., 2017; Varing et al., 2020), it was not used by any of the approaches mentioned above. This study provides a new promising wave breaking model by revisiting Rices' (1944) and Longuet-Higgins' (1957) statistical descriptions of Gaussian processes (i.e., for linear waves) to obtain the theoretical joint probability density between $c$ and $u$ ($p(c, u)$). We then model the probability of wave breaking ($P_b$) assuming a kinematic wave breaking criterion consistent with nonlinear waves, that is, a wave breaks if the fluid velocity at the wave crest is greater than the wave phase speed ($u > c$). This study focuses on analyzing dominant waves, defined as waves that have frequencies within $\pm30\%$ of the spectral peak frequency of the wind-sea (Banner et al., 2000). Future research will be dedicated to extend our efforts to broader wave scales. This study is organized as follows: Section 2 describes the proposed model, Section 3 presents three historical data sets used to evaluate the model, Section 4 presents the results, Section 5 discusses, and Section 6 concludes.

## 2. Definition of a Probabilistic Wave Breaking Model Based on Gaussian Field Theory

The kinematic wave breaking criterion $u/c = 1$ has been historically used as the onset of wave breaking for nonlinear, real waves (see Perlin et al. [2013] for a review). Recently, Barthelemy et al. (2018) and Derakhti et al. (2020) confirmed via numerical simulations that waves will inevitably start to break shortly after $u/c$ exceeds 0.85 in deep and shallow water. Further numerical simulations showed that wave breaking occurs when the maximum orbital velocity ($u_{max}$) equals $c$ somewhere along the wave profile and not necessarily at the wave crest (Varing et al., 2020). Although the relationship $u/c$ provides a solid physical background to
establish the onset of wave breaking, this approach has never been applied to spectral wave models because it requires phase-resolving the wave field. In the sections below, we circumvent this difficulty by defining a wave breaking probability model using the joint probability density between c and u corresponding to a given wave energy spectrum (E(f)). The efforts here are consistent with part of the recent work from Ardag and Resio (2020) in the sense that both works aim to solidify the use of the kinematic wave breaking criterion as the standard approach for modeling wave breaking.

## 2.1. Theoretical Derivation of the Joint Probability Density Distribution of Orbital Velocity at the Wave Crest and Phase Speed

Longuet-Higgins (1957) published a very complete work on the statistics of Gaussian wave fields. In particular, Longuet-Higgins (1957) studied the probability density of the speed of zero-crossings along a given line that is of interest in this study. In his study, the speed of zero-crossings were applied in particular to the zero-crossings of the space derivative of a Gaussian process, that is, the velocities of the local maxima in space (Longuet-Higgins, 1957, pp. 356–357). The present study describes how the same methodology can be extended to derive the joint density of the speed of space local maxima (or local crests) and simultaneous wave horizontal orbital velocity for a one-dimensional (1D) Gaussian sea state. For simplicity, this study follows the same notations as those of Longuet-Higgins (1957), and the reader is directed to Section 2.5 in Longuet-Higgins (1957) for further details.

As explained in Longuet-Higgins (1957), if $\xi_1(x,t)$ is a stationary-homogeneous process and we are interested in the points (e.g., in space) where this process crosses a level $x_1$, the joint distribution of the space derivative of $\xi_1$ noted $\xi_2$, with other related processes $\xi_3, \xi_4, ...$ at $\xi_1 = x_1$ is given by:

$$p(\xi_2, \xi_3, \xi_4, \dots)_{x_1} = \frac{|\xi_2| P(\xi_1, \xi_2, \xi_3, \xi_4, \dots)_{\xi_1=x_1}}{N_0(x_1)} \quad (1)$$

where $N_0(x_1)$ is the number of crossings of the level $x_1$ by $\xi_1$ (see Equation 2.2.5 in Longuet-Higgins [1957]). Here, we are interested in joint distributions at the local maxima in space of the wave elevation process $\xi_0$. Therefore, $\xi_1$ is the space derivative of the wave process and local maxima correspond to down-crossings of the zero level by $\xi_1 = \partial\xi_0/\partial x$.

$$\xi_1 = \frac{\partial \xi_0}{\partial x}, \quad \xi_2 = \frac{\partial^2 \xi_0}{\partial x^2} = \frac{\partial \xi_1}{\partial x}. \quad (2)$$

In the case of a Gaussian processes, $N_0(x_1)$ is:

$$N_0^{-}(x_1) = \frac{1}{2\pi} \sqrt{\frac{m_4}{m_2}} \exp\left(-\frac{x_1^2}{2m_2}\right), \quad N_0^{-}(0) = N_0^{-}(0) = \frac{1}{2\pi} \sqrt{\frac{m_4}{m_2}} \quad (3)$$

where $m_0, m_1, ..., m_i$ are the i-th wavenumber spectral moments and the minus sign indicates that we consider only down-crossings.

### 2.1.1. Speed of Local Maxima (Phase Speed)

Following Longuet-Higgins (1957), if we are interested in the speed c of the local maxima in space, that is, the speed of the down-crossings of $\xi_1$, we have:

$$c = -\frac{\partial \xi_1}{\partial t} / \partial x = -\frac{\xi_3}{\xi_2} \quad \text{with} \quad \xi_2 = \frac{\partial^2 \xi_0}{\partial x^2} \quad \text{and} \quad \xi_3 = \partial \xi_1 / \partial t. \quad (4)$$

Using Equation 1,
$$p(\xi_2, \xi_3)_0 = \frac{|\xi_2| p(\xi_1, \xi_2, \xi_3)|_{\xi_1=0}}{N_0^{-}} \qquad (5)$$

with $p(\xi_1, \xi_2, \xi_3)$ the point joint distribution of the three Gaussian processes $\frac{\partial^2 \xi_0}{\partial x}$, $\frac{\partial^2 \xi_0}{\partial x^2}$, $\frac{\partial^2 \xi_0}{\partial x \partial t}$:

$$p(\xi_1, \xi_2, \xi_3) = p(\xi_1) p(\xi_2, \xi_3) = \frac{e^{-\frac{\xi_1^2}{2m_2}} e^{-\frac{1}{2} \left[ \xi_2 \xi_3 \right] Q_c^{-1} \left[ \xi_2 \xi_3 \right]}}{2\pi \sqrt{m_2} \sqrt{(2\pi)^3 \det(Q_c)}} \qquad (6)$$

and covariance matrix:

$$Q = \begin{bmatrix} m_2 & 0 & 0 \\ 0 & m_4 & m_3' \\ 0 & m_3' & m_2'' \end{bmatrix} = \begin{bmatrix} m_2 & 0 \\ 0 & Q_c \end{bmatrix}. \qquad (7)$$

Note that following Longuet-Higgins (1957) notations, $m_i''$ indicates the mixed wavenumber-frequency i-th spectral moment, where the number of quotes indicates the order of the frequency spectral moment, for example,

$$m_3' = \int_0^\infty 2\pi f^1 k^3 E(k) dk, \qquad (8)$$

where $E(k)$ is a given wavenumber spectra.

Classically, to introduce $c$ in the joint density and obtain $p(c, \xi_3)_0$, we apply a change of variables

$$\xi_2 = -\frac{\xi_3}{c}, \quad \xi_3 = \xi_3 \qquad (9)$$

and after the integration of $p(c, \xi_3)_0$ over all the domain of definition of $\xi_3$, we obtain the distribution of $c$ (Longuet-Higgins, 1957, Equation 2.5.19):

$$p(c)_0 = \frac{m_4 m_2'' - m_3'^2}{2 \sqrt{m_4 \left( c^2 m_4 + 2cm_3' + m_2'' \right)^{3/2}}} \qquad (10)$$

Note that the sign on $c$ (or on $m_3'$) depends on the convention on the wave propagation direction. We have kept the convention used by Longuet-Higgins (1957) here.

### 2.1.2. Introducing the Orbital Velocity

As indicated in Equation 1, we can introduce in the formula a variable which represents the horizontal orbital velocity. For Gaussian waves the horizontal orbital velocity $u$ is defined as

$$u = \mathcal{H} \left( \frac{\partial \xi_0}{\partial t} \right) \qquad (11)$$

with $\mathcal{H}$ the Hilbert transform in time domain. Which means that

$$\xi_0 = \sum_i a_i \cos(k_i x - \omega_i t) \qquad (12)$$

is transformed in

$$u = \sum_i a_i \omega_i \cos(k_i x - \omega_i t), \qquad (13)$$
with $a_i$ the wave amplitude, $k_i$ the wavenumber and $\omega_i$ the angular wave frequency of the wave component $i$. As the Hilbert transform is a linear operator, $u$ is also Gaussian. As previously, at the local maxima we have:

$$p(\xi_2, \xi_3, u)_0 = \frac{|\xi_2| p(\xi_1, \xi_2, \xi_3, u)|_{\xi_1=0}}{N_0^{-}} \quad (14)$$

with a new covariance matrix for $\xi_1, \xi_2, \xi_3$, and $u$:

$$Q = \begin{bmatrix} m_2 & 0 & 0 & 0 \\ 0 & m_4 & m'_3 & m'_2 \\ 0 & m'_3 & m''_2 & m''_1 \\ m'_2 & m''_1 & m''_0 \end{bmatrix} = \begin{bmatrix} m_2 & 0 \\ 0 & Q_c \end{bmatrix}. \quad (15)$$

As previously, we can apply a similar change of variables

$$\xi_2 = -\frac{\xi_3}{c}, \quad \xi_3 = \xi_3, \quad u = u, \quad (16)$$

or the easier to deal with,

$$\xi_3 = -c\xi_2, \quad \xi_2 = \xi_2, \quad u = u \quad (17)$$

and integrate $p(c, \xi_2, u)_0$ over all the domain of definition of $\xi_2$. The result is more complicated, but again semi-analytical. The body of the integral has the form

$$e^{-\frac{1}{2} \left[ \xi(c) \xi_2^2 + \beta(c, u) \xi_2 + \alpha(u) \right]} \xi_2^2 \quad (18)$$

and its integration in $\xi_2$ on the down-crossings space $[-\infty, 0]$ gives

$$I(c, u) = \frac{\left( (2\phi^2 + 1)\sqrt{\pi}(\erf(\phi) + 1)e^{\phi^2} + 2\phi \right)}{\sqrt{2}\xi^{3/2}} e^{-\alpha/2} \quad (19)$$

with

$$\phi = \phi(c, u) = \frac{1}{2\sqrt{2}} \frac{\beta(c, u)}{\sqrt{\xi(c)}}, \quad \alpha = \alpha(u), \quad (20)$$

$$\Delta = \det(Q_c) = m'_3(m''_2 m''_1 - m''_3 m''_0) + m_4(m''_0 m''_2 - m''_1^2) + m'_2(m''_3 m''_1 - m''_2 m''_0), \quad (21)$$

$$\alpha(u) = \frac{m_4 m''_2 - m''_3 u^2}{\Delta}, \quad (22)$$

$$\beta(c, u) = 2 \frac{m'_3 m''_1 - m''_2 m''_2}{\Delta} u + 2 \frac{m_4 m''_1 - m''_2 m'_3}{\Delta} u c \quad (23)$$

and

$$\xi(c) = \frac{m''_0 m''_2 - m''_1^2}{\Delta} + 2 \frac{m'_3 m''_0 - m''_2 m'_1}{\Delta} c + \frac{m_4 m''_0 - m''_2^2}{\Delta} c^2. \quad (24)$$

The joint probability density of $(c, u)$ is then:
$$p(c, u) = \frac{1}{N_0^{-1} (2\pi)^2 \sqrt{m_2\Delta}} I(c, u) = \frac{I(c, u)}{2\pi\sqrt{m_4\Delta}}. \qquad (25)$$

Note again that the sign on c and u (or on $m_2'$ and $m_3'$) depends on the convention on the wave propagation direction and Longuet-Higgins (1957)'s convention is still used here. The coefficients $(\alpha, \beta, \xi)$ can be calculated numerically and $\Delta$ is the determinant of $Q_c$, the sub-matrix of $Q$, and after the inverse of $Q_c$ is calculated:

$$Q_c^{-1} = \begin{bmatrix} R & s \\ s' & r \end{bmatrix} \qquad (26)$$

we find

$$\alpha(u) = ru^2, \qquad (27)$$

$$\beta(c, u) = 2[1 \ c]su, \qquad (28)$$

$$\xi(c) = [1 \ c]R[1 \ c]. \qquad (29)$$

An example of the joint density of the couple (phase speed and horizontal particle velocity) at local maxima in space is shown in Figures 1a and 1b for a JONSWAP spectrum.

## 2.2. Modeling $P_b$ From $p(c, u)$

By using Equation 25 applied to the dominant spectral wave band (i.e., that contained in the interval [0.7$f_p$, 1.3$f_p$], where $f_p$ is the peak wave frequency), the probability of dominant wave breaking can be computed by integrating Equation 25 over all phase speeds and for orbital velocities over a threshold $Ac$, with $A$ a constant that will be defined in the next section:

$$P_b = \int_{u>Ac} \int_{0}^{\infty} p(c, u)dcdu. \qquad (30)$$

$P_b$ will be modeled following Equation 30 hereafter. Note that from the definitions in Equation 3, the proposed $P_b$ is defined as the number of breaking local maxima over the total number of local maxima. From the analysis of $p(c, u)$, we observed that spurious, nonmoving local maxima may exist around $c = 0$ and $u = 0$; therefore, to avoid artificially increasing $P_b$, we adopted a practical integration range of $c, u \in [0.05, +\infty]$ here. Note that this range may, however, only be valid for very narrow spectra. Further, we draw attention that, following from Equation 1, our $P_b$ model is defined in the space domain, whereas all the previous $P_b$ models and data are (at least partially) defined in the time domain (see Appendix A for details). For the very narrow spectral band used here, the differences between temporal and spatial definitions of $P_b$ are negligible. This is discussed further in Section 5.

Finally, the proposed model can be extended to accommodate two-dimensional (2D) spectra without changes on how $p(c, u)$ is calculated. This is done by applying an appropriated spreading function to any given 1D spectra (or directly inputting a directional spectra) and by recalculating the moments in Equation 8 to take directionality into account or, more explicitly,

$$m_i = \int_0^{2\pi} \int_0^{\infty} (f \cos\theta \cos\alpha + f \sin\theta \sin\alpha)^i E(f, \theta) df d\theta. \qquad (31)$$

An example considering the simplified cosine spreading law ($D(\theta) = \cos(\theta - \bar{\theta})^{2s}$) with $s = 20$, $\bar{\theta} = 0$ and $\alpha = 0$ applied to same JONSWAP spectrum shown in Figure 1a is shown in Figure 1c. Note that the differences in $p(c, u)$ between the 1D (Figure 1b) and the 2D (Figure 1d) spectra are negligible for the present assumptions. This relatively simple extension allows for the consideration of 2D wave spectral but we caution
<figure><img src="image_20.png" /><figcaption>**Figure 1.** Example of the application of the method. (a) JONSWAP spectrum for $H_{m0}$ = 15 m, $T_p$ = 10 s and shape parameter $\gamma_{js}$ = 10. (b) Obtained joint probability density between the wave phase speed and the horizontal particle velocity at wave crest calculated using Equation 25. Note that the joint probability density was computed using only the spectral energy between $0.7f_p$ and $1.3f_p$, that is, corresponding to the dominant wave band only. (c) Directional spectrum for the same parameters as in (a) and directional spreading $D(\theta) = \cos(\theta - \bar{\theta})^{2s}$ with $s = 20$ and $\bar{\theta} = 0$. (d) Obtained $p(c, u)$ considering only the spectral energy in the direction $\alpha = 0$.</figcaption></figure>

the reader that it may not be fully complete. A follow-up publication will be dedicated to include and assess
the effects of wave directionality in our method more rigorously.

**2.3. Definition of a Gaussian-Equivalent Nonlinear Wave Breaking Criterion**

The previously introduced joint probability density distribution $p(c, u)$ is based on Gaussian theory and therefore assumes that waves are linear. Breaking waves are, however, highly nonlinear. For real nonlinear waves, as detailed in the introduction, it is widely accepted that wave breaking starts when the water particle horizontal velocity at its crest ($u_{nl}$) reaches the wave phase speed ($c_{nl}$). A nonlinear wave breaking criterion can be thus be defined as $A_{nl} = u_{nl}/c_{nl} = 1$. Therefore, we assume that it is possible to obtain an equivalent kinematic criterion, $A_{lin} = \text{constant}$ that relates Gaussian waves to nonlinear waves.

Based on numerical experiments, Cokelet (1977) provided the potential and kinetic energy of a fully non-
linear regular wave in deep-water at the onset of wave breaking (see the last row of his Table A.0). Based on
his results, we define the kinematic criterion as the linear wave that has total energy equals to the nearly
breaking nonlinear regular wave computed by Cokelet (1977). Following Cokelet (1977), where *k*, *g* and *ρ*
are expressed as nondimensional variables, a deep-water wave at the breaking onset (see last row of his Ta-
ble A.0) has kinetic energy *T* = 3.827 × 10⁻² and potential energy *V* = 3.457 × 10⁻². The energy-equivalent
linear wave (denote with subscript *eq*) has, therefore, amplitude:
<figure><img src="image_21.png" /><figcaption>**Figure 2.** a) Example of joint probability density between *u* and *c* obtained from Equation 30. The colored lines indicate different values of *A* and the red dashed line shows *A* = *A*<sub>lin</sub> = 0.382. b) Possible values of *P*<sub>b</sub> for varying *A* calculated using the joint PDF from (a). The vertical dashed line shows *A* = *A*<sub>lin</sub> = 0.382. c) Obtained *P*<sub>b</sub> for varying *H*<sub>m0</sub> and *T*<sub>p</sub> and fixed *A* (0.382) and *γ*<sub>js</sub> (10). The dashed blue lines and marker indicate the *H*<sub>m0</sub> and *T*<sub>p</sub> values used in (a and b). Note that as in Figure 1, these results only consider dominant waves, that is, they were calculated from the spectrum between 0.7*f*<sub>p</sub> and 1.3*f*<sub>p</sub>.</figcaption></figure>

$$a_{eq} = \sqrt{2 \times E} = \sqrt{2 \times (V + T)} = 0.3817. \qquad (32)$$

For this particular case, the linear dispersion relation reads:

$$\omega^2 = gk = 1, \qquad (33)$$

the fluid velocity at crest of the energy-equivalent linear wave is:

$$u_{eq} = \omega a_{eq} = 0.3817, \qquad (34)$$

and the phase speed of the linear wave is:

$$c_{eq} = \sqrt{\frac{g}{k}} = 1. \qquad (35)$$

Given these constants, we obtain:

$$A_{lin} = \frac{u_{eq}}{c_{eq}} = \frac{0.3817}{1} = 0.3817. \qquad (36)$$

Following this approach, we define the correction coefficient $A = A_{lin} = 0.382$ that will be used as reference value hereafter for our tests. This result is consistent with recent findings from Ardag and Resio (2020) who reported from the re-analysis of Duncan's (1981) experimental results, a wave breaking threshold between 0.75 and 1.02 (see their Figure 1). Note, however, that these authors defined their wave breaking threshold as $u/c_g$, where $c_g$ is the group velocity and $u$ was obtained from linear wave theory. Replacing wave group velocity ($c_g$) by the wave phase speed (c) yields a range of possible values between 0.35 and 0.50, which is consistent with $A_{lin}$.

Figure 2 illustrates the sensitivity in wave breaking probability with changes in the wave breaking threshold $A$. For the given $p(c, u)$ in Figure 2a, letting $A$ vary from 0 to 1 resulted in an exponential increase in $P_b$ at $A \le 0.2$ (Figure 2b), which may be unrealistic. When setting $A = A_{lin} = 0.382$ and letting the significant wave height ($H_{m0}$) and wave peak period ($T_p$) vary in the definition of the JONSWAP spectrum, the results indicate that steeper waves are more probable to break, which is expected (Figure 2c). Finally, note that the wave breaking threshold $A$ might be sensitive to other wave and atmospheric parameters such as wave directionality or direct wind forcing (or, equivalently, wave age). In the next sections, the accuracy of our model is assessed using field observations and our results are compared with other parametric wave breaking formulations.
<figure><img src="image_22.png" /><figcaption>**Figure 3.** Field data. (a) Spectral data from TSG14. (b) $\Lambda(c)$ data from STG14. (c) Spectral data from SM13. (d) $\Lambda(c)$ data from SM13. The colored circular markers shown in (a and c) show the peak frequency ($f_p$) and the colored circular markers shown in (b and d) show the peak wave speed ($c_p$). The red dashed line in (b and d) shows the theoretical $c^{-6}$ decay predicted by Phillips (1985). In all plots, the color scale shows the wave age ($c_p/u_*$).</figcaption></figure>

### 3. Field Data

Three historical data sets were used to evaluate the present model. Further, six historical models (detailed in Appendix A) were chosen to contextualize our model in relation to other the state-of-the-art models. These historical models range from baseline models in which the only inputs are known environmental parameters (e.g., wind speed in Melville and Matusov [2002] or wave steepness in Banner et al., 2000) to fairly complex models that account for combinations of several phenomena (e.g., Romero, 2019).

#### 3.1. Thomson (2012) and Schwendeman et al. (2014) Data Set (TSG14)

The first data are from Thomson (2012) and Schwendeman et al. (2014), hereafter TSG14, and were collected in the Strait of Juan de Fuca, Washington. These data were collected by a gray scale video camera with a resolution of 640 × 480 pixels installed above the wheelhouse of Research Vessel R/V *Robertson* which recorded at an acquisition rate of 30 Hz (Schwendeman et al., 2014). These data were then projected into a metric coordinate grid with a resolution of 0.25 m (cross wave) and 0.075 m (along wave) using the method proposed by Holland et al. (1997) and were then used to obtain $\Lambda(c)$ using the spectral approach of Thomson and Jessup (2009). The data were collected in a (usually) fetch-limited region and for a young sea state; note, however, that the particular sea-states analyzed here may not be fetch-limited. Figure 3a shows the measured wave spectra, Figure 3b shows $\Lambda(c)$ distributions, and Table 1 shows a summary of these data. For these data, $P_b$ was calculated using the measured $\Lambda(c)$ distributions combined with the method described below in Equation 37. Additional information regarding the data collection is available from Thomson (2012) and Schwendeman et al. (2014).
<table><caption>**Table 1**

Data Summary for the Three Experiments Described in Sections 3.1, 3.2, and 3.3.</caption><thead><tr><th>Data set</th><th>Date [-]</th><th>Length [min]</th><th>H<sub>m0</sub> [m]</th><th>T<sub>p</sub> [s]</th><th>H<sub>p</sub> [m]</th><th>ε [-]</th><th>U<sub>10</sub> [ms<sup>-1</sup>]</th><th>u* [ms<sup>-1</sup>]</th><th>c<sub>p</sub> [ms<sup>-1</sup>]</th><th>Wave age [-]</th><th>P<sub>b</sub> [-]</th></tr></thead><tbody><tr><td>TSG14</td><td>14/02/2011 20:33</td><td>6.5</td><td>0.75</td><td>2.88</td><td>0.66</td><td>0.160</td><td>11.50</td><td>0.373</td><td>4.50</td><td>12.07</td><td>3.54E-03</td></tr><tr><td>TSG14</td><td>14/02/2011 20:58</td><td>5.1</td><td>0.75</td><td>2.96</td><td>0.66</td><td>0.152</td><td>12.55</td><td>0.417</td><td>4.62</td><td>11.08</td><td>9.57E-03</td></tr><tr><td>TSG14</td><td>14/02/2011 21:30</td><td>6.5</td><td>0.91</td><td>2.99</td><td>0.82</td><td>0.184</td><td>15.07</td><td>0.561</td><td>4.67</td><td>8.33</td><td>6.29E-02</td></tr><tr><td>TSG14</td><td>14/02/2011 21:44</td><td>8.5</td><td>1.09</td><td>3.17</td><td>1.00</td><td>0.200</td><td>15.73</td><td>0.599</td><td>4.94</td><td>8.25</td><td>1.01E-01</td></tr><tr><td>TSG14</td><td>14/02/2011 22:29</td><td>6</td><td>1.21</td><td>3.44</td><td>1.09</td><td>0.186</td><td>17.24</td><td>0.636</td><td>5.36</td><td>8.44</td><td>1.51E-01</td></tr><tr><td>TSG14</td><td>14/02/2011 22:37</td><td>4.8</td><td>1.37</td><td>3.53</td><td>1.24</td><td>0.199</td><td>18.01</td><td>0.660</td><td>5.52</td><td>8.36</td><td>7.61E-02</td></tr><tr><td>TSG14</td><td>15/02/2011 19:04</td><td>10</td><td>0.87</td><td>3.29</td><td>0.79</td><td>0.146</td><td>14.45</td><td>0.360</td><td>5.13</td><td>14.28</td><td>3.75E-03</td></tr><tr><td>TSG14</td><td>15/02/2011 19:19</td><td>6</td><td>0.90</td><td>3.31</td><td>0.81</td><td>0.149</td><td>13.11</td><td>0.477</td><td>5.17</td><td>10.85</td><td>4.05E-02</td></tr><tr><td>SM13</td><td>06/12/2010 21:59</td><td>10</td><td>0.61</td><td>3.51</td><td>0.52</td><td>0.085</td><td>6.46</td><td>0.205</td><td>5.48</td><td>26.68</td><td>7.96E-03</td></tr><tr><td>SM13</td><td>06/12/2010 23:00</td><td>10</td><td>0.61</td><td>3.33</td><td>0.54</td><td>0.097</td><td>7.55</td><td>0.342</td><td>5.20</td><td>15.22</td><td>1.95E-03</td></tr><tr><td>SM13</td><td>07/12/2010 00:00</td><td>10</td><td>0.73</td><td>3.45</td><td>0.66</td><td>0.112</td><td>8.62</td><td>0.319</td><td>5.38</td><td>16.85</td><td>3.24E-03</td></tr><tr><td>SM13</td><td>08/12/2010 00:00</td><td>10</td><td>0.34</td><td>2.04</td><td>0.23</td><td>0.110</td><td>5.24</td><td>0.160</td><td>3.19</td><td>19.96</td><td>1.65E-02</td></tr><tr><td>B00 (SO)</td><td>10/06/1992</td><td>5</td><td>9.20</td><td>13.46</td><td>8.02</td><td>0.089</td><td>19.80</td><td>0.835</td><td>21.01</td><td>25.17</td><td>2.70E-02</td></tr><tr><td>B00 (SO)</td><td>11/06/1992</td><td>9</td><td>4.20</td><td>12.04</td><td>3.66</td><td>0.051</td><td>16.00</td><td>0.626</td><td>18.78</td><td>30.02</td><td>0.00 E+00</td></tr><tr><td>B00 (BS)</td><td>1993</td><td>34-68</td><td>0.39</td><td>2.78</td><td>0.34</td><td>0.089</td><td>11.70</td><td>0.414</td><td>4.34</td><td>10.49</td><td>3.80E-02</td></tr><tr><td>B00 (BS)</td><td>1993</td><td>34-68</td><td>0.49</td><td>2.94</td><td>0.43</td><td>0.100</td><td>12.70</td><td>0.461</td><td>4.59</td><td>9.96</td><td>6.50E-02</td></tr><tr><td>B00 (BS)</td><td>1993</td><td>34-68</td><td>0.53</td><td>3.33</td><td>0.47</td><td>0.084</td><td>14.00</td><td>0.524</td><td>5.20</td><td>9.93</td><td>6.00E-02</td></tr><tr><td>B00 (BS)</td><td>1993</td><td>34-68</td><td>0.54</td><td>3.23</td><td>0.47</td><td>0.092</td><td>14.40</td><td>0.544</td><td>5.04</td><td>9.26</td><td>5.20E-02</td></tr><tr><td>B00 (BS)</td><td>1993</td><td>34-68</td><td>0.38</td><td>2.27</td><td>0.34</td><td>0.131</td><td>15.00</td><td>0.574</td><td>3.55</td><td>6.18</td><td>6.30E-02</td></tr><tr><td>B00 (BS)</td><td>1993</td><td>34-68</td><td>0.45</td><td>2.56</td><td>0.40</td><td>0.121</td><td>14.60</td><td>0.554</td><td>4.00</td><td>7.23</td><td>6.70E-02</td></tr><tr><td>B00 (BS)</td><td>1993</td><td>34-68</td><td>0.45</td><td>2.44</td><td>0.40</td><td>0.134</td><td>13.70</td><td>0.509</td><td>3.81</td><td>7.49</td><td>8.40E-02</td></tr><tr><td>B00 (BS)</td><td>1993</td><td>34-68</td><td>1.19</td><td>5.88</td><td>1.04</td><td>0.061</td><td>8.70</td><td>0.295</td><td>9.18</td><td>31.10</td><td>0.00 E+00</td></tr><tr><td>B00 (BS)</td><td>1993</td><td>34-68</td><td>1.32</td><td>6.24</td><td>1.15</td><td>0.060</td><td>11.20</td><td>0.391</td><td>9.74</td><td>24.91</td><td>0.00 E+00</td></tr><tr><td>B00 (BS)</td><td>1993</td><td>34-68</td><td>0.83</td><td>6.24</td><td>0.73</td><td>0.038</td><td>9.50</td><td>0.322</td><td>9.74</td><td>30.22</td><td>0.00 E+00</td></tr><tr><td>B00 (BS)</td><td>1993</td><td>34-68</td><td>0.89</td><td>5.88</td><td>0.78</td><td>0.045</td><td>10.70</td><td>0.368</td><td>9.18</td><td>24.91</td><td>0.00 E+00</td></tr><tr><td>B00 (BS)</td><td>1993</td><td>34-68</td><td>0.99</td><td>3.71</td><td>0.87</td><td>0.127</td><td>10.00</td><td>0.339</td><td>5.79</td><td>17.06</td><td>3.40E-02</td></tr><tr><td>B00 (BS)</td><td>1993</td><td>34-68</td><td>0.88</td><td>4.00</td><td>0.77</td><td>0.097</td><td>8.70</td><td>0.295</td><td>6.24</td><td>21.14</td><td>5.80E-02</td></tr></tbody></table>

Note. The parameters obtained from wave spectra were computed specifically for the dominant wave band shown in Figure 3 for TSG14 and SM13 cases. The wave height (Hp) and wave steepness (ε) parameters for dominant waves were calculated as per Banner et al. (2002) (see Section A1 for details). The wave age parameter was calculated as cp/u*.

## 3.2. Sutherland and Melville (2013) Data Set (SM13)

The second data set is from Sutherland and Melville (2013), hereafter SM13, and was collected using the Research Platform R/P FLIP during a 2-day field campaign in the Southern California Bight under the scope of the SoCal 2010 experiment (Sutherland & Melville, 2013). Here, we focus only on the visible im- agery collected by these authors to keep consistency with the previously presented data. Stereo video data were collected by a pair of video cameras mounted on the R/P FLIP for 10 minutes at the start of each hour and Λ(c) was obtained using a variation of the method of Kleiss and Melville (2011), that is, tracking the temporal evolution of breakers obtained via pixel intensity threshold. Figure 3c shows the measured wave spectra, Figure 3d shows Λ(c) distributions, and Table 1 shows a summary of these data. Note that because wave breaking was not observed for frequencies below 0.2 Hz and from numerical simulations (not shown) these waves corresponded to a cross-swell not forced by the wind, our analyses only consider waves in the frequency range 0.2 < f < 0.8 Hz. Additional information regarding the data collection is available from Sutherland and Melville (2013). For these and TSG14 data, P<sub>b</sub> was calculated using the measured Λ(c) dis- tributions combined with the formulas from Banner and Morison (2010):
$$P_b = \frac{\int_{c_0}^{c_1} c\Lambda(c)dc}{\int_{c_0}^{c_1} c\Pi(c)dc} \qquad (37)$$

where $c_0 = \frac{g}{2\pi} \frac{1}{1.3f_p}$, $c_1 = \frac{g}{2\pi} \frac{1}{0.7f_p}$, $\Pi(c) = \chi g/(2\pi c^3)$ and $\chi = 0.6$. The implication of this choice is discussed in further detail in Section 5.

### 3.3. Banner et al. (2000) Data Set (B00)

The third data set is from Banner et al. (2000), hereafter B00, and was collected in the Black Sea (BS), Lake Washington (LW), and the Southern Ocean (SO). These authors directly provide values for significant wave height $H_{m0}$, peak period ($T_p$), and the wave breaking probability in their Tables 1 (Black Sea, denoted as BS here) and 2 (Southern Ocean, denoted as SO here). The majority of the data were collected in the BS (13 data runs) and two data runs are from the SO. Given that the original spectral data were not published alongside their study, we approximate the observed spectra using the provided pairs $H_{m0}$, $T_p$ assuming a JONSWAP shape with $\gamma_{js} = 3.3$ (e.g., as previously done in Filipot et al. [2010]). Given that here we are only interested in a very narrow spectral band, the differences between observed and simulated spectra should be minimal. For more details regarding this data, refer to Banner et al. (2000).

## 4. Results

### 4.1. Comparison With Field Data

Figure 4 shows the comparison between estimated (or observed) (x-axis) and modeled (y-axis) values of $P_b$ for each model. In general, no model was able to closely reproduce the trends seen in the combined observed data, regardless of the underlying mathematical or physical formalism. Furthermore, orders of magnitude of difference between the models and, more worryingly, between the models and the measured data were observed. In general, models based on a wave steepness-derived wave breaking criterion (e.g., Banner et al., 2000, 2002) overestimated data derived from $\Lambda(c)$ while models based on $\Lambda(c)$ (e.g., Melville and Matusov, 2002; Sutherland and Melville, 2013) underestimated $P_b$ data that was not derived from $\Lambda(c)$ (i.e., B00 data). The model from Filipot et al. (2010) was found to be the most consistent model. From Figure 4g, the formulation presented in Section 2 with $A = A_{lin} = 0.382$ underestimated the observed $P_b$ for B00 and SM13 data (note that $P_b$ was too low to be displayed in the plot) but performed relatively well for the majority of TSG14 data. Using the mean absolute error (MAE) as a convenient metric to assess the models, it was found that the present model has errors in the same order of magnitude as the previous models. Given the spread in the results seen in Figure 4, no model could be considered a clear winner. For the discussion of these results, see Section 5.

### 4.2. Model Optimization

From the analysis of Figure 2, minor changes in $A$ can lead to major variations in $P_b$. Further, from the analysis of Figure 4, the proposed model underestimated $P_b$ for $A = A_{lin} = 0.382$, particularly for S13 and B00 data. Given that it is a common practice to optimize wave breaking models for particular data sets, we present two methods to do so using TSG14 data as an example. The same could be done for B00 and SM13 data but, for brevity, this is not done here. Given that the present model is not computationally expensive, the first approach consisted of varying $A$ from 0.1 to 0.5 in 0.001 intervals and finding the value of $A$ that resulted in the lowest squared error ($\sqrt{(p_{bi}^d - p_{bi}^m)^2}$, where the superscripts $d$ and $m$ indicate observed and modeled data, respectively) for each data run. Figure 5a shows the results of this procedure. The value $A = A_{opt} = 0.24$ was, on average, the optimal values for this particular data set. The second approach consisted of parameterizing the optimal value of $A$ for each data run as a function of a known environmental variable, in this example, the wave age $c_p/u^*$. (Figure 5b). The results of these two approaches are shown in Figures 5c and 5d, respectively. Both approaches considerably improved the model results from the baseline
<figure><img src="image_23.png" /><figcaption>**Figure 4.** Compassion between measured and computed $P_b$ for different models and data. (a) Banner et al. (2000), (b) Banner et al. (2002), (c) Melville and Matusov (2002), (d) Filipot et al. (2010), (e) Sutherland and Melville (2013), (f) Romero (2019), and (g) present model with $A = A_{lin} = 0.382$. The thick black line shows the linear regression between measured and modeled $P_b$ and the blue dashed line indicates the one-to-one correspondence in all panels. Data points with modeled $P_b < 10^{-5}$ or observed $P_b = 0$ are not shown in this plot. In all plots, $r_{xy}$ is Pearson's correlation coefficient and MAE indicates the MAE. Note. The logarithmic scale.</figcaption></figure>

model presented in Figure 4, with the parametric model (Figure 5d) performing slightly better when considering Pearson's correlation coefficient ($r_{xy}$) as a comparison metric.

**5. Discussion**

We have introduced a new model for obtaining the probability of wave breaking ($P_b$) for dominant waves based on the theoretical joint probability density distribution between wave phase speed (c) and horizontal orbital velocity at the wave crest (u) for unidirectional Gaussian wave fields. The present model has only one parameter for defining the wave breaking threshold (A), which makes it relatively easy to optimize for a given data set (as shown in Section 4.2). While the proposed model performed relatively well for one of
<figure><img src="image_24.png" /><figcaption>**Figure 5.** Results of the optimization procedures. (a) Optimization curves for each data record (colored lines) and the global averaged (black line). The vertical dashed line shows $A = A_{opt} = 0.24$. (b) Parametrization of $A$ as a function of $c_p/u$. The blue swath indicates the 95% confidence interval. For this particular case, $A = 0.008c_p/u* + 0.16$. Note. The logarithmic scale in (a, c, and d). In all plots, the color scale shows the wave age ($c_p/u*$). In (b–d) $r_{xy}$ is Pearson's correlation coefficient.</figcaption></figure>

the investigated data sets (TSG14), it greatly underestimated $P_b$ for the two other data sets (SM13 and B00). For the data investigated here, such underestimation did not result in a high MAE and, in fact, our model had one of the lowest MAE. Recent results of Barthelemy et al. (2018), Derakhti et al. (2020), and Varing et al. (2020) showed that waves with the horizontal fluid velocity that exceeds 0.85 times the phase velocity will inevitably break. These results suggest that the breaking threshold derived from Cokelet (1977) in Section 2.3 could be reduced by ≈15%. If we apply their findings to our case, we obtain $A = 0.382 \times 0.85 = 0.324$ which would help to reduce the underestimation of $P_b$, but not significantly. It is more probable that other environmental phenomena such as direct wind forcing, directional spreading, and long wave modulation, which are not accounted in our model, are the reason for such differences.

One of the most challenging aspects when assessing our model is, nevertheless, regarding the field data. The attribution of wave breaking occurrences to wave scales using time-series analysis, as done in Banner et al. (2000) or Filipot et al. (2010), is difficult because several wave scales can be present at the same time and space. This lead us to use $\Lambda(c)$ observations as well as data from Banner et al. (2000) to investigate our model. Different interpretations of how $\Lambda(c)dc$ is computed from field data can, however, generate orders of magnitude of difference in its moments (Banner et al., 2014; Gemmrich et al., 2013) and, consequently, in $P_b$. Next, it is difficult to relate the speed of the wave breaking front to the phase speed of the carrying wave because small, slower breaking waves could merely be traveling on top of longer, much faster waves. In particular, we believe that these wave breaking events can significantly contribute to the observed $\Lambda(c)dc$ distribution as they would have $c$ close to the peak wave phase speed (i.e., artificially increasing $P_b$). This wave breaking “sub-population” has not received much research interest because of its apparent small contribution to energy dissipation but, for our particular case, they directly impact model validation.
Further, relating $\Lambda(c)$ to $P_b$ is also challenging. Here, we adopted the convenient formula from Banner and Morison (2010). While this formula has some support from the literature (Ardhuin et al., 2010), the actual functional form of $\Pi(c)$ and the value for the constant $\chi$ (see Equation 37) are unknown and changes in these will lead to changes in $P_b$. The Gaussian framework developed in Section 2.1 provides an alternative method to obtain $\Pi(c)$ (e.g., from Equation 3), but this is beyond the scope of this introductory study and will be the focus of a future publication.

Finally, we would like to re-emphasize that our model is derived in the space domain, whereas $P_b$ data is (at least partially) obtained in the time domain. For the narrow spectral band investigated here, Monte Carlo simulations of linear waves indicate that the difference between $P_b$ modeled in space is less than 5% from $P_b$ modeled in time (not shown). Given all these complications and the fact that some historical models are being compared to data that was used to create them (e.g., Banner et al., 2000; Sutherland and Melville, 2013), we are unable to provide an accurate ranking of the existing models. Future research should focus, therefore, on obtaining $P_b$ data that is unambiguous and widely available. In this regard, and despite its own limitations, wave tank experiments could bring further insight into the statistics of dominant (or not) breaking waves. Such a data set would ultimately allow researchers to focus on models derived from physical and mathematical concepts (such as ours) rather than on empirical concepts.

## 6. Conclusion

We have presented a new statistical wave breaking model derived from Gaussian field theory that we have applied to obtain the probability of wave breaking for dominant, wind-sea waves. Although more mathematically complex than previous formulations, the present model relies on the ratio between the crest orbital velocity and the wave phase speed and uses only on a single free parameter, the wave breaking threshold A. Using theoretical results obtained by Cokelet (1977) for regular nearly breaking waves, we derived a wave breaking threshold to adapt our linear model to nonlinear waves. The present model has errors in the same order of magnitude as six other historical models when assessed using three field data sets. For a particular data set (TSG14), our model performed well, especially if the free-parameter A is fine-tuned. Additional observations are, however required, to further understand and quantify the dependence of A on environmental parameters that are not accounted for in our model (e.g., wind forcing, wave directionality, or modulation by long waves). Future research should be dedicated to collect more wave breaking observations in different and repeatable environmental conditions to provide reliable constraints for the optimization of the present and other wave breaking models. Still and although the research presented here is in early stages, the present model should be extendable to waves of any scale and, therefore, has the potential to be implemented in current state-of-the-art spectral wave models as a new wave breaking dissipation source term with relatively little effort.

## Appendix: Historic Parametric Wave Breaking Models

### A1. Banner et al. (2000)

Banner et al.'s (2000) is a popular model for calculating wave breaking probabilities for deep water, dominant waves. This model follows from observations and results from Donelan et al. (1972), Holthuijsen and Herbers (1986), and Banner and Tian (1998), who demonstrated the importance of the wave group modulation on the wave breaking onset. These authors conveniently obtained a parameterization for the probability of wave breaking ($P_b$) based solely on the spectral steepness of the dominant wave scale ($\epsilon_p$), assuming that their formulas would capture the influence of the wave group modulation on the wave breaking onset. Their formulation was derived using a data set of measurements collected in various environments ranging from lakes to open ocean conditions (Banner et al., 2000). From these observations, these authors were then able to obtain a wave breaking threshold behavior for the dominant waves as a function of the dominant spectral wave steepness given by:

$$ \epsilon_p = \frac{H_p k_p}{2} \qquad (A1) $$
in which $k_p$ is the wavenumber at peak frequency ($f_p$) and $H_p$ is the significant wave height of the dominant waves calculated as:

$$H_p = 4\sqrt{\int_{0.7f_p}^{1.3f_p} E(f)df} \qquad (A2)$$

where $E(f)$ is the wave energy spectra as a function of frequency. For their data, $P_b$ was then parameterized as a single equation with three free parameters ($p_1, p_2, p_3$):

$$P_b = p_1 + (c_p - p_2)^{p_3}, \qquad (A3)$$

For the available field data, Banner et al. (2000) found optimal values of $p_1 = 22$, $p_2 = 0.055$, and $p_3 = 2.01$. Note that hereafter, free parameters for the different models will be denoted as $p_n$ where $n$ is a sequential number.

## A2. Banner et al. (2002)

This study extended Banner et al. (2000) model to shorter wave scales (up to 2.48 times the peak wave frequency). From field data Banner et al. (2002) reported that the waves were breaking if the saturation spectrum $\sigma(f) = 2\pi^4 f^3 E(f)/2 g^2 = \sigma(k) = k^4 E(k)$ exceeded a threshold that was frequency dependent. These author's related this dependence to the directional spreading $\theta(k)$ which later led Banner and Morison (2010) to explicitly define the following empirical formulation:

$$P_b(k_c) = \mathcal{H}_h(\tilde{\sigma}(k_c) - p_1) \times p_2 \times (\tilde{\sigma}(k_c) - \tilde{\sigma}_t), \qquad (A4)$$

in which $\mathcal{H}_h$ is the Heaviside step function, $k_c$ is the central wavenumber for a given wavenumber range, $\tilde{\sigma}(k_c) = \sigma(k_c) / \theta(k_c)$ is the saturation spectrum normalized by the averaged directional spreading, $p_1 = 0.0045$ and $p_2 = 33$ are constants obtained from their observations. Following Banner et al. (2002), the directional spreading angle is calculated according to Hwang et al. (2000) (their Equation 19a):

$$\theta\left(\frac{k}{k_p}\right) = \begin{cases} 0.35 + 1.05\left(1 - \frac{k}{k_p}\right) & \text{if } \frac{k}{k_p} < 1.05 \\ 0.30 + 0.087\left(\frac{k}{k_p} - 1\right) & \text{if } 1.05 \le \frac{k}{k_p} < 5 \end{cases} \qquad (A5)$$

where $\theta$ is the directional spreading angle as a function of the wavenumber.

## A3. Filipot et al. (2010)

This method follows from the original studies of Le Mehaute (1962), Battjes and Janssen (1978), and Thornton and Guza (1983) and assumes that the probability distribution function (PDF) of breaking wave heights in the dominant wave scale is parameterized by its central frequency $f_c$ or, equivalently, by its representative phase speed $c(f_c)$ and the product between a Rayleigh PDF for the wave heights

$$P(H, f_c) = \frac{2H}{H_{rms}^2(f_c)} \exp\left[-\left(\frac{H}{H_{rms}(f_c)}\right)^2\right] \qquad (A6)$$

in which

$$H_r(f_c) = \frac{4}{\sqrt{2}} \sqrt{\int_0^\infty U_{fc}(f)E(f)df} \qquad (A7)$$

and
$$U_{fc} = 0.5 - 0.5\cos\left(\frac{\pi}{\delta}\left[\frac{f}{f_c} - 1 - \delta\right]\right) \quad (A8)$$

where $\delta$ is the bandwidth of a Hann window (in this study, $\delta = 0.6$), and a weighting function

$$W(H, f_c) = p_1 \left[ \frac{\beta_r}{\beta} \right]^2 \left\{ 1 - \exp \left[ - \left( \frac{\beta}{\bar{\beta}} \right)^{p_2} \right] \right\} \quad (A9)$$

in which $\beta = kH/\tanh(kh)$, and $p_1$ and $p_2$ are free parameters. To extend the formulation outside the shallow water domain, these authors replaced Thornton and Guza's (1983) breaking criterion based on the wave height ($H$) to water depth ($h$) ratio ($\gamma = H/h = 0.42$) with an adaptation of Miche's (1944) wave breaking parameter:

$$\beta_r = \frac{\bar{k}_r(f_c) H_r(f_c)}{\tanh(\bar{k}_r(f_c)h)} \quad (A10)$$

in which

$$\bar{k}_r(f_c) = \frac{\int_0^\infty U_{fc}(f)k(f)E(f)df}{\int_0^\infty U_{fc}(f)E(f)df} \quad (A11)$$

and

$$\tilde{\beta} = b(b_3 \tanh(kh)^3 - b_2 \tanh(kh)^2 + b_1 \tanh(kh) - b_0) \quad (A12)$$

in which $b = 0.48$, $b_3 = 1.0314$, $b_2 = 1.9958$, $b_1 = 1.5522$, and $b_0 = 0.1885$. In their model, the variable $\tilde{\beta}$ was obtained via numerical calculations of regular nearly breaking waves using the stream wave theory of Dean (1965). Finally, the wave breaking probability is obtained as:

$$P_b(f_c) = \int_0^\infty P(H, f_c)W(H, f_c) dH \le 1. \quad (A13)$$

To keep consistency with Section A2, $P_b$ will be only considered at the spectral peak; other definitions are, however, also possible.

### A4. Models Based on Phillips' (1985) $\Lambda(c)$

The major issue with the previous models is the difficulty to obtain reliable observations of the wave breaking probabilities as a spectral distribution solely from point measurements. Due to the presence of different wave scales at the time and location, it is indeed difficult to assign the breaking occurrence to a given wave frequency of wavenumber. To avoid this problem, Phillips (1985) proposed to use the speed of the breaking front as a proxy for the phase speed of the carrying wave. Phillips (1985) defined the parameter $\Lambda(c)dc$ as the “average total length per unit surface area of breaking fronts that have velocities in the range $c$ to $c + dc$” and then defined the following quantities:

$$L = \int \Lambda(c)dc \quad (A14)$$

and

$$R = \int c\Lambda(c)dc \quad (A15)$$

which represent the “total length of breaking fronts per unit area” (Equation A14) and “the total number of breaking waves of all scales passing a given point per unit time” (Equation A15). Assuming that Phillips (1985) assumptions hold, it is possible to obtain parametric models for $\Lambda$ from known variables (e.g., wind speed) and, consequently, for $P_b$ (see Equation 37).
### A4.1. Melville and Matusov (2002)

Melville and Matusov's (2002) model for $\Lambda(c)$ relies only on the wind speed measured at 10 m ($U_{10}$) to obtain $\Lambda(c)$. Following Melville and Matusov (2002) and using the explicit formula given by Reul and Chapron (2003), this parameterization is written as:

$$ \Lambda(c) = p_1 \left[ \frac{U_{10}}{10} \right]^3 10^{-4} \exp[-(p_2 c)] \quad (A16) $$

in which $p_1$ and $p_2$ are constants. For their data, Melville and Matusov (2002) found $p_1 = 3.3$ and $p_2 = 0.64$. As discussed by Reul and Chapron (2003), this formulation approaches Phillips's (1985) theoretical $c^{-6}$ but may overly estimate the amount of small breakers.

### A4.2. Sutherland and Melville (2013)

Sutherland and Melville (2013) used dimensional analysis to scale $\Lambda(c)$ and obtain a parameterization that is a function of the wind drag ($u_*$), peak wave phase speed ($c_p$), significant wave height ($H_s$), and three constants. From Sutherland and Melville's (2013) Equation 9 and their Figure 4, $\Lambda(c)$ is calculated as:

$$ \Lambda(c) = p_1 \frac{g}{c_p^3} \left(\frac{u_*}{c_p}\right)^{p_2} \left(\frac{c}{\sqrt{gH_s}}\left(\frac{gH_s}{c_p^2}\right)^{p_3}\right)^{-p_6} \quad (A17) $$

where $p_1 = 0.05$, $p_2 = 0.5$, and $p_3 = 0.1$ are constants obtained from the available data. Their formulation reproduces Phillips's (1985) $c^{-6}$ frequency dependency but does not have the typical roll-off at low $c$ as these authors chose to use infrared (other than visible) imagery to obtain and model their $\Lambda(c)$. This choice included the contribution of micro-scale breakers that do generate visible bubbles in their model, hence the difference.

### A4.3 Romero (2019)

Recently, Romero (2019) developed and implemented a new wave breaking parameterization in Wave-WatchIII which relies exclusively on $\Lambda(c)$. Different from previous parameterizations, Romero's (2019) takes into account both the modulations due to winds and long waves on $\Lambda(c)$. His model is fairly general but depends on six free parameters that needed to be laboriously obtained by comparing WaveWatchIII's significant wave height outputs with available measured significant wave heights from buoy data. In Romero's (2019) model, $\Lambda$ was modeled assuming that it is proportional to the crest lengths exceeding a slope threshold:

$$ \Lambda(f, \theta) = \left( \frac{2(2\pi)^2 p_1}{g} \right) f \exp\left[ -\left( \frac{p_2}{B(f, \theta)} \right) \right] M_{LW} M_W \quad (A18) $$

where $p_1 = 3.5 \times 10^{-5}$ and $p_2 = 5 \times 10^{-3}$ are constants to be obtained from the data, $M_{LW}$ is the modulation due to long waves, $M_W$ is the modulation due to winds and $B(f, \theta)$ is the directional wave breaking saturation spectra:

$$ B(f) = \int_{0}^{2\pi} B(f, \theta) d\theta = E(f) \left( \frac{2\pi f^5}{2g} \right). \quad (A19) $$

The modulation due to long waves is calculated according to Guimarães (2018):

$$ M_{LW} = \left[ 1 + p_3 \sqrt{\text{cmss}(E(f))} \cos^2(\theta - \hat{\theta}) \right]^{p_4} \quad (A20) $$

where $p_3 = 400$ and $p_4 = 3/2$ are also best-fit constants found by Romero (2019). The cumulative mean square slope (cmss) is defined as:
$$ \text{cmss} = \int_{0}^{\infty} E(f) \left( \frac{(2\pi)^4 f^4}{g^2} \right) df. \qquad (A21) $$

and

$$ \hat{\theta} = \tan \left( \frac{\int E(f, \theta) \sin(\theta) df d\theta}{\int E(f, \theta) \cos(\theta) df d\theta} \right) \qquad (A22) $$

The modulation due to the wind is computed as:

$$ M_w = \frac{\left( 1 + p_5 \max\left(1, \frac{f}{f_0}\right) \right)}{\left(1 + p_5\right)} \qquad (A23) $$

with

$$ f_0 = p_6 \frac{1}{u_*} \frac{g}{2\pi} \qquad (A24) $$

where $p_5 = 0.9$ is a constant related to the DIA algorithm and $p_6 = 3/28$ is yet another constant. Finally, the conversion from $\Lambda(f)$ to $\Lambda(c)$ is done using the relation $\Lambda(c)dc = \Lambda(f)df$ and the linear dispersion relation (see Romero's [2019]; Equations 17–23 for details).

## Data Availability Statement

All data used in this publication has been previously published by Banner et al. (2000), Sutherland and Melville (2013), and Schwendeman et al. (2014).

## References

Alsina, J. M., & Baldock, T. E. (2007). Improved representation of breaking wave energy dissipation in parametric wave transformation models. *Coastal Engineering*, 54, 765–769. https://doi.org/10.1016/j.coastaleng.2007.05.005

Ardag, D., & Resio, D. T. (2020). A new approach for modeling dissipation due to breaking in wind wave spectra. *Journal of Physical Oceanography*, 50(2), 439–454. https://doi.org/10.1175/jpo-d-19-0160.1

Ardhuin, F., Rogers, E., Babanin, A. V., Filipot, J.-F., Magne, R., Roland, A., et al. (2010). Semiempirical dissipation source functions for ocean waves. Part I: Definition, calibration, and validation. *Journal of Physical Oceanography*, 40(9), 1917–1941. https://doi.org/10.1175/2010JPO4324.1

Banner, M. L., Babanin, A. V., & Young, I. R. (2000). Breaking probability for dominant waves on the sea surface. *Journal of Physical Oceanography*, 30(12), 3145–3160. https://doi.org/10.1175/1520-0485(2000)030(3145:BPFDWO)2.0.CO;2

Banner, M. L., Gemmrich, J. R., & Farmer, D. M. (2002). Multiscale measurements of ocean wave breaking probability. *The Journal of Physical Oceanography*, 32(12), 3364–3375. https://doi.org/10.1175/1520-0485(2002)032(3364:MMOOWB)2.0.CO;210.1175/1520-0485(2002)032<3364:mmoowb>2.0.CO;2

Banner, M. L., & Morison, R. P. (2010). Refined source terms in wind wave models with explicit wave breaking prediction. Part I: Model framework and validation against field data. *Ocean Modelling*, 33(1–2), 177–189. http://doi.org/10.1016/j.ocemod.2010.01.002

Banner, M. L., & Tian, X. (1998). On the determination of the onset of breaking for modulating surface gravity water waves. *Journal of Fluid Mechanics*, 367, 107–137. https://doi.org/10.1017/s0022112098001517

Banner, M. L., Zappa, C. J., & Gemmrich, J. R. (2014). A note on the Phillips spectral framework for ocean whitecaps. *Journal of Physical Oceanography*, 44(7), 1727–1734. https://doi.org/10.1175/JPO-D-13-0126.1

Barthelemy, X., Banner, M. L., Peirson, W. L., Fedele, F., Allis, M., & Dias, F. (2018). On a unified breaking onset threshold for gravity waves in deep and intermediate depth water. *Journal of Fluid Mechanics*, 841, 463–488. https://doi.org/10.1017/jfm.2018.93

Battjes, J. A., & Janssen, J. (1978). Energy loss and set-up due to breaking of random waves. *Coastal Engineering*, 32(1), 569–587.

Chawla, A., & Kirby, J. T. (2002). Monochromatic and random wave breaking at blocking points. *Journal of Geophysical Research: Oceans*, 107(C7), 4–1. https://doi.org/10.1029/2001jc001042

Cokelet, E. D. (1977). Steep gravity waves in water of arbitrary uniform depth. *Philosophical Transactions of the Royal Society A: Mathematical. Physical and Engineering Sciences*, 286(1335), 183–230. https://doi.org/10.1098/rsta.1977.0113

Dean, R. G. (1965). Stream function representation of nonlinear ocean waves. *Journal of Geophysical Research*, 70(18), 4561–4572. https://doi.org/10.1029/jz070i018p04561

Derakhti, M., Kirby, J. T., Banner, M. L., Grilli, S. T., & Thomson, J. (2020). A unified breaking onset criterion for surface gravity water waves in arbitrary depth. *Journal of Geophysical Research: Oceans*, 125(2013), 1–28. https://doi.org/10.1029/2019jc015886

Donelan, M., Longuet-Higgins, M. S., & Turner, J. S. (1972). Periodicity in whitecaps. *Nature*, 239(5373), 449–451. https://doi.org/10.1038/239449a0

## Acknowledgments

This study benefited from France Energies Marines and State financing managed by the National Research Agency under the Investments for the Future program bearing the reference numbers ANR-10-IED-0006-14, ANR-10-IEED-0006-26 and ANR-10-IEED-0006-20 for the projects DiME, CARAVELE and DIMPACT. The authors thank Peter Sutherland and Jim Thompson for kindly sharing their data.
Duncan, J. H. (1981). An experimental investigation of breaking waves produced by a towed hydrofoil. *Proceedings of the Royal Society A: Mathematical. Physical and Engineering Sciences*, **377**(1770), 331–348. https://doi.org/10.1098/rspa.1981.0127

Eldeberky, Y., & Battjes, J. A. (1996). Spectral modeling of wave breaking: Application to Boussinesq equations. *Journal of Geophysical Research*, **101**, 1253–1264. https://doi.org/10.1029/95jc03219

Fady, J., Zappa, C. J., Banner, M. L., & Morison, R. P. (2013). Wave breaking in developing and mature seas. *Journal of Geophysical Research: Oceans*, **118**(9), 4542–4552. https://doi.org/10.1002/jgrc.20334

Filipot, J.-F., & Ardhuin, F. (2012). A unified spectral parameterization for wave breaking: From the deep ocean to the surf zone. *Journal of Geophysical Research*, **117**(C11), C00J08. https://doi.org/10.1029/2011JC007784

Filipot, J.-F., Ardhuin, F., & Babanin, A. V. (2010). A unified deep-to-shallow water wave-breaking probability parameterization. *Journal of Geophysical Research*, **115**(4), 1–15. https://doi.org/10.1029/2009JC005448

Filipot, J.-F., Guimaraes, P., Leckler, F., Hortsmann, J., Carrasco, R., Leroy, E., et al. (2019). La Jument lighthouse: A real-scale laboratory for the study of giant waves and their loading on marine structures. *Proceedings of the Royal Society A: Mathematical. Physical and Engineering Sciences*, **377**(2155), 20190008. https://doi.org/10.1098/rsta.2019.0008

Guimarães, P. V. (2018). *Sea surface and energy dissipation (Unpublished doctoral dissertation)*. Université de Bretagne Loire.

Holland, K. T., Holman, R. A., Lippmann, T. C., Stanley, J., Plant, N., & Plant, N. (1997). Practical use of video imagery in nearshore oceanographic field studies. *IEEE Journal of Oceanic Engineering*, **22**(1), 81–92. https://doi.org/10.1109/48.557542

Holthuijsen, L. H., & Herbers, T. H. C. (1986). Statistics of breaking waves observed as whitecaps in the open sea. *Journal of Physical Oceanography*, **16**(2), 290–297. https://doi.org/10.1175/1520-0485(1986)016<0290:sobwoa>2.0.co;2

Hwang, P. A., Wang, D. W., Walsh, E. J., Krabill, W. B., & Swift, R. N. (2000). Airborne measurements of the wavenumber spectra of ocean surface waves. Part II: Directional distribution. *Journal of Physical Oceanography*, **30**(11), 2768–2787. https://doi.org/10.1175/1520-0485(2001)031(2768:amotws)2.0.co;2

Janssen, T. T., & Battjes, J. A. (2007). A note on wave energy dissipation over steep beaches. *Coastal Engineering*, **54**(9), 711–716. https://doi.org/10.1016/j.coastaleng.2007.05.006

Kjeldsen, S. P., Vinje, T. P., Myrhaug, D. P., & Brdvig, P. P. (1980). Kinematics of deep water breaking waves. In *Paper presented at the off-shore technology conference*. https://doi.org/10.4043/3714-ms

Kleiss, J. M., & Melville, W. K. (2011). The analysis of sea surface imagery for whitecap kinematics. *Journal of Atmospheric and Oceanic Technology*, **28**(2), 219–243. https://doi.org/10.1175/2010JTECHO744.1

Kudryavtsev, V., Chapron, B., & Makin, V. (2014). Impact of wind waves on the air-sea fluxes: A coupled model. *Journal of Geophysical Research: Oceans*, **119**(2), 1217–1236. https://doi.org/10.1002/2013jc009412

Le Mehaute, B. (1962). On non-saturated breakers and the wave run-up. *Coastal Engineering Proceedings*, **1**(8), 6. https://doi.org/10.9753/icce.v8.6

Longuet-Higgins, M. S. (1957). The statistical analysis of a random, moving surface. *Philosophical Transactions of the Royal Society A: Mathematical, Physical & Engineering Sciences*, **249**(966), 321–387. https://doi.org/10.1098/rsta.1957.0002

Melville, W. K., & Matusov, P. (2002). Distribution of breaking waves at the ocean surface. *Nature*, **417**(6884), 58–63. https://doi.org/10.1038/417058a

Miche, A. (1944). Mouvements ondulatoires de la mer en profondeur croissante ou décroissante. *Annales des Ponts et Chaussées*, **114**, 42–78.

Perlin, M., Choi, W., & Tian, Z. (2013). Breaking waves in deep and intermediate waters. *Annual Review of Fluid Mechanics*, **45**(1), 115–145. https://doi.org/10.1146/annurev-fluid-011212-140721

Phillips, O. M. (1985). Spectral and statistical properties of the equilibrium range in wind-generated gravity waves. *Journal of Fluid Mechanics*, **156**, 505–531. https://doi.org/10.1017/S0022112085002221

Reul, N., & Chapron, B. (2003). A model of sea-foam thickness distribution for passive microwave remote sensing applications. *Journal of Geophysical Research*, **108**(10), 3321. https://doi.org/10.1029/2003JC001887

Rice, S. O. (1944). Mathematical analysis of random noise. *The Bell System Technical Journal*, **23**(3), 282–332. https://doi.org/10.1002/j.1538-7305.1944.tb00874.x

Romero, L. (2019). Distribution of surface wave breaking fronts. *Geophysical Research Letters*, **46**(17–18), 10463–10474. https://doi.org/10.1029/2019GL083408

Saket, A., Peirson, W. L., Banner, M. L., Barthelemy, X., & Allis, M. J. (2017). 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**, 642. https://doi.org/10.1017/jfm.2016.776

Schwendeman, M., Thomson, J., & Gemmrich, J. R. (2014). Wave breaking dissipation in a Young Wind Sea. *Journal of Physical Oceanography*, **44**(1), 104–127. https://doi.org/10.1175/JPO-D-12-0237.1

Sutherland, P., & Melville, W. K. (2013). Field measurements and scaling of ocean surface wave-breaking statistics. *Geophysical Research Letters*, **40**(12), 3074–3079. https://doi.org/10.1002/grl.50584

Thomson, J. (2012). Wave breaking dissipation observed with "SWIFT" drifters. *Journal of Atmospheric and Oceanic Technology*, **29**(12), 1866–1882. https://doi.org/10.1175/JTECH-D-12-00018.1

Thomson, J., & Jessup, A. T. (2009). A Fourier-based method for the distribution of breaking crests from video observations. *Journal of Atmospheric and Oceanic Technology*, **26**(8), 1663–1671. https://doi.org/10.1175/2009JTECHO622.1

Thornton, E. B., & Guza, R. T. (1983). Transformation of wave height distribution. *Journal of Geophysical Research*, **88**(C10), 5925–5938. https://doi.org/10.1029/jc088ic10p05925

Varing, A., Filipot, J.-f., Grilli, S., Duarte, R., Roeber, V., & Yates, M. (2020). A new kinematic breaking onset criterion for spilling and plunging breaking waves in shallow water. *Coastal Engineering*, **1–24**.

Zieger, S., Babanin, A. V., Erick Rogers, W., & Young, I. R. (2015). Observation-based source terms in the third-generation wave model WAVEWATCH. *Ocean Modelling*, **96**(1), 2–25. https://doi.org/10.1016/j.oocemod.2015.07.014
3.B Paper from Stringari & al. (2021a)
**scientific reports**

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

Check for updates

OPEN

# Deep neural networks for active wave breaking classification

Caio Eadi Stringari<fnref n="1" />,², Pedro Veras Guimarães<fnref n="1" />,²,³, Jean-François Filipot<fnref n="1" />,³, Fabien Leckler<fnref n="1" /> & Rui Duarte<fnref n="1" />

Wave breaking is an important process for energy dissipation in the open ocean and coastal seas. It drives beach morphodynamics, controls air-sea interactions, determines when ship and offshore structure operations can occur safely, and influences on the retrieval of ocean properties from satellites. Still, wave breaking lacks a proper physical understanding mainly due to scarce observational field data. Consequently, new methods and data are required to improve our current understanding of this process. In this paper we present a novel machine learning method to detect active wave breaking, that is, waves that are actively generating visible bubble entrainment in video imagery data. The present method is based on classical machine learning and deep learning techniques and is made freely available to the community alongside this publication. The results indicate that our best performing model had a balanced classification accuracy score of ≈90% when classifying active wave breaking in the test dataset. An example of a direct application of the method includes a statistical description of geometrical and kinematic properties of breaking waves. We expect that the present method and the associated dataset will be crucial for future research related to wave breaking in several areas of research, which include but are not limited to: improving operational forecast models, developing risk assessment and coastal management tools, and refining the retrieval of remotely sensed ocean properties.

Wave breaking is one of the most challenging water wave phenomena to investigate. Despite nearly two centuries of research, the current practical description of the process is still mostly empirical<fnref n="1" />⁻⁴. Precise wave breaking modelling requires explicit knowledge of the wave phase speed and the fluid velocity distribution on the wave crest, which can only be done by numerically solving the Navier-Stokes equations in a framework that currently is too computationally expensive for practical applications⁵. Consequently, current large-scale state-of-the-art wave models rely on statistical approaches (mostly the spectral approach) to represent the waves⁶⁻⁸. This type of models parameterizes wave breaking as a function of known parameters such as the wind speed⁶, the local wave height to water depth ratio⁷ or semi-empirical breaking wave height probability distributions²,⁹. Due to a lack of observed data, the constants involved in these models have been derived from limited datasets that may not adequately represent the natural environment. For example, a recent study has shown that popular surf zone parametric wave breaking models incorrectly represented the fraction of broken waves in their formulations with errors > 50%, despite these models being able to adequately represent surf zone energy dissipation possibly due to parameter tuning¹⁰. This paper aims to start addressing the data unavailability issue by providing a reliable and reproducible method to detect and track waves that are actively breaking in video imagery data.

Wave breaking directly affects several environmental phenomena. For example, wave breaking is considered as the main driver for air-sea exchanges¹¹, being often described as a function of Phillips'¹² Λ(c)dc parameter. This parameter is defined as the average total length per unit surface area of breaking fronts that have velocities in the range c to c + dc and its moments are assumed to correspond to different phenomena such as whitecap coverage (second moment), rate of air entrainment (third moment), and momentum flux (fourth moment). As previously identified¹³, different interpretations of how data is processed to obtain Λ(c)dc resulted in differences of 175% in the first moment and 300% in the fifth moment of Λ(c)¹³,¹⁴. This issue does not seem to limit the applicability of Phillips' framework. A recent study, for example, used a wave breaking parameterization fully based on Λ(c)dc to model whitecap coverage off the coast of California¹⁵. To the authors' knowledge, however, no study has assessed the impact that errors in wave breaking detection has on measured Λ(c) distributions and how such errors would impact on the conclusions that have been drawn from these models (for example, the assumption that Λ(c) is a simple function of wind speed¹¹). More importantly, in the case which new observations

<footnote n="1">¹France Energies Marines, 29280 Plouzané, France. ²PPGOceano, Federal University of Santa Catarina, Florianópolis 88040-900, Brazil. ³These authors contributed equally: Pedro V. Guimarães and Jean-François Filipot. ⁴email: Caio.Stringari@france-energies-marines.com</footnote>
diverge from Phillips' original theory (as we will later show in this paper), new and improved wave breaking models will certainly follow.

Further, breaking waves have a direct impact on the remote sensing of the ocean surface from satellites, aircraft, or platform mounted instruments. Wave breaking has been observed to significantly modulate measurements of backscatter off radars with large increases in backscatter being directly correlated with breaking waves¹⁶. Given that part of the data that will be used here (see Methods for details) coincides with the location of previous wave backscatter studies¹⁶, data derived from our method could lead to explicit formulations correlating wave breaking and radar backscatter. Such formulations could, in the future, lead to more reliable satellite-derived scatterometer data, which currently ignores¹⁷ or relies on empirical wave breaking models¹⁸. The foam generated by breaking waves at the sea surface also influences on the measurement of microwave brightness temperature signatures for moderate to high wind speeds¹⁹,²⁰. To date, no study has, however, systematically quantified such temperature signatures for large field datasets and most of our knowledge is empirical. Such limitation could be addressed by extending wave breaking detection and tracking methods to infrared wavelengths, which is directly correlated to water temperature²¹. Advancements in this regard could lead to improvements in global climate models given that the inclusion of wave-generated heat transfer has shown to reduce cold biases in such models for non-breaking waves²². Breaking waves are expected to transfer more heat than non-breaking waves²³ and, when considered, could improve climate models even further.

With the recent explosion in the usage of machine learning, particularly deep neural networks, it was only a matter of time before these techniques were adapted to wave research. Recent studies have shown that deep neural networks can accurately be used to classify different types of surf zone breakers²⁴, obtain wave heights from video data²⁵, track waves in the surf zone²⁶ and the shoreline position²⁷ and, when applied to pressure transducer data, distinguish between broken and unbroken waves¹⁰. In this paper, we describe a robust and extensible method to identify and track breaking waves in video imagery data using a combination of classic machine learning algorithms and deep learning. More importantly, we make the present dataset and code base fully available for future researchers who can either use it directly, re-train the models adding more training data, or expand the present method to other cases (for example, object segmentation). The development and standardization of a general framework to detect wave breaking in video imagery data should help to provide the wave breaking statistics database that is currently needed to assess, or further develop, our current understanding of several processes related to breaking waves. Ultimately, the efforts initiated here aim to lead to improvements in several areas of research such as wave modelling and forecasting, wave energy tracking and harvesting²⁸, ocean-atmosphere interaction¹⁵, remote sensing of the ocean¹⁶, and safety at sea and coasts²⁹.

## Method

**Model definition.** In this study, we have developed a novel method to detect active wave breaking in video imagery data. Similarly to previous studies, we exploit the fact that breaking wave crests generate characteristic white foam patches from bubble entrainment³⁰,³¹. Differently from the vast majority of previous methods, here we make a clear distinction between active wave breaking (that is, visible foam being actively generated by wave breaking) from passive foam and to all other non-wave breaking instances. To the authors' knowledge, only the methods of Mironov and Dulov³⁰ and Kleiss and Melville³² have previously attempted to distinguish between active wave breaking and passive foam. Both their methods start similarly to ours (see next paragraph) by identifying bright pixels in a given series of images via thresholding. Next, interconnected pixels are connected in space and time using a nearest neighbor search (Mironov and Dulov³⁰) or by extracting connected contours (Kleiss and Melville³²). In both previous methods, the separation between active and passive wave breaking is then done by empirically defining foam expansion direction and velocity limits. The present method aims to remove empirical steps in active wave breaking detection by using data-driven methods and deep neural networks instead of experiment-specific parameters.

Our method starts by naively identifying wave breaking candidates in a given input image. This was done by identifying bright (that is, white) pixels in the image using the local thresholding algorithm available from the OpenCV³³ library with its default parameters. Interconnected regions of bright pixels in the image were then clustered using DBSCAN³⁴ and a minimum enclosing ellipse³⁵ was fitted to each identified cluster. The DBSCAN step requires the user to define a minimum number of pixels to be clustered (default value is 10 for all data used here) and a maximum distance allowed between pixels (default value of 5 for all data used here). This step resulted in large quantities of bright pixel patches being detected. At this stage, it was not possible to determine whether a given identified cluster of pixels represented active wave breaking, passive foam, or any other instance of bright pixels (for example, birds, boats, coastal structures, sun glints or white clouds). Generically, this step can be thought of a feature extraction step and can be replaced by any other equivalent approach (for example, image feature extractors). To avoid exponential memory consumption growth generated by DBSCAN, the input image may be subdivided into blocks of regular size. This step was easily parallelized with performance increasing linearly with number of computing threads.

The second step of the method consisted of training a deep convolutional neural network to distinguish between active wave breaking and passive foam (and all other non-desired occurrences). This step can be understood as a binary classification step. Figure 1 shows a schematic representation of the present method. From the original images, randomly selected subsets of size 256 × 256 pixels centered on the fitted ellipses were extracted and manually labelled as either active wave breaking (label 1, or positive) or otherwise (label 0, or negative). The present training dataset was generated using raw video imagery data from Guimarães et al.³⁶ and consisted of 19,000 training and 1300 testing images (see Table 1). The training dataset is further split into 80% training and 20% validation data³⁷. The validation dataset is reserved for future hyper-parameter fine-tuning. Note that the present dataset has a class imbalance of ≈ 90% towards the negative label, that is, for each active wave breaking