# Generation of 3-hr Long-Crested Waves of Extreme Sea States With HOS-NWT Solver

Maxime Canard, Guillaume Ducrozet, Benjamin Bouscasse

► To cite this version:

Maxime Canard, Guillaume Ducrozet, Benjamin Bouscasse. Generation of 3-hr Long-Crested Waves of Extreme Sea States With HOS-NWT Solver. ASME 2020 39th International Conference on Ocean, Offshore and Arctic Engineering, Aug 2020, Virtual, United States. 10.1115/OMAE2020-18930. hal-03124367

HAL Id: hal-03124367

https://hal.science/hal-03124367

Submitted on 21 May 2021

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# GENERATION OF 3HR LONG-CRESTED WAVES OF EXTREME SEA STATES WITH HOS-NWT SOLVER

**Maxime Canard**

LHEEA Research Department
Ecole Centrale de Nantes, CNRS
Nantes, 44321
France

**Guillaume Ducrozet***

LHEEA Research Department
Ecole Centrale de Nantes, CNRS
Nantes, 44321
France

**Benjamin Bouscasse**

LHEEA Research Department
Ecole Centrale de Nantes, CNRS
Nantes, 44321
France

Email: guillaume.ducrozet@ec-nantes.fr

## ABSTRACT

The accurate control of wave fields generated in experiments and numerical simulations is of great interest for the ocean engineering community. In the context of wave-structure interactions, the recommended practices of classification societies are indeed based on the definition of a wave spectrum, that needs to be reproduced. The present work intends to address this problem from the numerical point of view, using a Numerical Wave Tank equipped with a wavemaker and an absorbing beach, based on the High-Order Spectral method (HOS-NWT). The challenging case of the generation of 3-hours long-crested extreme sea states is studied in details. An iterative procedure to reproduce a target wave spectrum at a given distance from the wavemaker is proposed. The quality of the sea state obtained is evaluated using several criteria defined from spectral quantities. A validation is first performed with a highly nonlinear but non-breaking sea-state. Statistical crest distributions obtained are compared with the Forristall and Huang distributions [1,2]. Then, the Gulf of Mexico 1,000 Year Return Period wave condition is generated. This corresponds to an extreme sea state with significant wave breaking occurrence. The numerical solver needs to be able to account for this phenomenon [3]. The Tian breaking model [4,5] is calibrated to realistically reproduce the dissipation due to breaking, with particular attention paid to the spatial discretization, enlightening its significant effect on breaking model actions. Consequences on the iterative correction process

are studied. The computed statistical quantities appear to be significantly different changing the spatial discretization, while the wave energy spectrum stands the same. It questions the relevance of the characterization of a sea state with the sole wave energy spectrum.

## INTRODUCTION

The design of structures at sea relies on the evaluation of their response facing a set of sea states they may encounter during their life-time. The classification societies define the adequate wave conditions that a specific marine structure should be able to withstand. The observations of rogue waves and their disastrous effects on offshore structures strengthen the need to statistically reproduce real waves conditions in model tests or numerical simulations. The common industry practice for accounting waves is to consider long-crested irregular waves for each design sea state considered and a typical storm duration of three hours. The key in the so-called stochastic approach is to consider a sufficient number of three-hour cases to evaluate correctly the extreme response expected during that duration. This corresponds to emulate enough sets of three-hour waves, to contain the spectral and statistical properties of a specified sea state, representative of the wave condition at sea. Each realization is generated using different sets of random phases inputs. The free surface elevation can be considered as a random quantity and the occurrences of its maxima can be studied, computing the crest

<footnote>*Address all correspondence to this author.</footnote>
heights probability of exceedance (POE). Both POE of single realizations (PDSR) and POE for the ensemble of realizations (PDER) can be studied depending the needs. The ensemble distribution is supposed to be statistically converged, allowing consistent comparisons with benchmark distributions. The control of those quantities stands crucial in the stochastic approach.

Nonetheless, the shape of the crest height POE varies with the nonlinearity of the input conditions. If the linear theory predict a Rayleigh distribution, the second order effects lead to higher crests. The theoretical Tayfun [6] and the parametric Forristall [7] distributions have been built to account for these effects. The latter is usually given as a reference distribution for PDER. As it has been developed using about 1000 three-hours realizations for each sea states, comparisons with PDSR should be carefully performed due to the vague correlation. The mathematical definition using Weibull shape follows. For a given wave crest height $\eta_c$,

$$P(\eta_c > \eta) = \exp\left(-\left[\frac{\eta}{\alpha_r H_s}\right]^\beta\right) \quad (1)$$

For long crested waves, $\alpha_r = 0.3536 + 0.2892S_1 + 0.106U_r$ and $\beta = 2 - 2.1597S_1 + 0.0956U_r^2$. The Ursell number and the integral steepness are defined as $U_r = \frac{H_s}{k_1^2 h^3}$ and $S_1 = \frac{2\pi H_s}{g T_1^2}$ with $T_1$ the mean wave period, $k_1$ the corresponding mean wave number, $h$ the water depth and $H_s$ the significant wave height.

A statistical departure from the second order theory is observed when generating sea states characterized by strong nonlinearity, resulting in higher crests at a given probability level [8,9]. The role of the modulational instability, assessed through the so-called Benjamin Feir Index (BFI) has been studied extensively in past years, focusing on its effect on the probabilities of occurrence of extreme waves, characterized by the kurtosis of the free surface elevation [10]. For long crested waves, narrow-banded and steep sea states lead to the development of high-order effects, increasing the number of extreme waves.

A reference distribution accounting for those phenomena is still needed. As the establishment of a consistent PDER needs a significant number of realizations, a fully experimental benchmark cannot be easily built, leading to the development of several semi empirical distributions using both simulated and experimental data. In [2], using fully nonlinear simulations together with calibration tests, Huang proposed new semi-empirical PDER and PDSR (mean, 99% upper bound, 99% lower bound) formulas of crest height distributions for long crested waves, giving a new benchmark to develop an efficient wave statistical quality-check procedure for wave basin model tests, especially in case of large and steep wave fields. Similarly to the Forristall distribution it uses the Weibull shape. A complete mathematical description of the coefficients can be found in [2].

Historically, the waves used in the design process were generated as the sum of linear contributions and the quality assessment of the generated sea state was performed measuring the spectral quantities. But nonlinear waves simulations are needed to accurately reproduce long-term statistics and improve the coupling with fully-nonlinear Computational Fluid Dynamics (CFD) solvers [11]. In this nonlinear context, normalized procedure and established criteria assessing the quality of the sea state are still missing, especially for very harsh conditions. The need of standardized wave input for CFD solvers is grounded in a more general effort to normalize the entire CFD modeling process for offshore applications. This is the main objective of the currently running Joint Industry Project (JIP) 'Reproducible CFD Modeling Practice for Offshore Applications' [12].

In the present study a particular attention is given to the extreme sea states. Indeed, the accurate reproduction of realistic severe conditions stands a crucial step in a sea keeping or station keeping test. This is awaited that the wave input containing the strongest events are needed when particular interest is put on the extreme responses. Specific procedures need to be developed when those specific wave conditions are generated in numerical or experimental wave tanks. Such harsh conditions exhibit very large steepness $\epsilon = \frac{H_s}{\lambda_p}$, with $\lambda_p$ standing for the peak wavelength. This induces complex physical processes during the wave generation and propagation, such as nonlinear wave interactions or wave breaking. The latter is not naturally solved in nonlinear wave models at use in ocean engineering [13]. However, wave breaking models have been recently implemented in such wave solvers to mimic the energy dissipation of this process, allowing the generation of the strongest wave conditions observed from in-situ measurements [5, 14]. In this context, when generating waves in a wave tank equipped with wavemaker, a spatial evolution of both spectral and statistical quantities is observed, depending of the distance from the wavemaker. They are the result of nonlinear wave interactions as well as wave breaking. An iterative procedure, compensating the latter energy changes has been developed to ensure the reproduction of a given sea state in a target zone.

In the framework of the JIP [12], a new set of criteria have been developed to assess the quality of the sea state reproduction. The present study uses a Numerical Wave Tank (NWT) based on the High Order Spectral method (HOS-NWT) implemented with a dedicated breaking model [5, 14, 15]. The main objective is to test the wave qualification criteria in extreme conditions, generating at a given position the target sea states through an iterative correction process. At first, for validation purpose, a highly nonlinear but non-breaking sea state is generated and analysed, comparing crest height POE with Huang distributions. The influence of changing the target distance from the wavemaker during the iterative process is discussed. The iterative procedure is then tested with a strong breaking case using the 1000-year Return
Period Gulf of Mexico (GOM) condition. A particular attention is given to the influence of breaking models depending of the spatial discretization. The spatial evolution of statistical and spectral quantities is also studied. Several set up leading to the same spectrum but different statistics are tested, questioning the relevance of the characterization of a sea state with the sole wave energy spectrum.

## NUMERICAL WAVE MODEL AND DATA ANALYSIS

This section is dedicated to the presentation of the numerical wave model at use in this study, together with some details relative to the methods used to process the data.

### High-Order Spectral - Numerical Wave Tank

A complete description of the solver can be found in [16]. We consider a 2D rectangular NWT with horizontal dimension $L_x$ and $h$ the water depth. A wavemaker is located on $x=0$. The walls stand perfectly reflective and a numerical absorbing beach is implemented at the end of the domain. The latter is designed so that no reflexion is expected to occur in the simulations. The fluid is assumed homogeneous, incompressible and inviscid. A velocity potential $\phi(x,z,t)$ is introduced verifying the Laplace equation in the entire domain, $\Delta\phi(x,z,t) = 0$. Dynamic and kinematic free surface conditions are expressed introducing the surface velocity potential $\tilde{\phi}(x,t) = \phi(x,z = \eta, t)$, with $\eta$ the free surface elevation. Those two conditions read, on $z = \eta(x,t)$

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

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

$\eta$ and $\tilde{\phi}$ are the two quantities of interest and $\nabla$ stands for the horizontal gradient. Time stepping is achieved solving Eqs. (2) and (3). The only unknown is the vertical velocity on the free surface $\frac{\partial \phi}{\partial z}$. This quantity is computed with the High-Order Spectral (HOS) method, using Taylor expansions and power series developments [17, 18]. The order of nonlinearity $M$ can be specified arbitrarily. Spatial discretization is performed through spectral decomposition, using wave-tank modes.

$$ \eta(x,t) = \sum_{i} \sum_{j} A_{ij}^{\eta}(t) \psi_{ij}(x) \quad (4) $$

$$ \tilde{\phi}(x,t) = \sum_{i} \sum_{j} A_{ij}^{\tilde{\phi}}(t) \psi_{ij}(x) \quad (5) $$

with $\psi_{ij}(x) = \cos(k_i^x x)$ and $k_i^x = \frac{i\pi}{L_x}$. The number of mode $N_x$ controls the spatial discretization.

Note that the original HOS model was dedicated to the propagation of a prescribed initial free surface elevation $\eta(x,t=0)$. The present NWT starts from rest and the waves are generated thanks to a wavemaker with a prescribed motion. An additional potential accounting for wavemaker motions is added, resulting in forcing terms modifying the free surface conditions.

$$ \phi = \phi_{spec} + \phi_{add} \quad (6) $$

$\phi_{add}$ and $\phi_{spec}$ verifying respectively the wavemaker and the free surface boundary conditions. We refer to [16] for more details about the HOS-NWT model, which is released open-source<fnref n="1" />.

### Breaking model

The previous HOS-NWT model is based on nonlinear potential flow formalism and consequently neglects vorticity and viscous phenomena. This enables an efficient and accurate wave propagation model. However, these assumptions prevent the use of such models in the case of breaking waves. To overcome this limitation, a methodology has been recently implemented in HOS-NWT to take into account those breaking events. Note that the breaking models implemented only emulate the wave spectrum modification due to breaking, dissipating energy in the relevant frequency range.

The Tian-Barthelemy model used for this study is detailed in [5, 15]. It is based on the definition of a breaking onset criteria, to detect the individual breaking event, prior to the effective breaking. This criterion compares the water particle velocity at top of the crest $U_x$ to the crest speed $C_x$, following Barthelemy's approach [19]

$$ B_x = \frac{U_x}{C_x} = \frac{\partial \phi}{\partial x} \frac{1}{C_x} > \text{threshold} \quad (7) $$

$C_x$ is computed through a complex procedure using Hilbert transforms. For the present study a threshold value of 0.85 is chosen as advised in [5]. The criterion is computed for each time step in the entire domain.

Once it is verified a dissipative term is added to the free surface (Eqs. (3) and (2)), following the Tian model [4]

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

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

<footnote n="1">¹https://github.com/LHEEA/HOS-NWT</footnote>
The eddy viscosity $v_{eddy}$ is defined as $v_{eddy} = \alpha \frac{H_{BR}L_{BR}}{T_{BR}}$ with $H_{BR}$, $T_{BR}$ and $L_{BR}$ depending on the wave geometry profile just prior to breaking. Their definition can be found in [4], the recommended value being $\alpha = 0.02$. $v_{eddy}$ stands the same on the entire spatio-temporal activation zone defined thanks to $T_{BR}$ and $L_{BR}$.

## Data analysis

Spectral and statistical analysis are performed using the free surface elevation time series obtained with numerical probes at several locations in the NWT. For the sake of brevity, only the main elements of the data analysis are recalled here, noting that it has been carried out following the JIP procedures [12].

Spectra are obtained from the Fourier transform of the temporal free surface elevation. They are computed for each numerical probe and smoothed by convolution, process dependent of the convolution window length. To gain consistency while changing the sea state conditions or the simulation time, smoothing have been adapted to the spectral width controlling the ratio $\delta f/v$ with $v$ the spectral width and $\delta f$ the frequency length of the convolution window. It is fixed for the present study to 0.4. Spectra obtained for each realization are then averaged and the spectral parameters are extracted.

Third- and fourth-order moments of the free surface probability density function (skewness and kurtosis) are computed along the wave-tank. Those are widely used to characterize the nonlinearity of the wave field. Skewness is an indicator of vertical asymmetry (dominated by second-order). Kurtosis is classical in the study of rogue waves as a simple measurement of the occurrence of unexpected large waves in a sea state. PDSR and PDER of the crest height are also computed at the target location. A crest height is here defined as the maximum free surface elevation between an upward and a downward zero crossing. The crest heights values are listed and sorted allowing POE computation. A 95% Jeffrey interval is computed for each value of the POE to evaluate the statistical uncertainty [20]. This is a binomial proportion confidence interval that can be evaluated assuming that for every given $\eta_c$ the number of wave crest verifying $\eta > \eta_c$ follows a binomial law $B(N, POE(\eta_c))$ with N the total number of wave crests. POE obtained are compared with Forristall [7] and Huang [2] distributions described in Introduction.

## CORRECTIVE PROCEDURE & QUALITY ASSESSMENT

The main objective of this study is to generate in a NWT a prescribed long crested wave spectrum at a given distance from the wavemaker. As a reminder, the classical procedure is to deduce the latter temporal motion from the wave spectrum and the linear transfer function of the wavemaker. The wave generation as well as the propagation are subject to nonlinear processes that modify the spectral shape at the target location and needs to be corrected. The most important phenomena are the modulational

instability as well as the wave breaking. Wave qualification criteria must be defined to assess the quality of the sea state after correction. Note that the complete procedure can be applied to experimental wave tanks.

## Iterative correction

For each sea state, the wavemaker motions are corrected through a linear iterative process. To characterize accurately the energy content of the wave field, each iteration requires 10 realizations, varying the random phases inputs. The first iteration uses the target wave spectrum to deduce the wavemaker motion. Then, the Fourier amplitudes of the input are corrected using the spectrum obtained at the target location. For each amplitude $A_{input}(f)$ a linear coefficient $C^1(f) = \frac{A_{target}(f)}{A_{simu}^1(f)}$ is applied (with $A_{target}(f)$ and $A_{simu}^1(f)$ the Fourier amplitude of the target and simulated spectrum). The correction is performed iteratively. The corrective coefficient of the step $n+1$ is $C^{n+1}(f) = \frac{A_{target}(f)}{A_{simu}^n(f)} C^n(f)$. To avoid noise, corrective coefficients are smoothed by convolution. The process ends when the simulated spectrum at target location matches the target sea state, as detailed in the next paragraph. Once the sea state is converged (using 10 realizations per iteration), the final data-set for this configuration is obtained with a total of 100 realizations.

## Wave qualification

The wave qualification criteria used in the present article have been developed in [12]. Both spectral and statistical quantities are checked using averaged values (mean on all 100 realizations).

The main criteria concerns the difference between the target and the simulated spectra obtained at the calibration wave gauge. This difference must be within 5% for the frequencies in the range $f \in [\frac{3}{4}f_p; \frac{3}{2}f_p]$ with $f_p$ the peak frequency of the corresponding spectrum. This frequency range is represented in Figure 1 between the 2 vertical dashed lines. This criteria is used as convergence criteria for the iterative correction process described in the previous section.

Complementary, the crest height distributions are checked at target location. PDSR crest distributions should lie between the 99% lower and upper bounds of the Huang distribution while the PDER should agree with the Huang ensemble distribution within 5% for POE values higher than $10^{-3}$.

Spatial stability of the sea states is also evaluated close to the target location. In this region, the spatial variations of $H_s$, peak period ($T_p$), kurtosis and skewness should be minimized.
<figure><img src="image_2.png" /><figcaption>FIGURE 1: PART OF THE WAVE SPECTRUM CONSIDERED FOR QUALITY ASSESSMENT.</figcaption></figure>

## PROBLEM SET-UP

### Wave conditions

The unidirectional sea states generated are represented by JONSWAP spectrum [1]

$$S(\omega) = \frac{\alpha_p g^2}{\omega^5} \exp \left[ -\frac{5}{4} \left( \frac{\omega}{\omega_p} \right)^{-4} \right] \gamma^{\exp[-(\omega-\omega_p)^2/(2\sigma^2\omega_p^2)]} \quad (10)$$

with $\omega$ the wave angular frequency, $\alpha_p$ the Phillips parameter, $\omega_p = 2\pi f_p$ the peak wave angular frequency, $\gamma$ the peak enhancement factor, $\sigma = 0.07$ for $\omega < \omega_p$ and $\sigma = 0.09$ for $\omega > \omega_p$.

Two target sea states have been generated. Their characteristic are described in Table 1. The Benjamin-Feir Index (BFI) is defined as $BFI = \mu_m \sqrt{2/\nu}$ with $\mu_m = k_m \sigma$, $k_m$ obtained from $\omega_m = \frac{m_1}{m_0}$, $\sigma$ the standard deviation, and $\nu$ the spectral width. The BFI is an important parameter regarding the stability of a sea state [10].

Regarding the nonlinearity of the wave fields, characterized by the steepness $\epsilon$, Case 1 is the mildest condition. Breaking events should not occur, allowing the verification of the proposed procedure, free from the influence of the breaking model. Case 3 corresponds to the Gulf of Mexico 1000-year Return Period sea state. It is expected that this extreme case involves severe breaking. Several target locations and spatial discretizations have been tested (the target location is controlled by $X_t$ the position from wavemaker). The simulations presented in this article are summarized in Table 2.

## Numerical set-up

The simulations have been run at 1/100 model scale (Froude similitude) to easily compare the results with possible experi-

<table><caption>TABLE 1: SEA STATES CHARACTERISTICS (ORIGINAL FULL SCALE).</caption><thead><tr><th></th><th>Case 1</th><th>Case 3</th></tr></thead><tbody><tr><td>γ</td><td>1.0</td><td>2.6</td></tr><tr><td>H<sub>s</sub></td><td>6m</td><td>17m</td></tr><tr><td>T<sub>p</sub></td><td>12.25s</td><td>15.5s</td></tr><tr><td>λ<sub>p</sub></td><td>234m</td><td>375m</td></tr><tr><td>ε = H<sub>s</sub>/λ<sub>p</sub></td><td>2.5%</td><td>4.5%</td></tr><tr><td>BFI</td><td>0.3</td><td>0.7</td></tr></tbody></table>

<table><caption>TABLE 2: WAVES GENERATED.</caption><thead><tr><th>Case</th><th>N<sub>x</sub> (k<sub>max</sub>/k<sub>p</sub>)</th><th>X<sub>t</sub></th><th>α</th></tr></thead><tbody><tr><td>1</td><td>1024 (25)</td><td>7m (3λ<sub>p</sub>)</td><td>/</td></tr><tr><td>1</td><td>1024 (25)</td><td>15m (6,2λ<sub>p</sub>)</td><td>/</td></tr><tr><td>3</td><td>400 (15)</td><td>7m (1,9λ<sub>p</sub>)</td><td>0,04</td></tr><tr><td>3</td><td>256 (10)</td><td>7m (1,9λ<sub>p</sub>)</td><td>0,04</td></tr><tr><td>3</td><td>256 (10)</td><td>17m (4,5λ<sub>p</sub>)</td><td>0,04</td></tr><tr><td>3</td><td>184 (7)</td><td>7m (1,9λ<sub>p</sub>)</td><td>0,04</td></tr><tr><td>3</td><td>184 (7)</td><td>17m (4,5λ<sub>p</sub>)</td><td>0,04</td></tr></tbody></table>

mental measurements conducted in the ECN facilities. The dimensions of the numerical domain are those of the ECN wave tank, length $L_x = 50 m$, depth $h = 5 m$. ECN wavemaker geometry is also reproduced and numerical probes are placed along the wave tank every meter.

The spatial discretization of the HOS-NWT is controlled by $N_x$, the number of points in the domain. In this pseudo-spectral model, $N_x$ is also linked to the maximum wave number $k_{max} = \frac{N_x \pi}{L_x}$ solved. For irregular sea states, the discretization at use is chosen with respect to the ratio $k_{max}/k_p$, with typical values in the range 20–30 for non-breaking cases. In the results presented here, it varies from 7 to 25 to enlighten its effect on the breaking model.

The order of non linearity $M$ of the HOS method is fixed to 3 in this study. This ensures an accurate and efficient solution of the problem. The simulation time is 3 hours real scale (1080s model scale), corresponding to approximately 750 waves at peak period for the sea states considered. The effective simulation
time is slightly larger to take into account the transitory stage (the NWT starting at rest at $t = 0$). The numerical time step is adjusted during the course of the simulation thanks to the use of an adaptative Runge-Kutta scheme, specifying a tolerance value.

## RESULTS AND DISCUSSION

This section presents an overview of the results obtained for the reproduction in a NWT of two different spectra at different locations in the tank. At first, the case of an irregular sea state without breaking is studied, followed by a particular focus on the case of an extreme sea state with significant breaking.

### Non-breaking case

The iterative procedure detailed in previous section is first applied to the irregular nonlinear sea state described as Case 1 in Tab. 1. This wave field is assumed free from breaking events and two different locations for the reproduction in the wave tank are tested, $X_t = 3\lambda_p$ and $X_t = 6\lambda_p$. The spatial discretization chosen is $k_{max} = 25k_p$ allowing an accurate resolution for frequencies lower than $2.5f_p$.

Independently of the target location the corrective process converges with only one iteration. Figure 2 gathers the spectrum and the associated shape criteria obtained at the target location. As the results stand independent of $X_t$, only the case $X_t = 3\lambda_p$ is presented. The wave qualification criterion on the spectral shape is verified, the simulated spectrum follows the target sea state at the target location within 5% in the relevant frequency range. Even without correction the spectrum deviations stand below 10%. This is a direct consequence of the fact that, i) the sea state considered is not highly unstable (BFI = 0.3) and ii) the propagating distance is limited ($X_t = [3; 6]\lambda_p$). For the two target locations, the corrective coefficients stand close to one in all the frequency range. Subsequently the sea state in the entire domain is nearly independent of the target location. The maximum spectrum deviation in the range $f \in [\frac{3}{4}f_p; \frac{3}{2}f_p]$ computed all along the wave tank is presented in Fig. 3 for the two configurations. The target location is pointed out with the red marker, indicating where the iterative procedure is conducted. As expected, this is the location of the minimal value of the criteria. It is worth noting that the spectrum shape is correct over the entire computational domain, except for a short transition zone between $x = 0$ and $x = 2\lambda_p$ probably altered by the presence of wavemaker evanescent modes. Consequently the spectral parameters ($H_s, T_p$) are stable along the wave tank, as required by the wave qualification criteria.

The 100 3-hours single realizations crest height distributions at target location are presented on top of Fig. 4 for $X_t = 3\lambda_p$. Large variability of the statistics are clearly noticed, even for simulations of significant duration (3-hours). The vague correlation with ensemble benchmark realizations such as Forristall

<figure><img src="image_3.png" /><figcaption>FIGURE 2: CASE 1 - SPECTRUM SHAPE VERIFICATION - $X_t = 3\lambda_p$.</figcaption></figure>

stands evident. The 99% bounds built by Huang [2] appear as appropriate analysis tools, surrounding most of single realizations. The ensemble crest distribution at target location are compared with Forristall and Huang POER distributions in Fig. 4 and 5 for target locations $X_t = 3\lambda_p$ and $6\lambda_p$ respectively. As expected a departure from the second order theory is observed for POE values below $10^{-1}$, following Huang PDER distribution. Strong events are then underestimated by the Forristall distribution, testifying the presence of nonlinear processes that could be due to the modulational instabilities (BFI = 0.3). Nonetheless
<figure><img src="image_4.png" /><figcaption>FIGURE 3: CASE 1 SPECRTUM MAXIMUM DEVIATION ALONG THE WAVETANK.</figcaption></figure>

shapes of probability tails diverge from Huang distribution for POE values below $3.10^{-3}$. This phenomenon, observed for the two target locations tested, is stronger for $X_t = 6\lambda_p$. In this configuration, the Huang semi-empirical distribution underestimates the POE obtained with direct simulation in a wave tank configuration.

Complementary to the study on probability distributions, Fig. 6 presents the spatial evolution of the skewness and kurtosis of the free surface elevation. The kurtosis evolution along the wave tank shows an increasing transition zone before reaching a standstill after $x = 5\lambda_p$. Skewness, measuring in some area second order effects on statistics, stands similar in the entire domain. Those are expected results in the context of nonlinear waves [10] which are not influenced by the location of the target in this configuration.

Even if from a spectral point of view the generated sea state reach the wave qualification criteria, the control of statistical quantities appears to be more complex. Their ergodicity in space, assumed building reference distributions [2] is not verified. Variation are noticed along the wave-tank. Further studies need to be conducted in this respect. The spatial evolution of the statistics could explain the differences observed between the Huang distribution and the POE obtained.

## Breaking case

In the present section, the case 3 of Tab. 1 is generated at $X_t = 1.9\lambda_p$ ($X = 700m$ real scale) using Tian model to account for breaking waves [4]. The spatial discretization taken is $k_{max}/k_p = 15$ which stands relatively coarse. Note that the value used for case 1 was $k_{max}/k_p = 25$. In this extreme condition the robustness of the breaking model is not perfect,

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

<figure><img src="image_6.png" /><figcaption>FIGURE 4: CASE 1 - DISTRIBUTIONS OF WAVE CRESTS AT TARGET LOCATION $X_t = 3\lambda_p$.</figcaption></figure>

preventing the use of higher spatial discretizations. Indeed, as a consequence of the spectral method used in HOS-NWT, the high frequency terms of Eqs. (2) and (3) are subject to numerical instabilities enhanced by i) the maximum wave number $k_{max}$ and ii) the nonlinearity of the wave field. The Tian breaking model has been built to overcome this phenomenon. Nonetheless a lack of robustness is observed using an accurate spatial resolution together with strong breaking conditions. Future works will be conducted to avoid this issue. The eddy viscosity coefficient (Eq. (9)) is fixed to match the decrease of significant wave height with respect to the propagating distance with some reference experimental tests [21]. In order to test the sensitivity of the procedure, numerical configurations using $k_{max}/k_p = 10$
<figure><img src="image_7.png" /><figcaption>FIGURE 5: CASE 1 - ENSEMBLE OF WAVE CRESTS AT TARGET LOCATION $X_t = 6\lambda_p$.</figcaption></figure>

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

<figure><img src="image_9.png" /><figcaption>FIGURE 6: CASE 1 STATISTICAL QUANTITIES ALONG WAVE TANK.</figcaption></figure>

and $k_{max}/k_p = 7$ have been generated. Note that such coarse discretizations strongly decrease the quality of the simulated physics, the numerical solution not being converged.

After 4 iterations the spectra generated at target location reach the target sea state. Deviations with the latter are below 2.5% in the relevant frequency range, verifying the associated criterion. Figure 7 presents the target spectrum together with

<figure><img src="image_10.png" /><figcaption>FIGURE 7: CASE 3 - SPECTRA AT TARGET LOCATION $X_t = 1.9\lambda_p$, $k_{max} = 15k_p$.</figcaption></figure>

the one obtained without correction and the final one for the reference discretization $k_{max} = 15k_p$. It is obvious that the initial spectrum suffers from a great lack of energy, as a result of the effect of the wave breaking model. However, after 4 iterations, the verification criteria are fulfilled. Note that the procedure also lead to the target spectrum for the coarser discretizations. Being able to control the spectrum shape is consequently independent of the quality of the numerical solution.

Spatial evolution of the spectral quantities is presented in Figs. 8 and 9. For the reference discretization $k_{max} = 15k_p$, the significant wave height strongly decreases along the wave tank (Fig. 8). This spectral evolution stands crucial in the physics of the wave field studied and its realistic reproduction is necessary. As a result, the spatial stability of the wave spectrum cannot be verified and should not be requested in breaking conditions. Maximum deviations from target spectrum presented in Fig. 9 constrain to one $\lambda_p$ the length of the zone where the spectrum shape reaches the demanded criterion. The additional energy generated through the corrective compensation process only balances breaking dissipation at a given position in the NWT. The target extreme wave condition is consequently unstable and can only be reached at a specific location.

The degradation of the discretization appears to minimize the energy dissipation. With $k_{max}/k_p = 7$ the phenomenon is quasi non-existent. Tian breaking model action is clearly altered, reminding that its dissipation rate has been calibrated for the reference discretization $k_{max} = 15k_p$. The spatial discretization naturally filters frequencies larger than $f = \sqrt{\frac{k_{max}}{k_p}} f_p$, modifying the physics of the solution. A significant part of the extreme events leading to breaking are then already erased due to the high energy cut-off. Tian breaking model is consequently less active
<figure><img src="image_11.png" /><figcaption>FIGURE 8: CASE 3 - SIGNIFICANT WAVE HEIGHT ALONG THE WAVETANK - $X_t = 1.9\lambda_p$.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>FIGURE 9: CASE 3 - SPECRTUM MAXIMUM DEVIATION ALONG THE WAVETANK - $X_t = 1.9\lambda_p$.</figcaption></figure>

using coarse discretizations, preventing energy dissipation along the NWT. As a result the spatial stability of the wave spectrum increases with coarse discretizations. The zone verifying targeted spectral wave condition covers nearly $5\lambda_p$ for $k_{max}/k_p = 7$ (Fig. 9). It is important to note that this stability is totally unrealistic. The quality of such a wave field is strongly questioned.

Statistical quantities should also be affected by the extreme wave condition and the induced modulational instability. Fig. 10 presents the PDSR and PDER crest distributions computed at target location, for the reference discretization $k_{max}/k_p = 15$. Per construction the breaking model is active on the most extreme events of the sea state. Consequently crest heights distributions computed at target location stand lower than those observed for non-breaking condition (Case 1), even if steepness and BFI are

larger. A significant part of PDSR lies out of Huang 99% bounds, unsatisfying the associated wave qualification criterion. Complementary the PDER follows Huang distribution for POE values above $10^{-1}$ before reaching the Forristall predictions.

PDER obtained at target location with coarser discretizations are presented in Fig. 11 and 12. Distributions lie above the PDER computed for $k_{max}/k_p = 15$. With the coarser discretization $k_{max}/k_p = 7k_p$, the PDER nearly reaches the -5% Huang distribution, verifying the ensemble distribution criterion. The number of extreme events appears then to increase while $k_{max}/k_p$ is decreasing, testifying again the coarse discretization alienating effect on breaking model actions. Note that those PDER, varying from one configuration to another, are associated with the same wave energy spectrum. The latter is therefore insufficient to fully characterize a sea state.

Figure 13 presents the spatial evolution of the skewness and kurtosis of the free surface elevation along the wave tank for this breaking sea-state. The skewness stands larger than for Case 1, as a result of the larger steepness. This is indeed known as the primary parameter influencing the skewness in an irregular sea state [10]. The evolution of the kurtosis also appears enhanced in this configuration, with larger departure from Gaussian theory at the end of the tank. However, the present target location is fixed closed to the wavemaker ($X_t = 1.9\lambda_p$). Crest height PDER computed at $x = 9\lambda_p$ for the discretization $k_{max}/k_p = 15$ stands closer to Huang distribution (cf Fig 14). Note that $x = 9\lambda_p$ is far from the target position. The zone considered implies large spectrum deviation from the target sea state (Fig. 9). Huang distribution is nearly reached but at the expense of the spectrum quality. Larger target location should be generated in order to stand closer to both spectral and statistical criteria. Nonetheless stronger conditions emerge from the larger correction coefficients needed. Consequently the insufficient robustness of the breaking model regarding the induced numerical instabilities prevents the generation of those sea states without decreasing the spatial discretization $k_{max}/k_p = 15$, already relatively coarse. Finally, it is reminded that the previous study inducing non-breaking sea states questioned the possible relevance of the Huang distribution in the present configuration of a NWT. This enlightens the need of experiments able to extract the statistical features of those kind of sea states.

## CONCLUSION

The main objective of the study was to propose a procedure based on an iterative process targeting a wave spectrum at a target location, using spectral and statistical criteria to assess the quality of the sea state. For mild conditions (case 1) the spectra generated easily verify the criteria. Nonetheless, in the present NWT configuration the statistics are not ergodic in space. Depending on the target location the crest height PDER appears to be above Huang benchmark. Therefore the statistical criteria based on this
<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /><figcaption>**FIGURE 10:** CASE 3 - DISTRIBUTIONS OF WAVE CRESTS AT TARGET LOCATION - $X_t = 1.9\lambda_p - k_{max} = 15k_p$.</figcaption></figure>

distribution are not totally satisfied. Controlling the spectrum is here not enough to control the statistics.

For the extreme sea state the procedure stands more complex due to breaking effects. The Tian model is assumed to reproduce realistic energy dissipation along the NWT. But using this model, it is difficult to control the target spectrum, the crest distribution and the target location simultaneously, as demanded by the wave qualification criteria. This can be approximately performed using coarse spatial discretizations. But the dissipation due to breaking is then unrealistic. The relevance of such a set-up is questioned.

The iterative procedure has been successfully implemented to reproduce wave spectrum at target location in the NWT. Nonetheless the control of statistical quantities is not achieved.

<figure><img src="image_15.png" /><figcaption>**FIGURE 11:** CASE 3 - WAVE CRESTS ENSEMBLE DISTRIBUTION AT TARGET LOCATION - $X_t = 1.9\lambda_p - k_{max} = 10k_p$.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**FIGURE 12:** CASE 3 - WAVE CRESTS ENSEMBLE DISTRIBUTION AT TARGET LOCATION - $X_t = 1.9\lambda_p - k_{max} = 7k_p$.</figcaption></figure>

The Huang distribution fails to totally capture the statistical behaviour of the simulated configurations. Consequently, several parts of the process have to be questioned. The statistical ergodicity in space assumed building semi-empirical distributions is not observed. Therefore, the relevance of such benchmarks in the procedure varies with the target location and the nonlinearity of the wave field. The need of extensive experimental results in this configuration of a wave tank is obvious. It appears necessary to assess the relevance of the Huang distribution as well as to validate the numerical simulations (including the breaking model). The latter has demonstrated some limitations associated
<figure><img src="image_17.png" /><figcaption>**FIGURE 13:** CASE 3 STATISTICAL QUANTITIES ALONG WAVETANK - $X_t = 1.9\lambda_p$.</figcaption></figure>

to its robustness. This is also the subject of ongoing work.

## ACKNOWLEDGMENT

This work has been performed in the framework of the Chaire Hydrodynamique et Structure Marines CENTRALE NANTES - BUREAU VERITAS. The input for wave cases and the methodology is inspired by the currently running Joint Industry Project (JIP) 'Reproducible CFD Modeling Practice for Offshore Applications'.

## REFERENCES

[1] Komen, G. J., Cavaleri, L., Donelan, M., Hasselmann, K., Hasselmann, S., and Janssen, P., 1996. "Dynamics and modelling of ocean waves". *Dynamics and Modelling of Ocean Waves*, by GJ Komen and L. Cavaleri and M. Donelan and K. Hasselmann and S. Hasselmann and PAEM Janssen, pp. 554. ISBN 0521577810. Cambridge, UK: Cambridge University Press, August 1996., p. 554.

[2] Huang, Z. J., and Zhang, Y., 2018. "Semi-empirical single realization and ensemble crest distributions of long-crest nonlinear waves". In ASME 2018 37th International Conference on Ocean, Offshore and Arctic Engineering, American Society of Mechanical Engineers Digital Collection.

<figure><img src="image_18.png" /><figcaption>**FIGURE 14:** CASE 3 - WAVE CRESTS ENSEMBLE DISTRIBUTION AT $X=9\lambda_p - X_t = 1.9\lambda_p - k_{max} = 15k_p$.</figcaption></figure>

[3] Baquet, A., Kim, J., and Huang, Z. J., 2017. "Numerical modeling using CFD and potential wave theory for three-hour nonlinear irregular wave simulations". In ASME 2017 36th International Conference on Ocean, Offshore and Arctic Engineering, American Society of Mechanical Engineers Digital Collection.

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

[5] Seiffert, B. R., and Ducrozet, G., 2018. "Simulation of breaking waves using the high-order spectral method with laboratory experiments: wave-breaking energy dissipation". *Ocean Dynamics*, **68**(1), pp. 65-89.

[6] Tayfun, M. A., 1990. "Distribution of large wave heights". *Journal of waterway, port, coastal, and ocean engineering*, **116**(6), pp. 686-707.

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

[8] Petrova, P., Cherneva, Z., and Soares, C. G., 2006. "Distribution of crest heights in sea states with abnormal waves". *Applied Ocean Research*, **28**(4), pp. 235-245.

[9] Onorato, M., Osborne, A. R., Serio, M., Cavaleri, L., Brandini, C., and Stansberg, C., 2006. "Extreme waves, modulational instability and second order theory: wave flume experiments on irregular waves". *European Journal of Mechanics-B/Fluids*, **25**(5), pp. 586-601.
[10] Onorato, M., Cavaleri, L., Fouques, S., Gramstad, O., Janssen, P. A., Monbaliu, J., Osborne, A. R., Pakozdi, C., Serio, M., Stansberg, C., et al., 2009. "Statistical properties of mechanically generated surface gravity waves: a laboratory experiment in a three-dimensional wave basin". *Journal of Fluid Mechanics*, **627**, pp. 235-257.

[11] Choi, Y., Bouscasse, B., Seng, S., Ducrozet, G., Gentaz, L., and Ferrant, P., 2018. "Generation of regular and irregular waves in navier-stokes cfd solvers by matching with the nonlinear potential wave solution at the boundaries". In ASME 2018 37th International Conference on Ocean, Offshore and Arctic Engineering, American Society of Mechanical Engineers Digital Collection.

[12] Workgroup, N. P., 2019. Year 1 report. Tech. rep., JIP on Reproducible CFD Modeling Practice for Offshore Applications.

[13] Ducrozet, G., Bonnefoy, F., and Perignon, Y., 2017. "Appli-cability and limitations of highly non-linear potential flow solvers in the context of water waves". *Ocean Engineering*, **142**, pp. 233-244.

[14] Seiffert, B. R., and Ducrozet, G. "A comparative study of wave breaking models in a high-order spectral model". In ASME 2017 36th International Conference on Ocean, Offshore and Arctic Engineering, American Society of Mechanical Engineers Digital Collection.

[15] Seiffert, B. R., Ducrozet, G., and Bonnefoy, F., 2017. "Simulation of breaking waves using the high-order spectral method with laboratory experiments: Wave-breaking onset". *Ocean Modelling*, **119**, pp. 94-104.

[16] Ducrozet, G., Bonnefoy, F., Le Touzé, D., and Ferrant, P., 2012. "A modified high-order spectral method for wave-maker modeling in a numerical wave tank". *European Journal of Mechanics-B/Fluids*, **34**, pp. 19-34.

[17] Dommermuth, D. G., and Yue, D. K., 1987. "A high-order spectral method for the study of nonlinear gravity waves". *Journal of Fluid Mechanics*, **184**, pp. 267-288.

[18] West, B. J., Brueckner, K. A., Janda, R. S., Milder, D. M., and Milton, R. L., 1987. "A new numerical method for surface hydrodynamics". *Journal of Geophysical Research: Oceans*, **92**(C11), pp. 11803-11824.

[19] Barthelemy, X., Banner, M., Peirson, W., Fedele, F., Al-lis, M., and Dias, F., 2018. "On a unified breaking onset threshold for gravity waves in deep and intermediate depth water". *Journal of Fluid Mechanics*, **841**, pp. 463-488.

[20] Brown, L. D., Cai, T. T., and DasGupta, A., 2001. "Interval estimation for a binomial proportion". *Statistical science*, pp. 101-117.

[21] Derbanne, Q., Menard, O., Darquier, M., and Frechou, D., 2009. "Génération, propagation et dissipation de la houle en bassin d'essai de très grande longueur". In Proc. of ATMA 2019.