# Computing local crest phase speed in multi-directional sea states through the Spatial Hilbert Transform

Ying Wang, Guillaume Ducrozet

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

## To cite this version:

Ying Wang, Guillaume Ducrozet. Computing local crest phase speed in multi-directional sea states through the Spatial Hilbert Transform. Ocean Engineering, 2025, 331, pp.121357. 10.1016/j.oceaneng.2025.121357. hal-05053652

HAL Id: hal-05053652

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

Submitted on 2 May 2025

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

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

<figure><img src="image_2.png" /></figure>
<figure><img src="image_3.png" /></figure>

Research paper

# Computing local crest phase speed in multi-directional sea states through the Spatial Hilbert Transform

Ying Wang*, Guillaume Ducrozet ID

Nantes Université, École Centrale Nantes, CNRS, LHEEA, UMR 6598, F-44000, Nantes, France

## ARTICLE INFO

### Keywords:

- Wave breaking

- Wave phase speed

- Hilbert transform

- Short-crested sea states

- HOS

- SHTM

## ABSTRACT

The detection of wave breaking plays a fundamental role in the modelling of this complex physical phenomenon within nonlinear potential flow solvers. Recent studies have demonstrated that utilizing the ratio of fluid velocity to local wave crest phase speed enables accurate identification of individual breaking waves in a wave field. However, evaluating the local wave crest phase speed poses great challenges, especially for multi-directional sea states. Based on the multi-dimensional Hilbert transform in space, we correlated the local phase function and amplitude components to a mixed local amplitude (envelope) and phase function. Then, associated with these slow-varying assumptions of the local quantities, two methodologies, the linear spatial Hilbert transform method and the spatial Hilbert transform method, are proposed to compute the local wave crest phase speed in multi-directional sea states. These approaches are first compared to a classical crest-tracking method for short-crested waves and then investigated in detail by accounting for the existence of nonlinearity, directionality, and spectrum bandwidth. The present work shows that the proposed spatial Hilbert transform method allows for an efficient and accurate evaluation of the local wave crest phase speed in multi-directional sea states, particularly for strongly nonlinear configurations.

## 1. Introduction

In the applications of ocean engineering studies, the modelling of environmental conditions over large spatial and temporal scales is required. This is usually achieved using nonlinear potential flow solvers that are limited to non-breaking waves. However, in real sea states, breaking waves are not only prevalent but also play a significant role, which gives rise to the requirement to extend the application range of nonlinear potential flow solvers to be able to include breaking waves. This involves two parts: firstly, the detection of individual breaking events, and secondly the energy dissipation due to wave breaking. For the detection, different types of breaking onsets have been developed, and a comprehensive review can be found in Perlin et al. (2013), in which the pros and cons of each type of breaking criterion are systematically discussed. Thanks to Barthelemy et al. (2018) who proposed a breaking onset by using the quotient of the local energy flux velocity and the local crest phase speed, which eventually reduces to kinematic breaking onset on the free surface, i.e. $B = ||u||/||c||$ where $||·||$ represents the norm of the (corresponding) vector. They argued that the wave evolves to breaking when the breaking parameter $B$ at wave crest exceeds a specific value 0.85, which has been extensively validated by Saket et al. (2017), Derakhti et al. (2020) and Seiffert et al. (2017). It is

worth noting that the 'local' character refers to quantities locally defined both in space and time.

For such a kinematic breaking onset, evaluating the local wave crest phase speed plays a crucial role. Its accurate and efficient computation is the key point of Barthelemy's breaking onset. However, it poses a great challenge, especially in multi-directional sea states. Intuitively, one can use the so-called Crest Tracking Method (CTM) to compute it, which directly tracks the wave crest location at every time instant; the local wave crest phase speed is then computed by taking the temporal derivative of the wave crest locations collected. This can be considered the classical method, implemented for instance by Derakhti et al. (2020), Mohanlal et al. (2023), Barthelemy et al. (2018), Touboul and Banner (2021), Derakhti et al. (2018) and Seiffert et al. (2017). However, while the CTM is applicable and straightforward for uni-directional waves, it becomes tricky for complex multi-directional wave fields. For stability, some numerical treatments, such as splines to spatially smooth the free surface and temporal filters to smooth the crest phase speed itself, are required but can be computationally expensive in multi-directional sea states.

Unlike conventional CTM, the methods based on the use of the Hilbert Transform (HT) to compute the local crest phase speed do not need to track the wave crest location. Instead, they retain the local

<footnote>* Corresponding author.
E-mail address: ying.wang@ec-nantes.fr (Y. Wang).</footnote>
character of the wave celerity, which is more flexible and easy to set up to complex sea states. The basic idea of the Hilbert Transform Methods (HTMs) is to take the Temporal Hilbert Transform (THT) or Spatial Hilbert Transform (SHT) of the surface elevation to construct the local phase function, and then use different strategies to find the local crest phase speed. Originally, Melville (1983) proposed to use the THT to post-process of the experimental data in order to compute the local phase function in uni-directional sea states, and then to take the spatio-temporal derivative of this local phase function to find the local angular frequency $\omega(x, t)$ and wavenumber $k(x, t)$ to evaluate the local wave phase speed, $c(x, t) = \omega(x, t)/k(x, t)$. This will be referred to as the Temporal Hilbert Transform Method (THTM), which was widely used to assess the local wave phase speed for post-processing experimental data (Hwang et al., 1989; Wu and Nepf, 2002). Interestingly, Stansell and MacFarlane (2002) proposed to utilize the partial HT in time to find the local angular frequency and the partial HT in space to calculate the local wavenumber, and then to compute the corresponding local wave phase speed. A different approach was adopted by Kurnia and Van Groesen (2014), Seiffert et al. (2017) and Craciuinescu and Christou (2020), which involved taking the SHT of free surface elevation along the direction of propagation to determine the local wavenumber, associated with the linear dispersion relation to obtain the local wave crest phase speed in unidirectional sea states.

All the previous studies regarding the use of SHT are limited to unidirectional waves, while those using THT are confined to post-processing. In the current work, based on the multi-dimensional SHT, we correlated the local phase function and amplitude components to a mixed local amplitude and phase function. Associated with slow-varying assumptions of the local quantities, i.e. local wavenumber $k(x, t)$ and local angular frequency $\omega(x, t)$, Whitham's modulation theory (Whitham, 1965) is derived through the proposed SHT expression in multi-directional sea states. Then, combined with the $k$ condition (a condition that involves taking the spatial derivative of the local phase function to compute the local wavenumber) of Whitham's modulation theory and the linear dispersion relation, the Linear Spatial Hilbert Transform Method (LSHTM) is generalized to multi-directional sea states. By fully implementing the $k-\omega$ conditions to the mixed local phase function, the Spatial Hilbert Transform Method (SHTM) is proposed to compute the local wave phase speed in multi-directional sea states.

For validation, linear and nonlinear regular waves with and without propagating angles, for which an analytical solution of the crest phase speed is available, are initially investigated. Then, the role of the nonlinearity, spectrum bandwidth, and directionality of the wave in the aforementioned methodologies for local crest phase speed evaluation are studied via four frequency focusing waves in unidirectional sea states, and four frequency and directional focusing waves in multi-directional sea states. In addition, through uni- and multi-directional irregular waves, the correctness of the proposed mixed local amplitude and phase function is validated by using the SHT-constructed surface elevation to reproduce the original surface elevation.

The rest of the paper is organized as follows. In Section 2, Whitham's modulation theory (Whitham, 1965) is first presented, followed by the fundamental mathematical definition of the multi-dimensional SHT. Then, we correlate the locally defined amplitude and phase function components to a mixed local amplitude and phase function and discuss the role played by the slow-varying assumption of the local quantities in their evaluation, based on which the LSHTM and SHTM were deduced. The CTM for directional sea states is presented as a reference. In Section 3, the key point of the numerical models we used is briefly introduced, and the forthcoming application of the current work is briefly discussed. In Section 4, linear and nonlinear regular waves with and without propagating angles, $\beta = 0^{\circ}$ and $\beta = 45^{\circ}$, are studied, followed by frequency focusing waves, and lastly by combined frequency and directional focusing waves. The paper closes with a brief conclusion and some perspectives for future work.

## 2. Methodologies

In the modelling of wave breaking in a strongly nonlinear potential-flow model, the free surface elevation $z = \eta(x, t)$ is known at each time instant $t$. $x = (x, y)$ represents the horizontal coordinates, while $z$ stands for the vertical coordinate, with $z = 0$ located at the mean undisturbed free surface. Different methodologies have been developed to evaluate the local wave phase speed from this spatio-temporal field.

### 2.1. Multi-directional Hilbert transform methods

Based on the slow-varying assumption of the locally defined wavenumber $k(x, t)$, angular frequency $\omega(x, t)$, and amplitude $A(x, t)$, Whitham's modulation theory, developed by Whitham (1965) and Lighthill (1978), is first presented. For our purpose of instantaneously computing the local wave crest phase speed at every instant in which the free surface elevation is available, the temporal variable will be treated as a constant for now. Then, the multi-dimensional SHT and the so-called complex signals are presented. Based on this, we developed a novel mixed local amplitude and phase function, which allows us to uniquely define the local phase function and amplitude of the free surface elevation, and bridge them with Whitham's theory. The role of the wave spectrum bandwidth in the evaluation of local $(k, \omega)$ is discussed, and two multi-directional HTMs are proposed to evaluate the local crest phase speed.

#### 2.1.1. Whitham's modulation theory

Whitham (1962) first proposed the concepts of local wavenumber, and angular frequency, which are slowly varying functions in time $t$ and space $x$, and defined via a phase function $\varphi(x, t)$ in the uni-directional wave. Whitham (1965) systematically introduced the local $(k, \omega, A)$ as slowly varying functions of time and space, and mathematically developed the so-called modulation theory, which lays the cornerstone of this theory. In the same paper, he explicitly explained what the slow variation of the local $(k, \omega, A)$ is and how to determine it. It indicates that two scales are involved: the smaller scale is determined by the wavelength and period of the oscillations, during which the local quantities remain roughly constant, and the larger scale is determined by the $x$ and $t$ themselves, over which those local quantities change appreciably. Besides, two treatises, written by Whitham (1974) and Lighthill (1978), provide comprehensive details and insightful discussion of this modulation theory. This modulation theory, generalised to multi-directional sea states, allows us to define a local phase function, $\varphi(x, t)$. The value of this quantity at every crest is an even multiple of $\pi$, and at troughs an odd multiple of $\pi$. The water surface displacement $\eta(x, t)$ is therefore given by

$$ \eta(x, t) = \mathfrak{R}\{A(x, t) \exp[j\varphi(x, t)]\} \qquad (2.1) $$

where $A(x, t)$ is a positive, slowly varying amplitude, $\mathfrak{R}[\cdot]$ is the function to take the real part, and $j$ represents the unit imaginary number, $j = \sqrt{-1}$. Correspondingly, based on the slow-varying assumption of the local $(k, \omega)$, the phase function should satisfy

$$ \nabla\varphi(x, t) = \mathbf{k}(x, t) \qquad (2.2) $$

$$ \frac{\partial\varphi(x, t)}{\partial t} = -\omega(x, t) \qquad (2.3) $$

which will be named as $k-\omega$ conditions, and $\nabla$ represents the horizontal gradient, $\nabla = (\partial/\partial x, \partial/\partial y)$. The wave propagation direction is therefore determined by

$$ \mathbf{n}(x, t) = \mathbf{k}(x, t)/k(x, t) \qquad (2.4) $$

where $k(x, t)$ is the local wavenumber, obtained by taking the modulus of the local wavenumber vector

$$ k(x, t) = \|\mathbf{k}(x, t)\| = \sqrt{(k_x(x, t))^2 + (k_y(x, t))^2} $$
Additionally, from Eqs. (2.2) and (2.3), a relationship between wavenumber and angular frequency exists

$$ \frac{\partial k(\mathbf{x}, t)}{\partial t} + \nabla\omega(\mathbf{x}, t) = \mathbf{0} \qquad (2.5) $$

which can be interpreted as an equation of continuity for the local phase. Whitham's modulation equation is widely used to describe waves in various fields such as water waves (Chu and Mei, 1970; Esler and Pearce, 2011), nonlinear optics (El et al., 2007), meteorology (Porter and Smyth, 2002), oeanography (Smyth and Holloway, 1988), etc. More aspects related to the application of Whitham's modulation theory can be found in a comprehensive review (Minzoni and Smyth, 2016).

### 2.1.2. Multi-dimensional spatial Hilbert transform analysis

The Hilbert Transform (HT) is not new, and dates from David Hilbert's work (Hilbert, 1905) in which he introduced the HT to solve a special case of the Riemann-Hilbert problem. Several progresses have been made since Hilbert's original work: Gabor (1946) proposed the analytical signal, Poularikas and Grigoryan (2018) constructed complex signals for multi-dimensional signals, Huang et al. (1998) introduced the so-called Hilbert-Huang transform, etc. Nowadays, the HT is widely used in different fields, such as biomedicine (Chavez et al., 2006; Raut et al., 2021; Schäfer et al., 1998), signal processing (Potamianos and Maragos, 1994), fluid mechanics (Huang et al., 1999; Krechetnikov, 2024), acoustics (Feldman, 2011), etc. It is mainly applied to temporal signals, which are one-dimensional. What is proposed in the present work is to apply this methodology to 2D signals, which in our case are spatial coordinates of the free surface elevation. From a signal processing point of view, in short-crested seas, the surface elevation $\eta(\mathbf{x}, t)$ is a three-dimensional signal, i.e. two spatial dimensions and one temporal dimension. As aforementioned, we have surface elevation information in space at each instant from which we want to extract features, so for our application we fix the temporal variable, and the three-dimensional signal is therefore reduced to two spatial dimensions. The partial HT in $x$ and $y$, and 2D mixed Spatial Hilbert Transforms (SHT) of the wave elevation $\eta(\mathbf{x}, t)$ are defined by

$$ \mathbb{H}_x[\eta(\mathbf{x}, t)] = \frac{1}{\pi} \text{p.v.} \int_{-\infty}^{\infty} \frac{\eta(x', y, t)}{x - x'} dx' \qquad (2.6) $$

$$ \mathbb{H}_y[\eta(\mathbf{x}, t)] = \frac{1}{\pi} \text{p.v.} \int_{-\infty}^{\infty} \frac{\eta(x, y', t)}{y - y'} dy' \qquad (2.7) $$

$$ \mathbb{H}_{xy}[\eta(\mathbf{x}, t)] = \frac{1}{\pi^2} \text{p.v.} \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \frac{\eta(x', y', t)}{(x - x')(y - y')} dx' dy' \qquad (2.8) $$

where p.v. is the Cauchy principal value of the integral, $\mathbb{H}_{xy}[\eta(\mathbf{x}, t)]$ denotes the mixed 2D SHT, and $\mathbb{H}_x[\eta(\mathbf{x}, t)]$ and $\mathbb{H}_y[\eta(\mathbf{x}, t)]$ are the partial SHTs with respect to $x$ and $y$, respectively. For the two-dimensional signals, two complex signals are constructed:

$$ \Psi_1(\mathbf{x}, t) = \eta(\mathbf{x}, t) - \mathbb{H}_{xy}[\eta(\mathbf{x}, t)] + j(\mathbb{H}_x[\eta(\mathbf{x}, t)] + \mathbb{H}_y[\eta(\mathbf{x}, t)]) \qquad (2.9) $$

$$ \Psi_2(\mathbf{x}, t) = \eta(\mathbf{x}, t) + \mathbb{H}_{xy}[\eta(\mathbf{x}, t)] + j(\mathbb{H}_x[\eta(\mathbf{x}, t)] - \mathbb{H}_y[\eta(\mathbf{x}, t)]) \qquad (2.10) $$

In polar form, Eqs. (2.9) to (2.10) can be rewritten as

$$ \Psi_1(\mathbf{x}, t) = A_1(\mathbf{x}, t) \exp\{j\phi_1(\mathbf{x}, t)\} \qquad (2.11) $$

$$ \Psi_2(\mathbf{x}, t) = A_2(\mathbf{x}, t) \exp\{j\phi_2(\mathbf{x}, t)\} \qquad (2.12) $$

where $A_1(\mathbf{x}, t)$ and $A_2(\mathbf{x}, t)$ are the local amplitude components defined by

$$ A_1(\mathbf{x}, t) = \sqrt{(\eta(\mathbf{x}, t) - \mathbb{H}_{xy}[\eta(\mathbf{x}, t)])^2 + (\mathbb{H}_x[\eta(\mathbf{x}, t)] + \mathbb{H}_y[\eta(\mathbf{x}, t)])^2} \qquad (2.13) $$

$$ A_2(\mathbf{x}, t) = \sqrt{(\eta(\mathbf{x}, t) + \mathbb{H}_{xy}[\eta(\mathbf{x}, t)])^2 + (\mathbb{H}_x[\eta(\mathbf{x}, t)] - \mathbb{H}_y[\eta(\mathbf{x}, t)])^2} \qquad (2.14) $$

$\phi_1(\mathbf{x}, t)$ and $\phi_2(\mathbf{x}, t)$ are the local phase function components defined by

$$ \phi_1(\mathbf{x}, t) = \arctan \left( \frac{\mathbb{H}_x[\eta(\mathbf{x}, t)] + \mathbb{H}_y[\eta(\mathbf{x}, t)]}{\eta(\mathbf{x}, t) - \mathbb{H}_{xy}[\eta(\mathbf{x}, t)]} \right) \qquad (2.15) $$

$$ \phi_2(\mathbf{x}, t) = \arctan \left( \frac{\mathbb{H}_x[\eta(\mathbf{x}, t)] - \mathbb{H}_y[\eta(\mathbf{x}, t)]}{\eta(\mathbf{x}, t) + \mathbb{H}_{xy}[\eta(\mathbf{x}, t)]} \right) \qquad (2.16) $$

From Eqs. (2.11) to (2.16), one finds

$$ \eta(\mathbf{x}, t) - \mathbb{H}_{xy}[\eta(\mathbf{x}, t)] = A_1(\mathbf{x}, t) \cos\{\phi_1(\mathbf{x}, t)\} \qquad (2.17) $$

$$ \eta(\mathbf{x}, t) + \mathbb{H}_{xy}[\eta(\mathbf{x}, t)] = A_2(\mathbf{x}, t) \cos\{\phi_2(\mathbf{x}, t)\} \qquad (2.18) $$

With Eqs. (2.17) and (2.18), it is possible to reconstruct the original free surface elevation with local amplitude and phase components

$$ \eta(\mathbf{x}, t) = \frac{1}{2} \Re [A_1(\mathbf{x}, t) \exp\{j\phi_1(\mathbf{x}, t)\} + A_2(\mathbf{x}, t) \exp\{j\phi_2(\mathbf{x}, t)\}] \qquad (2.19) $$

$$ \eta(\mathbf{x}, t) = A_{12}(\mathbf{x}, t) \cos\{\phi_{12}(\mathbf{x}, t)\} \qquad (2.20) $$

where $A_{12}(\mathbf{x}, t)$ represents a mixed local amplitude modulation, i.e. wave envelope, and $\phi_{12}(\mathbf{x}, t)$ a new mixed local phase function, which is respectively defined by

$$ A_{12}(\mathbf{x}, t) = \frac{1}{2} \sqrt{(A_1(\mathbf{x}, t))^2 + (A_2(\mathbf{x}, t))^2 + 2A_1(\mathbf{x}, t)A_2(\mathbf{x}, t) \cos(\phi_1(\mathbf{x}, t) - \phi_2(\mathbf{x}, t))} \qquad (2.21) $$

$$ \phi_{12}(\mathbf{x}, t) = \arctan \left( \frac{A_1(\mathbf{x}, t) \sin(\phi_1(\mathbf{x}, t)) + A_2(\mathbf{x}, t) \sin(\phi_2(\mathbf{x}, t))}{A_1(\mathbf{x}, t) \cos(\phi_1(\mathbf{x}, t)) + A_2(\mathbf{x}, t) \cos(\phi_2(\mathbf{x}, t))} \right) \qquad (2.22) $$

Theoretically, we have to point out that as the proposed mixed local amplitude and phase functions are constructed from the original 2D signal (the free surface elevation in the current work), and its corresponding SHT, all the information contained in the original signal is also present in this representation. Indeed, from the construction process of Eq. (2.20), one can easily find that there is no slow-varying assumption made, which underlines that Eqs. (2.20) and (2.1) do not need any assumption on the bandwidth of the signal and are consequently generic. In the next short subsection, we will discuss the role of slow variation in the evaluation of the local quantities in detail.

### 2.1.3. The role of bandwidth in the evaluation of local $k-\omega$

The slow-varying property of quantities is related to the frequency content of the signal and it will hold as soon as the signal is narrow-banded. For water waves, it is reasonable to write the mixed local phase function in the following form

$$ \phi_{12}(\mathbf{x}, t) = \mathbf{k}(\mathbf{x}, t) \cdot \mathbf{x} - \omega(\mathbf{x}, t)t \qquad (2.23) $$

Taking the spatial gradient and temporal derivative of Eq. (2.23) yields

$$ \nabla\phi_{12}(\mathbf{x}, t) = \mathbf{k}(\mathbf{x}, t) + \underbrace{[\nabla\mathbf{k}(\mathbf{x}, t)] \cdot \mathbf{x} - \nabla\omega(\mathbf{x}, t)t}_{\text{additional terms}} \qquad (2.24) $$

$$ \frac{\partial\phi_{12}(\mathbf{x}, t)}{\partial t} = -\omega(\mathbf{x}, t) + \underbrace{\frac{\partial\mathbf{k}(\mathbf{x}, t)}{\partial t} \cdot \mathbf{x} - \frac{\partial\omega(\mathbf{x}, t)}{\partial t}t}_{\text{additional terms}} \qquad (2.25) $$

where $[\nabla\mathbf{k}(\mathbf{x}, t)]$ is a second order tensor. Apparently, if the signal is broadband, i.e. the local quantities are not slow-varying with $x$ and $t$, the additional terms in Eqs. (2.24) and (2.25) can become non-negligible. However, those additional terms can be removed provided the signal is narrow-band, namely the local quantities are slow-varying in space and time.

To use the $k-\omega$ conditions to derive the local wavenumber and angular frequency as well as the conservation equation of waves as shown in Eq. (2.5), the slow variation will be made to $k-\omega$ such that the additional terms in Eqs. (2.24) and (2.25) can be removed. However, for the current work, the slow-varying assumption of the local amplitude is not necessary, so the mixed local amplitude will exactly hold for the wave envelope without any slow-varying assumption as will be illustrated in Section 4.
### 2.1.4. Crest phase speed evaluation through multi-dimensional SHT

Through the multi-dimensional SHT, we deduced a mixed local amplitude and phase function in 2.1.2, which strictly constructs the wave envelope and local phase function of the surface elevation, as discussed in 2.1.3. Regarding the objectives of the current work, i.e. obtaining the local crest phase speed at each time instant from the instantaneous free surface elevation, two strategies have been developed based on slow-varying assumptions of local $(k, \omega)$, which will be individually illustrated below. The first one makes use of the spatial information (phase) solely at the time step $t$ to obtain the local crest phase speed and consequently needs the use of the linear dispersion relation; the second one uses the information of the phase at several time steps up to $t$.

**Linear Spatial Hilbert Transform Method (LSHTM):** Implementing the $k$ condition, i.e. Eq. (2.2), to the mixed local phase function Eq. (2.22) to compute the local wavenumber vector $k(x, t)$ associated with the linear dispersion relation, the local phase speed $c(x, t)$ can be determined:

$$c_{LSHTM}(x, t) = \frac{\omega_{lin}(x, t)}{k(x, t)} \mathbf{n} \qquad (2.26)$$

$$\omega_{lin}^2 = gk \tanh kh \qquad (2.27)$$

$$\text{with} \quad \begin{cases} k(x, t) = \frac{\partial \varphi}{\partial t}(x, t) \\ \mathbf{n}(x, t) = \frac{\varphi(x, t)}{k(x, t)} \end{cases} \qquad (2.28)$$

where $h$ is the still water depth and $\mathbf{n}(x, t)$ represents the local wave propagation direction. The local wavenumber vector $k(x, t)$ can be explicitly deduced from the analytic expression of the mixed local phase function $\varphi_{12}(x, t)$ in Eq. (2.22). The derivation and analytic expression of $k(x, t)$ is provided in Appendix A.

**Spatial Hilbert Transform Method (SHTM):** The local wavenumber vector $k(x, t)$ and angular frequency $\omega(x, t)$ are directly computed by applying the $k-\omega$ conditions to the local phase function $\varphi(x, t)$. Then the phase speed is computed as follows:

$$c_{SHTM}(x, t) = \frac{\omega(x, t)}{k(x, t)} \mathbf{n}(x, t) \qquad (2.29)$$

$$\text{with} \quad \begin{cases} \omega(x, t) = - \frac{\partial \varphi}{\partial t}(x, t) \\ k(x, t) = \frac{\partial \varphi}{\partial n}(x, t) \\ \mathbf{n}(x, t) = \frac{\nabla \varphi(x, t)}{k(x, t)} \end{cases} \qquad (2.30)$$

In Eq. (2.30), the temporal derivative of the mixed local phase function is involved in evaluating the angular frequency $\omega(x, t)$. In the present work, it is numerically computed with a fourth-order backward Lagrange difference scheme as described in Appendix B. Similarly to LSHTM, the wavenumber vector $k(x, t)$ is computed analytically as described in Appendix A.

### 2.2. Crest tracking method

Lastly, the methodology adopted in the present work for the Crest Tracking Method (CTM) in short-crested seas is discussed. It will be used as a reference to characterize the properties of the proposed methodologies based on Spatial Hilbert Transforms. The CTM procedure is separated into two parts: first identifying the individual 3D waves, and second tracking the wave crest location in time.

For the identification of individual 3D waves in a field $\eta(x, t)$, two basic concepts are introduced. Firstly, for each $y$-location a wave-by-wave analysis along the $x$-direction allows the definition of a set of 2D (unidirectional) waves and their corresponding wave crests. Secondly, the 3D (short-crested) waves are assembled by a set of 2D waves in the transverse $y$-direction based on a dedicated assembled condition detailed afterwards. As an illustration, the red dashed lines in Fig. 1(a) represent a set of 2D wave slices, and Fig. 1(b) displays wave crest lines of a set of assembled 3D waves, where the green lines represent the crest lines of the assembled 3D waves. Further details on the identification of the 3D

wave, which is complex from the numerical implementation view, are detailed in Appendix C.

Identifying the 3D waves is the first step to studying short-crested seas, which makes it possible to temporally track the 3D crest line evolution. If the location of the crest is given by $(x_c(t), y_c(t))$, correspondingly, the local crest phase speed components in the $x$ and $y$ direction are computed via taking the temporal derivative of the crest location.

$$c_x = \frac{dx_c}{dt}, \quad c_y = \frac{dy_c}{dt} \qquad (2.31)$$

Subsequently, taking the modulus of the local wave phase speed $c$ gives the normal wave celerity.

$$c_{CTM}(x_c, t) = \sqrt{(c_x)^2 + (c_y)^2} = \sqrt{(dx_c/dt)^2 + (dy_c/dt)^2} \qquad (2.32)$$

Numerically, we employ the least square method in time to perform the temporal derivative of the successive individual crest locations to evaluate the local crest phase speed. In the least square fitting approach, it is necessary to relate the information of wave crest locations from one time step to the other, so the concepts of mother wave and wave family are introduced. For simplicity, a 2D wave is used to demonstrate those concepts. As shown in Fig. 2(a), if the wave crest location of wave $i-1$, $x_c(t_{i-1})$ is located in the range of $(x_u(t_i), x_c(t_i))$, wave $i-1$ (the notation is exactly the same as the time instant index) is regarded as the mother wave of wave $i$. Besides, if more than one wave satisfies this condition, we choose the one closest to $x_c(t)$. The wave family, as illustrated in Fig. 2(b), is nothing but the temporal series of the waves that satisfy the condition of the mother wave at a series of consecutive temporal points. The temporal series of the wave crest locations, therefore, are collected as

$$\{x_c(t_{i-n}), \dots, x_c(t_{i-1}), x_c(t_i)\}$$

Thus, the local crest phase speed component in $x$ is calculated by performing the least square fitting in time. In the $y$-direction, the same least square fitting is applied. Additionally, to determine the number of wave family members we used to perform the least squares fitting, a time threshold $T_{th}$, is introduced to define the minimum time window for the least squares fitting. This defines the minimum duration of an individual wave before computing its crest phase speed. This is taken as $T_{th} = T_c/16$, with $T_c$ the characteristic period of the wavefield studied.

## 3. Implementation in nonlinear wave models

To assess the different crest phase speed methodologies, nonlinear wave models are used for the description of the spatio-temporal free surface elevation $\eta(x, t)$. It is of particular interest in such models to have an efficient and accurate way of estimating the local crest phase speed for the development of a breaking model in short-crested seas to dissipate energy loss due to wave breaking. However, the proposed methodologies are general and can be used in other configurations. The nonlinear wave models we used are two open-source softwares, HOS-ocean<fnref n="1" /> for open ocean and HOS-NWT<fnref n="2" /> for Numerical Wave Tank (NWT) (Ducrozet, 2007; Ducrozet et al., 2012, 2016), released by the Laboratory in Hydrodynamics, Energetics and Atmospheric Environment (LHEEA) of École Centrale de Nantes and CNRS. More details of these two High-Order Spectral (HOS) models can be found in the previously cited articles. Here, we will briefly introduce them.

Assuming that the fluid is inviscid, incompressible, and irrotational, the velocity potential $\phi(x, z, t)$ can be introduced to satisfy the governing Laplace equation. Following (Zakharov, 1968), a surface velocity potential quantity $\phi_s(x, z = \eta(x, t), t)$ is introduced, leading to the following kinematic and dynamic free surface conditions for HOS-ocean

$$\frac{\partial \eta}{\partial t} = (1 + |\nabla \eta|^2) W - \nabla \phi_s \cdot \nabla \eta \qquad (3.1)$$

<footnote n="1">¹ https://gitlab.com/lheea/HOS-Ocean</footnote>

<footnote n="2">² https://gitlab.com/lheea/HOS-NWT</footnote>
<figure><img src="image_4.png" /><figcaption>Fig. 1. Schematic view of 2D waves in the longitudinal direction and the wave crest lines of assembled 3D waves. (a) The red dashed lines are the longitudinal slices, consisting of several 2D waves. (b) The green lines are the wave crest lines and the green circles are the crests of the unidirectional waves (i.e. the longitudinal slice of the directional wave field). In the longitudinal direction, every 3D wave only contains one 2D wave, but in the transversal direction, it is composed of a cluster of 2D waves.</figcaption></figure>

$$ \frac{\partial \phi_s}{\partial t} = -g\eta - \frac{1}{2} |\nabla \phi_s|^2 + \frac{1}{2} (1 + |\nabla \eta|^2) W^2 \quad (3.2) $$

expressed on $z = \eta(x, t)$ with $\nabla = (\partial/\partial x, \partial/\partial y)$ the horizontal gradient. In those equations $W(x, t) = \frac{\partial \phi}{\partial x}(x, \eta, t)$ is the vertical velocity at the free surface which can be evaluated by employing the so-called HOS scheme. The lateral boundary is assumed to be periodic in both the x- and y-direction, and various initial conditions-from regular wave to directional irregular sea states-are accessible within HOS-ocean (Ducrozet et al., 2017). Once the initial condition is specified and $W$ obtained, these two surface quantities $(\eta, \phi_s)$ can be time-marched by adopting the 4th order Runge-Kutta Cash-Karp scheme with adaptive step size (Cash and Karp, 1990). In the HOS-NWT model, the free surface conditions are similar to Eqs. (3.1) and (3.2) but an additional potential is introduced to take into account the presence of the wavemaker and the associated boundary condition. In this NWT, all other side walls are

treated as fixed, and consequently opposite to the wavemaker, a numerical beach is necessary to prevent wave reflections. All the physical features of an experimental wave tank are reproduced in this NWT, allowing for a straightforward comparison to dedicated experiments. Note that with HOS-NWT, the initial conditions correspond to a wave tank at rest, the waves being generated by the prescribed wavemaker motion.

The HOS models use a pseudo-spectral formalism and the discretization is consequently qualified through the largest wavenumber in Fourier space, $k_{xmax} = N_x \pi / L_x$, $k_{ymax} = N_y \pi / L_y$ for HOS-NWT and $k_{xmax} = (N_x - 1) \pi / L_x$, $k_{ymax} = (N_y - 1) \pi / L_y$ for HOS-ocean, in terms of the number of grid points ($N_x, N_y$) and size of the computational domain ($L_x, L_y$) in x and y-direction, respectively. Complementary, the Nodes Per characteristic Wavelength (NPW) is introduced to estimate the density of the spatial discretization. The accuracy of the time-stepping scheme is
<figure><img src="image_5.png" /><figcaption>Fig. 2. Demonstration of the definition of the mother wave and the wave family in 2D. The suffix of the time points represents the name of the waves, the wave name at $t_{i-1}$ is $i-1$, for instance. The red dots are the wave crests, $x_c(t_i)$, $x_u(t_i)$ and $x_d(t_i)$ are the wave crest locations, zero up- and down-crossing points of wave $i$, the index $n$ is the number of the members of the wave family.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 3. Visualization of the contour and the crest phase speed of the linear regular wave. The solid green lines are the wave crest line, and the red arrows are the wave phase speed at the wave crest which is computed by the SHTM.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 4. Comparison of CTM, LSHTM, SHTM, as well as the analytical solution in linear regular wave configuration with $\beta = 0°$. The crest phase speed is output from a single wave family, and the original HOS grid, NPW = 21, is used for computation. The Refined CTM results are obtained with NPW = 2133 on the refined grid.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Fig. 5. Comparison of CTM, LSHTM, SHTM as well as the analytical solution in a linear regular wave configured with $\beta = 45^{\circ}$. The crest phase speed is output from a single wave family, and the original HOS grid, NPW = 11 in the direction of wave propagation, is used for computation.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 6. Comparison of CTM, LSHTM, SHTM as well as the reference solution in the nonlinear regular wave configuration with $\beta = 0^{\circ}$. The crest phase speed is output from a single wave family and the original HOS grid, NPW = 32, is used for computation. The Refined CTM results are obtained with NPW = 1600 on the refined grid.</figcaption></figure>

controlled by a tolerance toler, and the HOS order mHOS governs the nonlinearity and accuracy of the HOS schemes.

HOS pseudo-spectral models usually need a limited spatial discretization for accurate wave propagation, yet this is not adequate for CTM to find a satisfactory result. Then, to make the free surface smoother, a refined grid system, represented by $N_{xref}$ and $N_{yref}$ in x and y-direction, can be easily introduced by adding zero-padding points in the Fourier space under the framework of the HOS model (Li et al., 2021). Note that such a refined grid system only aims to obtain a smoother spatial description of the surface quantities (i.e. it does not change the numerical solution within HOS). The implementation of SHT-based methodologies is straightforward within HOS models in which the estimation of the Hilbert transforms is readily available thanks to the pseudo-spectral description of the surface quantities. However, the temporal derivative, involved in SHTM, needs to store surface elevation in several time steps, whereas the LSHTM does not.

## 4. Results

In this section, four different wave fields are investigated in detail to evaluate the effect of the directionality, nonlinearity as well as bandwidth of the wave spectrum on the different approaches used to compute the local crest phase speed described in Section 2. Firstly, linear and non-linear regular waves, configured with propagating angles of $\beta = 0^{\circ}$ and $\beta = 45^{\circ}$, are simulated within the HOS-ocean solver. Frequency focusing waves, configured in HOS-NWT, are separately investigated, followed by combined frequency and directional focusing waves. The accuracy of estimation of the crest phase speed is estimated by the relative error defined as

$$ \text{Relative error} = \left| \frac{c_n^{\text{num}} - c_n^{\text{ref}}}{c_n^{\text{ref}}} \right| \qquad (4.1) $$

where the $c_n^{\text{ref}}$ is the reference value.

<figure><img src="image_10.png" /><figcaption>Fig. 7. Comparison of CTM, LSHTM, SHTM as well as the reference solution of the nonlinear regular wave configured with $\beta = 45^{\circ}$, the crest phase speed is output from a single wave family, and the original HOS grid, NPW = 32, is used for computation.</figcaption></figure>
<table><caption>**Table 1**
Configuration of the linear regular waves.</caption><thead><tr><th>β</th><th>a (m)</th><th>k (m<sup>-1</sup>)</th><th>f (Hz)</th><th>c<sup>ana</sup><sub>n</sub> (m/s)</th><th>mHOS</th><th>toler</th><th>NPW</th></tr></thead><tbody><tr><td>0°</td><td>0.02</td><td>4.7124</td><td>1.0821</td><td>1.4428</td><td>1</td><td>10<sup>-5</sup></td><td>21</td></tr><tr><td>45°</td><td>0.02</td><td>6.6643</td><td>1.2869</td><td>1.2132</td><td>1</td><td>10<sup>-5</sup></td><td>11</td></tr></tbody></table>

## 4.1. Linear regular wave

As an initial verification, two linear regular waves are studied, whose characteristics were provided in Table 1 with *a* the wave amplitude, *k* the wavenumber (in the direction of wave propagation), *f* the frequency, $c_n^{ana}$ the analytical value of the phase speed which will be used as the reference, and mHOS the HOS order. The first one propagates along the *x*-direction with an angle of $\beta = 0°$, and the second one with a propagating angle of $\beta = 45°$, as shown in Fig. 3(a,b). The length and width of the computational domain are set as $L_x = 8.0$ m, and $L_y = 8.0$ m, and the water depth is set as infinite. This ensures that the free surface elevation exhibits periodicity in *x* and *y* directions for both configurations.

For the propagating angle of $\beta = 0°$, a comparison of the crest phase speed obtained by CTM, LSHTM, SHTM as well as the analytical solution, is presented in Fig. 4. The results from LSHTM and SHTM are superposed in the first figure. Note that all time plots using SHTM are obtained by evaluating the quantities at successive time steps. In this unidirectional wave configuration, the crest phase speed is output from

a single crest propagating throughout the whole domain, i.e. a single wave family. Compared to the analytical solution, the results of LSHTM and SHTM, obtained with the original HOS discretization (NPW = 21), show very good agreement. In detail, the LSHTM and SHTM overlap with the theoretical value as shown in Fig. 4(a), and the relative errors are less than $10^{-7}$ as displayed in Fig. 4(b). The equivalence of LSHTM and SHTM is expected in this linear configuration. Conversely, large discrepancies between theoretical and CTM results are observed. This is attributed to the relatively coarse mesh adopted to represent the free surface in the HOS simulations. Indeed, the evaluation of the crest phase speed with CTM has the ability to converge toward the analytical solution if a refined grid is used, see the refined CTM results obtained with NPW = 2133 in Fig. 4 from which relative errors less than $10^{-2}$ are accomplished. This is one possible limitation of the CTM, which exacerbates the computational effort for practical applications.

In the same manner as $\beta = 0°$, for the linear configuration of the oblique regular wave, $\beta = 45°$ for the current case, a comparison of the crest phase speed, computed by the aforementioned approaches and analytical solution is drawn in Fig. 5 where perfect agreement (overlapped) is observed for all except CTM. The oscillations appearing in the CTM are for the same reason as the previous case $\beta = 0°$ which we have discussed. Additionally, it is interesting to note that the LSHTM and SHTM can compute the wave phase speed in the whole domain (and not only at the crest location). Fig. 3 displays the wave phase speed at the crest computed by the SHTM.

It is interesting to present a more detailed discussion on why the SHTM and LSHTM achieve perfect agreement in the computation of the

<figure><img src="image_11.png" /><figcaption>**Fig. 8.** Illustration of the complex signal $\Psi_1$ at the theoretical focusing time $t_f$ and location $x_f$ for case FFW-1.</figcaption></figure>
linear regular wave phase speed. Such wave propagating with $\beta$ angle can be expressed as

$$ \eta_{LRW}(\mathbf{x}, t) = a \cos[k(x \cos \beta + y \sin \beta) - \omega t] \quad (4.2) $$

From Eq. (4.2), the analytical normal phase speed is given by $c_n^{ana} = \omega/k$. Computing SHTs with Eq. (4.2) gives

$$ \mathbb{H}_x[\eta_{LRW}(\mathbf{x}, t)] = a \sin[k(x \cos \beta + y \sin \beta) - \omega t] \quad (4.3) $$

$$ \mathbb{H}_y[\eta_{LRW}(\mathbf{x}, t)] = a \sin[k(x \cos \beta + y \sin \beta) - \omega t] \quad (4.4) $$

$$ \mathbb{H}_{xy}[\eta_{LRW}(\mathbf{x}, t)] = -a \cos[k(x \cos \beta + y \sin \beta) - \omega t] \quad (4.5) $$

Inserting Eqs. (4.2)–(4.4) into Eqs. (2.9) and (2.10) yields

$$ \Psi_1(\mathbf{x}, t) = 2a \cos[k(x \cos \beta + y \sin \beta) - \omega t] + j(2a \sin[k(x \cos \beta + y \sin \beta) - \omega t]) $$

$$ \Psi_2(\mathbf{x}, t) = 0 $$

From those two complex signals, it is possible to construct the local amplitude and phase components using Eqs. (2.13)–(2.16) as

$$ A_1(\mathbf{x}, t) = 2a \qquad ; \quad A_2(\mathbf{x}, t) = 0 \\ \varphi_1(\mathbf{x}, t) = k(x \cos \beta + y \sin \beta) - \omega t ; \quad \varphi_2(\mathbf{x}, t) = k(x \cos \beta + y \sin \beta) - \omega t $$

Then, the mixed local amplitude and phase functions read

$$ A_{12}(\mathbf{x}, t) = a \quad (4.6) $$

$$ \varphi_{12}(\mathbf{x}, t) = k(x \cos \beta + y \sin \beta) - \omega t \quad (4.7) $$

From the definition of those mixed functions, the SHT reconstructed surface elevation yields

$$ \eta_{SHT}(\mathbf{x}, t) = a \cos[k(x \cos \beta + y \sin \beta) - \omega t] $$

which is rigorously equivalent to the original surface elevation we have given in Eq. (4.2). Substituting the mixed local phase function Eq. (4.7) into Eqs. (2.26)–(2.28), we find that the LSHTM-computed normal phase speed is $c_{LSHTM} = \omega/k$, which strictly corresponds to the analytical solution. Correspondingly, if we insert Eq. (4.7) into Eqs. (2.29) and (2.30), the SHTM-computed normal phase speed is $c_{SHTM} = \omega/k$, which is also rigorously the analytical solution.

## 4.2. Nonlinear regular wave

To investigate how the nonlinearity of the wave impacts the behaviours of the aforementioned methodologies, two nonlinear regular waves, configured with propagating angles of $\beta = 0^\circ$ and $\beta = 45^\circ$, were simulated within the HOS-ocean solver. As a comparison, the exact solution of the studied nonlinear regular waves is given by the stream function theory, which allows us to describe steep nonlinear regular waves. More details about the numerical implementation of the stream function theory used in this work (CN-Stream open-source model) can be found in Ducrozet et al. (2019).

The computational domain for $\beta = 0^\circ$ is set as $L_x = 8.0\lambda_0$ and $L_y = 8.0\lambda_0$ in the longitudinal and transverse directions respectively, where $\lambda_0 = 6.28 \text{ m}$ ($k_0 = 1.0 \text{ m}^{-1}$) is the wavelength in the direction of wave propagation. The water depth is set as infinite. Regarding the nonlinear oblique wave at $\beta = 45^\circ$, the wavelength is kept the same and the size of the computational domain is adjusted to ensure a periodicity in both x and y directions. The detailed configurations are listed in Table 2, in which $f_0$ is the frequency of the wave component.

A comparison between CTM, LSHTM, SHTM and reference results was performed and is displayed in Fig. 6(a) for $\beta = 0^\circ$. As expected,

<figure><img src="image_12.png" /><figcaption>Fig. 9. Illustration of the complex signal $\Psi_1$ at the theoretical focusing time $t_f$ and location $x_f$ for case FFW-3.</figcaption></figure>
<table><caption>**Table 2**
Configuration of the nonlinear regular waves.</caption><thead><tr><th>β</th><th>a (m)</th><th>k<sub>0</sub> (m<sup>-1</sup>)</th><th>f<sub>0</sub> (Hz)</th><th>ε<sub>n</sub><sup>ref</sup> (m/s)</th><th>mHOS</th><th>toler</th><th>NPW</th></tr></thead><tbody><tr><td>0°</td><td>0.1</td><td>1.0</td><td>0.501</td><td>3.1478</td><td>7</td><td>10<sup>-8</sup></td><td>32</td></tr><tr><td>45°</td><td>0.1</td><td>1.0</td><td>0.501</td><td>3.1478</td><td>7</td><td>10<sup>-8</sup></td><td>32</td></tr></tbody></table>

the LSHTM is not able to capture the nonlinear information of the crest phase speed since it employs the linear dispersion relation to evaluate the angular frequency, inducing a noticeable discrepancy. On the contrary, the SHTM shows remarkable agreement with the reference results, overlapping with the reference value as shown in Fig. 6(a) and with maximum relative errors of less than 0.02% as shown in Fig. 6(b). Similarly to the linear case, large discrepancies between CTM and reference solution are seen when using the computational grid (NPW = 32) during the solution of the wave propagation. This can be solved by using a refined grid (NPW = 1600) as shown in Fig. 6(a, b). After configuring a propagating angle, the SHTM, LSHTM and CTM show the same conclusion as $β = 0°$, and a comparison is drawn in Fig. 7.

### 4.3. Frequency focusing wave

To see how the bandwidth of the wave spectrum and the nonlinearity, represented by the wave steepness, influence the local crest phase speed evaluation, the uni-directional frequency focusing waves studied

by Seiffert et al. (2017) and Seiffert and Ducrozet (2018) with HOS-NWT were revisited here. The theoretical linear free surface elevation of the wave is a superposition of different sinusoidal elementary solutions based on the linear theory.

$$ \eta_{lin}(x, t) = \sum_{n=1}^{N_f} a_n \cos(k_n x - \omega_n t - \theta_n) \quad (4.8) $$

where $a_n$ is the amplitude, $k_n$ represents the wavenumber, $\omega_n$ denotes the angular frequency defined by the frequency $f_n$ (i.e. $\omega_n = 2\pi f_n$), $\theta_n$ is the phase of each sinusoidal wave component, $N_f$ is the number of wave components, and the steepness of each wave component is given by $\epsilon_n = k_n a_n$. The surface elevation at the wavemaker is expressed as follows to achieve energy focusing at a specific location $x_f$ and time $t_f$

$$ \eta(0, t) = \sum_{n=1}^{N_f} a_n \cos[-\omega_n t - (k_n x_f - \omega_n t_f)] \quad (4.9) $$

Using Eq. (4.9) with the adequate wavemaker transfer function ensures the definition of the wavemaker motion that will provide the focusing at location $x_f$ and time $t_f$ in the wave tank. The focusing is expected to be perfect in a linear framework, while nonlinearities are expected to produce a shift of the focusing time and location. Following (Seiffert et al., 2017), the computational domain is set as $L_x = 100$ m, $h = 2.9$ m, configured by a constant steepness $\epsilon_n$, a central frequency $f_c$, a frequency band $\delta f$ as shown in Table 3, and the focusing location and time as

<figure><img src="image_13.png" /><figcaption>**Fig. 10.** Reconstruction of the original surface elevation signal on the theoretical focusing time $t_f$ and location $x_f$ for case FFW-1.</figcaption></figure>
<table><caption>**Table 3**
Configuration of the frequency focusing waves.</caption><thead><tr><th>Cases</th><th>f<sub>c</sub> (Hz)</th><th>δf (Hz)</th><th>ε<sub>n</sub></th><th>mHOS</th><th>NPW</th></tr></thead><tbody><tr><td>FFW-1</td><td>0.6</td><td>0.54</td><td>0.0048</td><td>1</td><td>44</td></tr><tr><td>FFW-2</td><td>0.6</td><td>1.08</td><td>0.0048</td><td>1</td><td>44</td></tr><tr><td>FFW-3</td><td>0.6</td><td>0.54</td><td>0.0048</td><td>5</td><td>44</td></tr><tr><td>FFW-4</td><td>0.6</td><td>0.54</td><td>0.0051</td><td>5</td><td>44</td></tr></tbody></table>

$(x_f, t_f) = (22.0 \text{ m}, 50.0 \text{ s})$ are specified. Cases FFW-1 and FFW-2 assume a linear solver with different frequency bandwidths, while cases FFW-3 and FFW-4 include nonlinearities in the Numerical Wave Tank with wave fields of different steepnesses.

For uni-directional seas, Eqs. (2.9) and (2.10) are strictly equivalent since $H_y[\eta(x, t)]$ and $H_{xy}[\eta(x, t)]$ are zero. Note that for unidirectional configuration, the real part of $\Psi_1$ is simply the free surface elevation and the imaginary part is its SHT. The complex signals $\Psi_1$ are displayed in Fig. 8 for FFW-1 and Fig. 9 for FFW-3 as an illustration. As a preliminary validation, the free surface elevation obtained through the SHT expression (Eq. (2.20)), labeled $\eta_{SHT}(x, t)$, is compared to the original signal obtained from the numerical simulations $\eta(x, t)$. The comparison is achieved in both space and time for linear and nonlinear configurations, cases FFW-1 and FFW-3, as depicted in Figs. 10 and 11, respectively. The SHT expression perfectly reconstructs the original surface elevation signal for both linear and nonlinear setups, which validates the correctness of the SHT expression. From the proposed SHTM, one

can recover, in addition to the phase that will be used for the local crest phase speed evaluation, the envelope of the wave signal $A(x, t)$, which is also depicted in Figs. 10 and 11(a).

Since multiple frequency components are contained in the current irregular wave fields, finding the local crest phase speed's analytical solution poses a great challenge, particularly when nonlinearity is included. Thus, the results from the CTM with a refined grid (NPW = 4443) will be considered as the reference, and the proposed methodologies based on Hilbert Transforms will be validated against this CTM solution.

A comparison of the approaches mentioned above is drawn in Fig. 12(a) for the linear configured case FFW-1, in which the SHTM overestimates the local crest phase speed, and the LSHTM overestimates it even more. Quantitatively, the relative error of SHTM, estimated by Eq. (4.1) is 13.6% at the real focusing location and time, and that of the LSHTM is 19.1%. The discrepancies between SHTMs and the reference CTM are primarily attributed to the effect of the spectrum bandwidth, as discussed in Section 2.1.3, which is an intrinsic limitation of HTMs, and one that would not be improved by employing a refined grid. Since the $k-\omega$ conditions are only satisfied when the slow-varying assumption of the local wavenumber and angular frequency is made, the accuracy of the SHTM and LSHTM will decrease when the spectrum bandwidth starts to play a role. In addition, there is one point that needs to be further clarified in Fig. 12(a): at first glance, the LSHTM results are close to the SHTM results for the linear irregular wave investigated, but this is misleading. As shown in Section 2.1.4, the difference between LSHTM

<figure><img src="image_14.png" /><figcaption>**Fig. 11.** Reconstruction of the original surface elevation signal on the theoretical focusing time $t_f$ and location $x_f$ for case FFW-3.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>Fig. 12. Local crest phase speed, computed by CTM, LSHTM, SHTM for all cases listed in Table 3; $t'_f$ represents the real focusing time.</figcaption></figure>

and SHTM is that they adopt different ways to compute the local angular frequency after obtaining the same single local wave number. SHTM directly uses the $\omega$-condition on the local phase function to evaluate the local angular frequency, while LSHTM uses the linear dispersion relation on this local wavenumber. The SHT-based methodologies define a single characteristic local wavenumber, but its content is physically made of different wavelengths. This means that using the linear dispersion relation will be less accurate when the frequency content broadens

and energy spreads toward values different from the characteristic local wavenumber.

The case FFW-2, depicted in Fig. 12(b), exhibits a bandwidth that is doubled compared to case FFW-1 ($\delta f = 1.08$ Hz). This leads to larger discrepancies for both SHTM and LSHTM, the latter being more impacted. Thus, for the linear configured irregular wave, the bandwidth of the wave spectrum plays a fundamental role in the accuracy of the SHTM and LSHTM, i.e. the accuracy of the SHTM will decrease

<figure><img src="image_16.png" /><figcaption>Fig. 13. Illustration of the complex signals, $\Psi_1$ and $\Psi_2$ for the nonlinear HOS configuration. Top and middle: the contour of real and imaginary parts of the complex signal (a) $\Psi_1$ and (b) $\Psi_2$ on the theoretical focusing time $t_f$. Bottom: the temporal variation of the complex signal (a) $\Psi_1$ and (b) $\Psi_2$ at the theoretical location $x_f$.</figcaption></figure>
when broadening the bandwidth, and the LSHTM will decrease even more.

Interestingly, as shown in Fig. 12(c), when we configure the same case as the highly nonlinear one, i.e. case FFW-3, an encouraging agreement is achieved for SHTM, whose relative errors are estimated to be around 6.7% at the real focusing location and time $t'_f$, whereas the relative errors of the LSHTM are 27.3%. This is mainly because of the effect

of nonlinearity. Further, if we increase the steepness from $\epsilon_n = 0.0048$ to $\epsilon_n = 0.0051$, higher accuracy is observed for the SHTM as presented in Fig. 12(d) in which the relative error at the real focusing location and time is 4.5%. However, such an improvement with the level of nonlinearity is not observed for LSHTM since it uses linear dispersion which fails to capture part of the nonlinear information. If we compare the results of cases FFW-1, FFW-2, and FFW-3 presented in Fig. 12(a,c,d), it can be

<figure><img src="image_17.png" /><figcaption>Fig. 14. Reconstruction of the original surface elevation signal on the theoretical focusing time $t_f$ and location $x_f$ for the nonlinear HOS configuration.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>Fig. 15. Illustration of the crest phase speed vector (red arrows) at the theoretical focusing time $t_f$, computed by SHTM. The green lines represent the crest line.</figcaption></figure>

concluded that the local crest phase speed appears to be driven more by the nonlinearity than by the frequency bandwidth, which results in an accurate estimate of local crest phase speed with SHTM for highly nonlinear waves. It should be pointed out, however, that the question of how the entangled bandwidth of the wave spectrum and nonlinearity quantitatively influence the local crest phase velocity evaluated by SHTM remains open.

## 4.4. Frequency and directional focusing wave

In this subsection, more general frequency and directional energy focusing waves will be considered within the HOS-NWT solver. Following Wu and Nepf (2002), the wavemaker motion is specified through the surface elevation at the wave paddles.

$$ \eta(0, y_l, t) = \sum_{n=1}^{N_f} a_n \cos \left[ -k_n(|y_l| \sin(\beta_l) + x_f \cos(\beta_l)) + \omega_n t_f \right] \quad (4.10) $$

where the propagation angle, $\beta_l$, is determined by

$$ \beta_l = \arctan \left( \frac{|y_l - y_f|}{x_f} \right) \quad (4.11) $$

and the amplitude of each component is determined by a so-called gain factor G

$$ a_n = G/k_n \quad (4.12) $$

The computational domain is given as $L_x = 11.0$ m, $L_y = 4.0$ m, $h = 0.6$ m, and the theoretical focusing location and time are fixed as $x_f = (3.3$ m, 2.0 m), $t_f = 12.0$ s. The number of wave components is set as $N_f = 32$, the gain factor G is defined through the global focusing steepness, $S = k_c \sum_{n=1}^{N_f} a_n$, as

$$ G = S \left( k_c \sum_{n=1}^{N_f} \frac{1}{k_n} \right) \quad (4.13) $$

and the frequencies are equispaced across a bandwidth of $\delta f$ with a central frequency $f_c$. The detailed configurations are presented in Table 4.

The complex signals, constructed by Eqs. (2.9) and (2.10), are displayed in Fig. 13 for an illustration. As a verification, employing Eq. (2.20) to reproduce the original surface elevation signal at the theoretical focusing time and location is initially checked for the nonlinear configured case FDFW-3. The original surface elevation is excellently reproduced at the theoretical focusing point, as depicted in Fig. 14(a)

<table><caption>Table 4
Configuration of the frequency and directional focusing waves.</caption><thead><tr><th>Cases</th><th>f<sub>c</sub> (Hz)</th><th>δf (Hz)</th><th>S</th><th>mHOS</th><th>NPW</th></tr></thead><tbody><tr><td>FDFW-1</td><td>1.08</td><td>0.789</td><td>0.21</td><td>1</td><td>61</td></tr><tr><td>FDFW-2</td><td>1.08</td><td>1.578</td><td>0.21</td><td>1</td><td>61</td></tr><tr><td>FDFW-3</td><td>1.08</td><td>0.789</td><td>0.21</td><td>7</td><td>61</td></tr><tr><td>FDFW-4</td><td>1.08</td><td>0.789</td><td>0.26</td><td>7</td><td>61</td></tr></tbody></table>

where the local amplitude stands for the wave envelope. This reconstruction at any time instant (here we only select the theoretical focusing time) holds in space, as shown in Fig. 14(b). The original surface elevation, plotted by the blue contour, overlaps with the surface computed with Eq. (2.20), represented by the red dashed wire-frame. Interestingly, the mixed local amplitude still retains its envelope property as presented in Fig. 14(b) in which the local mixed amplitude, plotted by the transparent green surface with dash-dotted lines, covers the wave surface. Fig. 15 displays the crest phase speed vector, calculated by the SHTM, at the wave crest line on the theoretical focusing time for case FDFW-3. Visually, the directionality is correct. Note that the SHTM can compute the phase velocity at each spatial point, displayed here at the crests for each y-transect to validate the direction evaluation.

A more detailed comparison of the local crest phase speed of the four cases listed in Table 4 is displayed in Fig. 16, where the duration plotted is a single whole wave family for CTM. For case FDFW-1, the LSHTM and SHTM show a comparable ability to capture the local crest phase speed and are close to the CTM results (which are obtained through a very fine refined grid, NPW = 3138); the relative errors at the real focusing time and location are 7.19% and 8.57%, respectively. As discussed in Section 4.3, the discrepancies between SHTM and the reference CTM result will be amplified if we broaden the bandwidth of the spectrum, i.e. case FDFW-2, as demonstrated in Fig. 16(b), and this discrepancy is more severe for LSHTM. When the nonlinearity grows with the numerical configuration within HOS-NWT, the discrepancies, led by the additional terms in Eqs. (2.23) and (2.24), do not show significant improvement in SHTM since the nonlinearity of the wave is not too high. If a breaking wave, as reported by Wu and Nepf (2002), is configured, i.e. FDFW-4, the SHTM achieves excellent precision, as shown in Fig. 16(d). This confirms the conclusion reached in Section 4.3 that the local crest phase speed is dominated more by the nonlinearity than by the frequency bandwidth. In detail, the relative error is 1.7% for SHTM, and 36.4% for LSHTM. Regarding case FDFW-4, one point that needs to be further clarified in Fig. 16(d) is that the numerical simulation stops after the real focusing time since the wave breaking occurs, which subsequently leads to numerical instability.
<figure><img src="image_19.png" /><figcaption>Fig. 16. Local crest phase speed, computed by CTM, LSHTM, SHTM for all cases listed in Table 4. $t'_f$ represents the real focusing time.</figcaption></figure>

## 5. Conclusion

The present work correlated the local phase function and amplitude components, obtained by the multi-dimensional Hilbert transform, to a mixed local amplitude and phase function, which builds a way to find the envelope of multi-directional waves from spatial information. Whitham's modulation theory is derived through the newly developed local amplitude and phase function in multi-directional sea states by introducing slow-varying assumptions to the local wavenumber and angular frequency. Two methodologies, LSHTM and SHTM, are proposed to compute the local crest phase speed in multi-directional sea states.

The proposed LSHTM and SHTM were validated through two linear regular waves with propagating angles $\beta = 0^\circ$ and $\beta = 45^\circ$. Then, two nonlinear regular waves with propagating angles $\beta = 0^\circ$ and $\beta = 45^\circ$ were studied, through which we illustrate that the proposed SHTM is accurate for the evaluation of the nonlinear crest phase speed, whereas the LSHTM loses some nonlinear information about the crest phase speed due to its reliance on the linear dispersion relation.

As for linear irregular waves, the bandwidth of the wave spectrum plays a negative role in the accuracy of these two HTMs, i.e. a wider spectrum bandwidth leads to lower precision since the additional terms in Eqs. (2.24) and (2.25) would not be negligible when the bandwidth is broad. However, according to our study, for extreme nonlinear irregular waves, the local crest phase speed appears to be driven more by the nonlinearity than by the frequency bandwidth, which results in an accurate estimate of the local crest phase speed with SHTM for highly nonlinear waves. From our interest in detecting wave breaking on large spatio-temporal scales (up to 100 km$^2$ in space, up to 1 h in time) for ocean engineering applications, the proposed SHTM allows us to benefit from its accuracy, efficiency, stability as well as applicability when compared to the conventional

CTM in multi-directional sea states. The advantages of SHTM are summarised as:

* **Computational effort:** the SHTM is shown to be more efficient than the CTM. In all cases studied in the present work, CTM requires to have a refined mesh in order to produce a smooth evaluation of the crest speed. In unidirectional sea states, the utilization of splines would not lead to very large computational burden, but it does in multi-directional sea states, particularly taking into account that small time steps also need to be adopted for CTM. Taking the linear regular wave as an example, the wave field is converged with the HOS model when the Nodes Per characteristic Wave-length (NPW) is equal to 21. SHTM manages to compute crest phase speed with high accuracy with this mesh. However, large oscillations are observed when using CTM with the same mesh, as shown in Fig. 4(a). To achieve an accurate estimate of the crest phase speed, an extremely refined grid (NPW = 2133), more than one hundred times denser than the original mesh, is necessary as shown in Fig. 4(b).

* **Applicability.** In complex multidirectional sea states, the proposed SHTM is straightforward to implement. However, tracking multiple crests in a directional wave field is very complex to implement from a numerical point of view.

* **Stability:** The SHTM is strictly built on the Spatial Hilbert Transform, which can be easily obtained by spatial FFTs, and no numerical tricks such as filters were necessary. On the contrary, if we do not use splines, which are computationally expensive, or filters, the CTM evaluation will lead to large oscillations.

Thus, our present work makes the evaluation of the breaking onset parameter of Barthelemy et al. (2018) practical, which paves the way for its utilization in ocean/coastal engineering practice.
## Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## CRediT authorship contribution statement

Ying Wang: Writing – original draft, Visualization, Validation, Software, Methodology, Investigation, Conceptualization; Guillaume Ducrozet: Writing – review & editing, Supervision, Software, Project administration, Methodology, Funding acquisition, Conceptualization

## Acknowledgement

The work presented in this paper was partly supported by WASANO project funded by the French National Research Agency (ANR) as part of the “Investissements d'Avenir” Programme (ANR-16-IDEX-0007). Ying Wang is funded by the Chinese Scholarship Council (No. 202206680015).

## Appendix A. Derivation of the local wavenumber

Computing the local wavenumber by an explicit derivative expression of the proposed mixed local phase function is deduced here. Since the procedure is quite complex, the numerator part and denominator part inside the brace of the arctan function of the mixed local phase function are separately rewritten as two functions $N(x, t)$ and $D(x, t)$ for simplicity, respectively.

$$ \varphi_{12}(\mathbf{x}, t) = \arctan \left( \frac{\mathcal{N}(\mathbf{x}, t)}{D(\mathbf{x}, t)} \right) \qquad (A.1) $$

$$ \mathcal{N}(\mathbf{x}, t) = A_1(\mathbf{x}, t) \sin(\varphi_1(\mathbf{x}, t)) + A_2(\mathbf{x}, t) \sin(\varphi_2(\mathbf{x}, t)) \qquad (A.2) $$

$$ D(\mathbf{x}, t) = A_1(\mathbf{x}, t) \cos(\varphi_1(\mathbf{x}, t)) + A_2(\mathbf{x}, t) \cos(\varphi_2(\mathbf{x}, t)) \qquad (A.3) $$

The local wavenumber is computed by taking the gradient of the mixed local phase function.

$$ \mathbf{k}(\mathbf{x}, t) = \nabla \varphi_{12}(\mathbf{x}, t) = \frac{\nabla \mathcal{N}(\mathbf{x}, t) D(\mathbf{x}, t) - \mathcal{N}(\mathbf{x}, t) \nabla D(\mathbf{x}, t)}{\mathcal{N}^2(\mathbf{x}, t) + D^2(\mathbf{x}, t)} \qquad (A.4) $$

with

$$ \nabla \mathcal{N}(\mathbf{x}, t) = \nabla A_1(\mathbf{x}, t) \sin(\varphi_1(\mathbf{x}, t)) A_1(\mathbf{x}, t) \cos(\varphi_1(\mathbf{x}, t)) \nabla \varphi_1(\mathbf{x}, t) + \nabla A_2(\mathbf{x}, t) \sin(\varphi_2(\mathbf{x}, t)) + A_2(\mathbf{x}, t) \cos(\varphi_2(\mathbf{x}, t)) \nabla \varphi_2(\mathbf{x}, t) \qquad (A.5) $$

$$ \nabla D(\mathbf{x}, t) = \nabla A_1(\mathbf{x}, t) \cos(\varphi_1(\mathbf{x}, t)) - A_1(\mathbf{x}, t) \sin(\varphi_1(\mathbf{x}, t)) \nabla \varphi_1(\mathbf{x}, t) + \nabla A_2(\mathbf{x}, t) \cos(\varphi_2(\mathbf{x}, t)) - A_2(\mathbf{x}, t) \sin(\varphi_2(\mathbf{x}, t)) \nabla \varphi_2(\mathbf{x}, t) \qquad (A.6) $$

where

$$ \cos(\varphi_1(\mathbf{x}, t)) = \frac{1}{A_1(\mathbf{x}, t)} (\eta(\mathbf{x}, t) - \mathbb{H}_{xy}[\eta(\mathbf{x}, t)]) \qquad (A.7) $$

$$ \sin(\varphi_1(\mathbf{x}, t)) = \frac{1}{A_1(\mathbf{x}, t)} (\mathbb{H}_x[\eta(\mathbf{x}, t)] + \mathbb{H}_y[\eta(\mathbf{x}, t)]) \qquad (A.8) $$

$$ \cos(\varphi_2(\mathbf{x}, t)) = \frac{1}{A_2(\mathbf{x}, t)} (\eta(\mathbf{x}, t) + \mathbb{H}_{xy}[\eta(\mathbf{x}, t)]) \qquad (A.9) $$

$$ \sin(\varphi_2(\mathbf{x}, t)) = \frac{1}{A_2(\mathbf{x}, t)} (\mathbb{H}_x[\eta(\mathbf{x}, t)] - \mathbb{H}_y[\eta(\mathbf{x}, t)]) \qquad (A.10) $$

and

$$ \nabla \varphi_1(\mathbf{x}, t) = \frac{1}{A_1^2(\mathbf{x}, t)} \left[ \begin{aligned} & (\nabla \mathbb{H}_x[\eta(\mathbf{x}, t)] + \nabla \mathbb{H}_y[\eta(\mathbf{x}, t)])(\eta(\mathbf{x}, t) - \mathbb{H}_{xy}[\eta(\mathbf{x}, t)]) - \\ & (\mathbb{H}_x[\eta(\mathbf{x}, t)] + \mathbb{H}_y[\eta(\mathbf{x}, t)])(\nabla \eta(\mathbf{x}, t) + \nabla \mathbb{H}_{xy}[\eta(\mathbf{x}, t)]) \end{aligned} \right] \qquad (A.11) $$

$$ \nabla \varphi_2(\mathbf{x}, t) = \frac{1}{A_2^2(\mathbf{x}, t)} \left[ \begin{aligned} & (\nabla \mathbb{H}_x[\eta(\mathbf{x}, t)] - \nabla \mathbb{H}_y[\eta(\mathbf{x}, t)])(\eta(\mathbf{x}, t) + \mathbb{H}_{xy}[\eta(\mathbf{x}, t)]) - \\ & (\mathbb{H}_x[\eta(\mathbf{x}, t)] - \mathbb{H}_y[\eta(\mathbf{x}, t)])(\nabla \eta(\mathbf{x}, t) + \nabla \mathbb{H}_{xy}[\eta(\mathbf{x}, t)]) \end{aligned} \right] \qquad (A.12) $$

$$ \nabla A_1(\mathbf{x}, t) = \frac{1}{A_1(\mathbf{x}, t)} \left[ (\eta(\mathbf{x}, t) - \mathbb{H}_{xy}[\eta(\mathbf{x}, t)])(\nabla \eta(\mathbf{x}, t) - \nabla \mathbb{H}_{xy}[\eta(\mathbf{x}, t)]) + (\mathbb{H}_x[\eta(\mathbf{x}, t)] + \mathbb{H}_y[\eta(\mathbf{x}, t)])(\nabla \mathbb{H}_x[\eta(\mathbf{x}, t)] + \nabla \mathbb{H}_y[\eta(\mathbf{x}, t)]) \right] \qquad (A.13) $$

$$ \nabla A_2(\mathbf{x}, t) = \frac{1}{A_2(\mathbf{x}, t)} \left[ (\eta(\mathbf{x}, t) + \mathbb{H}_{xy}[\eta(\mathbf{x}, t)])(\nabla \eta(\mathbf{x}, t) + \nabla \mathbb{H}_{xy}[\eta(\mathbf{x}, t)]) + (\mathbb{H}_x[\eta(\mathbf{x}, t)] - \mathbb{H}_y[\eta(\mathbf{x}, t)])(\nabla \mathbb{H}_x[\eta(\mathbf{x}, t)] - \nabla \mathbb{H}_y[\eta(\mathbf{x}, t)]) \right] \qquad (A.14) $$

## Appendix B. Derivation of the local angular frequency

The fourth-order Lagrange interpolation of the local phase function in time is given by

$$ \mathcal{L}_{\varphi}(\mathbf{x}, t) = \sum_{i=1}^{4} \varphi(\mathbf{x}, t_i) l_i(t) \qquad (B.1) $$

in which the spatial variables are treated as constant; $l_i(t)$ is the Lagrange basis polynomial defined as

$$ l_i(t) = \prod_{\substack{j=1 \\ j \neq i}}^{4} \frac{t - t_j}{t_i - t_j} \qquad (B.2) $$

Thus, the angular frequency is figured out by taking the temporal derivative of the Lagrange interpolation of the local phase function Eq. (B.1).

$$ \omega(\mathbf{x}, t) = \frac{\partial \mathcal{L}_{\varphi}(\mathbf{x}, t)}{\partial t} = \sum_{i=1}^{4} \varphi(\mathbf{x}, t_i) \frac{dl_i(t)}{dt} \qquad (B.3) $$

where the derivative of the Lagrange basis polynomial gives

$$ \frac{dl_i(t)}{dt} = \sum_{\substack{j=1 \\ j \neq i}}^{4} \left( \frac{1}{t_i - t_j} \prod_{\substack{k=1 \\ k \neq i \\ k \neq j}}^{4} \frac{t - t_k}{t_i - t_k} \right) \qquad (B.4) $$

## Appendix C. Distinguishing the 3D wave

To distinguish one 3D wave from another, two basic concepts are introduced. Firstly, it is assumed that every slice $i_2$ of the free surface in the longitudinal direction contains a set of 2D waves and their corresponding wave crests. Secondly, the 3D waves are assembled by a set of 2D waves in the transverse direction based on an assembled condition, which will be illustrated below. Therefore, the algorithm that distinguishes 3D waves can be divided into two parts: the first part involves finding all the 2D waves at every slice, while the second part entails assembling 3D waves using the identified 2D waves.

**Part 1:** Find all the waves and corresponding zero up-crossing and down-crossing points in every longitudinal slice $i_2$ (i.e. every profile with a specific refined grid number in the transversal direction as shown in Fig. 1(a)), then save the information of crest and corresponding zero up-crossing and down-crossing points. In the 2D slice, the wave was distinguished by the consecutive zero up-crossing and down-crossing points. The illustration of the algorithm is presented via a flowchart in Fig. C.17. Additionally, before assembling 3D waves, all 2D waves are initially marked as unselected. Each individual 2D wave carries information about a pair of consecutive up and down-crossing points, including the location and wave elevation of the wave crest.

**Part 2:** All the 2D waves and wave crests have been identified and saved in part 1. In this part, the identified 2D waves will be assembled into 3D waves via the assembled condition. As shown in Fig. C.18, the so-called assembled condition is geometrically defined as follows: the crest of wave $j_2$ in slice $i_{2jump}$ located in the range of the consecutive zero up and down-crossing points of the wave $j_{1jump}$ in slice $i_{1jump}$ - 1, i.e. the red point in slice $i_{2jump}$ located in the shadowed zone. To assemble the 3D waves, start with an unselected 2D wave $j_1$ at slice $i_2$, initialize $j_{1jump} = j_1$ and $i_{2jump} = i_2 + 1$, then iterate 2D waves $j_2$ at slice $i_{2jump}$
<figure><img src="image_20.png" /><figcaption>**Fig. C.17.** Flowchart illustrating the algorithm to find the 2D waves at every longitudinal slice $i_2$. $\eta_{i1}$ represents the wave elevation at point $i_1$ with a specific slice $i_2$, $N_{xref}$ and $N_{yref}$ are the refined grid point numbers in the x- and y-direction, $n_{crest}$ represents the 2D wave crest number (i.e. the 2D wave number) of slice $i_2$. Y and N denote yes and no, respectively.</figcaption></figure>

and save all 2D waves $j_2$ that satisfy the assembled condition. If there is more than one 2D wave crest that satisfies the assembled condition in slice $i_{2jump}$, only the one among all saved $j_2$ closest to $j_{1jump}$ at slice $i_{2jump} - 1$ is collected to construct the 3D wave, and the collected $j_2$ is naturally marked as selected, but the others are not. This algorithm for constructing 3D waves is demonstrated via a flowchart in Fig. C.19. It is

worth noting that the saved 3D wave, denoted as *iwv* (where *iwv* represents the number of assembled 3D waves), contains information about a sequence of consecutive up- and down-crossing points that establish the edges of the 3D wave *iwv*. Additionally, it includes details about the location and wave elevation of the 2D waves forming the 3D wave crest line of the 3D wave *iwv*.

<figure><img src="image_21.png" /><figcaption>**Fig. C.18.** Illustration of the condition of assembling 3D wave by 2D waves. The consecutive green points are the zero up-crossing and down-crossing points, and the shadowed area is the transversal extended region of the zero up-crossing and down-crossing points of the 2D wave $j_{1jump}$ in slice $i_{2jump} - 1$. The red point in slice $i_{2jump}$ is the wave crest of the 2D wave $j_2$.</figcaption></figure>
<figure><img src="image_22.png" /><figcaption>**Fig. C.19.** Flowchart of the algorithm of assembling 3D waves by 2D waves. The so-called slice loop is the loop of the refined grid point in the y-direction; the crest loop iterates the 2D wave in slice $i_2$; *Is selected* ? in the decision box means "Is the 2D wave that is being looped collected by a 3D wave?"; $iwv$ absorb 2D wave $j_2$ means that the 2D wave $j_2$ is selected to construct a 3D wave. The notations used in the diagram are listed: $iwv$ is the 3D wave numbering; $n_{crest1}$ and $n_{crest2}$ are the number of 2D waves in slices $i_2$ and $i_{2jump}$; $N_{yref}$ denotes the number of refined grid points in the y-direction; $i_{2jump}$ and $j_{1jump}$ are the temporary indices that are used to shift slice and 2D wave along the transversal direction; $n_s$ represents the number of 2D waves that satisfy the assembled condition in slice $i_{2jump}$. If there is more than one 2D wave crest that satisfies the assembled condition in slice $i_{2jump}$, only the one among all saved $j_2$ closest to $j_{1jump}$ is collected to construct the 3D wave. It is naturally marked as selected, but the others are not. The saved 3D wave $iwv$ contains the information of a sequence of consecutive up and down-crossing points, forming the edges of $iwv$, as well as information about a set of the location and wave elevation of the 2D waves, creating the 3D wave crest line of $iwv$.</figcaption></figure>

## References

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. J. Fluid Mech. 841, 463–488.

Cash, J.R., Karp, A.H., 1990. A variable order Runge-Kutta method for initial value problems with rapidly varying right-hand sides. ACM Trans. Math. Softw. (TOMS) 16 (3), 201–222.

Chavez, M., Besserve, M., Adam, C., Martinerie, J., 2006. Towards a proper estimation of phase synchronization from time series. J. Neurosci. Methods 154 (1–2), 149–160.

Chu, V.H., Mei, C.C., 1970. On slowly-varying stokes waves. J. Fluid Mech. 41 (4), 873–887.

Craciunescu, C.C., Christou, M., 2020. On the calculation of wavenumber from measured time traces. Appl. Ocean Res. 98, 102115.

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. J. Geophys. Res. 125 (7), e2019JC015886.

Ducrozet, G., 2007. Modélisation des processus non-linéaires de génération et de propagation d'états de mer par une approche spectrale. Ph.D. thesis. Université de Nantes; Ecole Centrale de Nantes (ECN).

Ducrozet, G., Bonnefoy, F., Le Touzé, D., Ferrant, P., 2012. A modified high-order spectral method for wavemaker modeling in a numerical wave tank. Eur. J. Mech. - B/Fluids 34, 19–34.

Ducrozet, G., Bonnefoy, F., Le Touzé, D., Ferrant, P., 2016. Hos-ocean: open-source solver for nonlinear waves in open ocean based on high-order spectral method. Comput. Phys. Commun. 203, 245–254.

Ducrozet, G., Bonnefoy, F., Perignon, Y., 2017. Applicability and limitations of highly non-linear potential flow solvers in the context of water waves. Ocean Eng. 142, 233–244.

Ducrozet, G., Bouscasse, B., Gouin, M., Ferrant, P., Bonnefoy, F., 2019. Cn-stream: Open-source library for nonlinear regular waves using stream function theory. arXiv preprint arXiv:1901.10577

El, G.A., Gammal, A., Khamis, E.G., Kraenkel, R.A., Kamchatnov, A.M., 2007. Theory of optical dispersive shock waves in photorefractive media. Phys. Rev. A 76 (5), 053813.

Esler, J.G., Pearce, J.D., 2011. Dispersive dam-break and lock-exchange flows in a two-layer fluid. J. Fluid Mech. 667, 555–585.
Feldman, M., 2011. Hilbert transform in vibration analysis. Mech. Syst. Signal Process. 25 (3), 735–802.

Gabor, D., 1946. Theory of communication. Part 1: the analysis of information. J. Inst. Electr. Eng. - Part III 93 (26), 429–441.

Hilbert, D., 1905. Über eine Anwendung der integralgleichungen auf ein Problem der funktionentheorie.

Huang, N.E., Shen, Z., Long, S.R., 1999. A new view of nonlinear water waves: the Hilbert spectrum. Annu. Rev. Fluid Mech. 31 (1), 417–457.

Huang, N.E., Shen, Z., Long, S.R., Wu, M.C., Shih, H.H., Zheng, Q., Yen, N.-C., Tung, C.C., Liu, H.H., 1998. The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis. Proc. R. Soc. Lond. Ser. A 454 (1971), 903–995.

Hwang, P.A., Xu, D., Wu, J., 1989. Breaking of wind-generated waves: measurements and characteristics. J. Fluid Mech. 202, 177–200.

Krechetnikov, R., 2024. Hilbert transform view of water-wave theory. J. Fluid Mech. 986, R3.

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

Li, Z., Bouscasse, B., Ducrozet, G., Gentaz, L., Le Touzé, D., Ferrant, P., 2021. Spectral wave explicit Navier–Stokes equations for wave-structure interactions using two-phase computational fluid dynamics solvers. Ocean Eng. 221, 108513.

Lighthill, J., 1978. Waves in Fluids. Cambridge University Press.

Melville, W.K., 1983. Wave modulation and breakdown. J. Fluid Mech. 128, 489–506.

Minzoni, A.A., Smyth, N.F., 2016. Modulation theory, dispersive shock waves and gerald beresford whitham. Physica D 333, 6–10.

Mohanlal, S., Harris, J.C., Yates, M.L., Grilli, S.T., 2023. Unified depth-limited wave breaking detection and dissipation in fully nonlinear potential flow models. Coastal Eng. 183, 104316.

Perlin, M., Choi, W., Tian, Z., 2013. Breaking waves in deep and intermediate waters. Annu. Rev. Fluid Mech. 45, 115–145.

Porter, A., Smyth, N.F., 2002. Modelling the morning glory of the gulf of carpentaria. J. Fluid Mech. 454, 1–20.

Potamianos, A., Maragos, P., 1994. A comparison of the energy operator and the Hilbert transform approach to signal and speech demodulation. Signal Process. 37 (1), 95–120.

Poularikas, A.D., Grigoryan, A.M., 2018. Transforms and Applications Handbook. CRC Press.

Raut, R.V., Snyder, A.Z., Mitra, A., Yellin, D., Fujii, N., Malach, R., Raichle, M.E., 2021. Global waves synchronize the brain’s functional systems with fluctuating arousal. Sci. Adv. 7 (30), eabf2709.

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. J. Fluid Mech. 811, 642–658.

Schäfer, C., Rosenblum, M.G., Kurths, J., Abel, H.-H., 1998. Heartbeat synchronized with ventilation. Nature 392 (6673), 239–240.

Seiffert, B.R., Ducrozet, G., 2018. Simulation of breaking waves using the high-order spectral method with laboratory experiments: wave-breaking energy dissipation. Ocean Dyn. 68 (1), 65–89.

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.

Smyth, N.F., Holloway, P.E., 1988. Hydraulic jump and undular bore formation on a shelf break. J. Phys. Oceanogr. 18 (7), 947–962.

Stansell, P., MacFarlane, C., 2002. Experimental investigation of wave breaking criteria based on wave phase speeds. J. Phys. Oceanogr. 32 (5), 1269–1283.

Touboul, J., Banner, M.L., 2021. On the breaking inception of unsteady water wave packets evolving in the presence of constant vorticity. J. Fluid Mech. 915, A16.

Whitham, G.B., 1962. Mass, momentum and energy flux in water waves. J. Fluid Mech. 12 (1), 135–147.

Whitham, G.B., 1965. Non-linear dispersive waves. Proc. R. Soc. Lond. Ser. A 283 (1393), 238–261.

Whitham, G.B., 1974. Linear and Nonlinear Waves. John Wiley & Sons.

Wu, C.H., Nepf, H.M., 2002. Breaking criteria and energy losses for three-dimensional wave breaking. J. Geophys. Res. 107 (C10), 41–1–41–18.

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