Statistics of Extreme Waves in Coastal Waters: Large Scale Experiments and Advanced Numerical Simulations

The formation mechanism of extreme waves in the coastal areas is still an open contemporary problem in fluid mechanics and ocean engineering. Previous studies have shown that the transition of water depth from a deeper to a shallower zone increases the occurrence probability of large waves. Indeed, more efforts are required to improve the understanding of extreme wave statistics variations in such conditions. To achieve this goal, large scale experiments of unidirectional irregular waves propagating over a variable bottom profile considering different transition water depths were performed. The validation of two highly nonlinear numerical models was performed for one representative case. The collected data were examined and interpreted by using spectral or bispectral analysis as well as statistical analysis. The higher probability of occurrence of large waves was confirmed by the statistical distributions built from the measured free surface elevation time series as well as by the local maximum values of skewness and kurtosis around the end of the slope. Strong second-order nonlinear effects were highlighted as waves propagate into the shallower region. A significant amount of wave energy was transmitted to low-frequency modes. Based on the experimental data, we conclude that the formation of extreme waves is mainly related to the second-order effect, which is also responsible for the generation of long waves. It is shown that higher-order nonlinearities are negligible in these sets of experiments. Several existing models for wave height distributions were compared and analysed. It appears that the generalised Boccotti’s distribution can predict the exceedance of large wave heights with good confidence.


Introduction
Extreme wave, also known as freak wave or rogue wave, refers in oceanography to large water wave with crest-to-trough wave height H exceeding twice the significant wave height H s in the wave field, or with wave crest height η c higher than 1.25H s [1].In a Gaussian sea state, wave heights H follow a Rayleigh distribution when the wave field is assumed to be narrow-banded.In such cases, large waves fulfilling the criteria H/H s > 2 are not so unusual, occurring approximately once every 3000 waves.For instance, if the average wave period T ave = 15 s, it implies that the observer could probably encounter one of such waves in 12 h.What makes this particular field of research interesting is the fact that the extreme waves may occur not only in the energetic storm sea state but also in a calm sea state, making these waves outstanding and exceptional compared to the surrounding waves.These extreme waves are supposed to be very rare basing on Rayleigh distribution model, whereas they seem to have a larger occurrence in the real world, as suggested in the following scientific studies and reports [2][3][4][5].
The possible mechanisms of the formation of extreme waves are summarised and discussed in recent reviews [6,7].One of the well-known mechanisms of the generation of freak waves is the so-called Benjamin-Feir (or modulational) instability [8,9], which is frequently studied within the framework of the nonlinear Schrödinger equation (NLS).The NLS equation can describe the evolution of the narrow-banded weakly nonlinear waves [10][11][12][13].More recent reviews show that real sea states are more complex and one needs to consider three-dimensionality [14], dissipation [15][16][17] and breaking [18] to achieve more accurate modeling of such wave fields.For constant water depth, the modulation instability of unidirectional waves should disappear with the relative water depth kh < 1.363 (k denotes the wave number and h is the water depth) [8,19].However, for the multidirectional sea states, this is no longer the case.The three-dimensionality not only affects the modulation instability but also the statistical parameters [20,21]; the sea states are either focusing or defocusing depending on the spreading angle.Although we are aware of the essential importance of short-crest waves, the topic of this paper is restricted to one horizontal propagation direction case since the uneven bottom effects on freak waves are not fully understood yet in such configuration.Recently, it has been emphasised that the formation of extreme waves is more related to the second/third-order non-resonant or bound harmonic waves than modulation instability, especially in the finite water depth case where the instabilities are further attenuated [22].
The evidence of the link between a water depth transition and a higher probability of the occurrence of extreme waves has been shown both experimentally [23,24], and numerically [25][26][27].There are also studies dealing with extreme wave statistics in coastal areas with different shapes of variable bathymetry.Katsardi et al. (2013) [28] conducted experimental tests with unidirectional waves propagating over mild bed slopes (1:100 and 1:250) including breaking zones, and made extensive comparisons of wave height distributions.Nonlinear transformation of irregular waves propagating over sloping bottoms (1:15, 1:30 and 1:45) is discussed using wavelet-based bicoherence in [29].Ma et al. (2014) [30] studied spatial variations of skewness, kurtosis and groupiness factor for irregular waves over a bar (1:20) in shallow water.The obliquity effects on the statistical parameters are discussed in [21] using 3D simulations with a High-Order Spectral (HOS) model.
In this context, the objectives of the present study were three-fold: (1) studying experimentally the statistical distribution of (extreme) wave heights considering variable bottom coastal conditions in a uni-directional water wave flume; (2) comparing the collected temporal wave profiles and corresponding statistical distributions with the results from two advanced numerical models; and (3) comparing the measured and simulated wave height distributions with a number of existing statistical models.Regarding Objective 1, the spatial evolution of unidirectional irregular sea states due to the presence of a constant bottom gradient was studied in the large scale facility of Tainan Hydraulics Laboratory (THL) in Taiwan.This large facility allowed us to investigate the full life cycle of extreme waves, from the generation to the degeneration.Several cases with different experimental conditions were tested, while one representative case was investigated for the validation of numerical nonlinear models.In addition to the overall evolution of the sea state, particular attention was paid to the sloping bottom zone as well as the two depth transition regions.Regarding Objective 2, two accurate numerical models were adopted to simulate the experiments, both based on the fully nonlinear potential flow theory: a higher-order Boussinesq-type model based on the model in [31] and the highly nonlinear and dispersive model, whispers3D, using a spectral approximation of the potential in the vertical direction [32,33].As shown below, the corresponding results indicate that the increase of the probability of occurrence of freak waves was clearly observed in experiments and well captured in the simulations, over a region where the water depth is relatively mildly varying.Both measured and simulated wave profiles were studied by using statistical, spectral and bispectral analysis.Finally, Objective 3 was devoted to the statistical distributions of measured and simulated wave heights comparing with six existing statistical distribution models.
The present article is structured as follows.In Section 2, the experimental set-up and the tested sea state conditions are introduced.The mathematical models and their numerical implementations are presented in Section 3. Several signal processing techniques are adopted to interpret the obtained data: from the measurements as well as the simulations, including spectral analysis (Section 4), bispectral analysis (Section 5), nonlinear wave parameters analysis (Section 6), and statistical analysis of the wave height distribution (Section 7).Conclusively, the main results are summarised in Section 8.

Experimental Set-Up and Test Conditions
The  The first probe, located at x = 3 m, was used as a reference for the input signal generated by the wavemaker, and then 7 additional wave probes were used to measure the wave evolution in the deeper area.Probes 9-12 were used to record the wave evolution over the slope.Several probes were distributed in the shallower region to observe the post-dynamics of extreme waves appearing around the end of the bottom slope.
Unidirectional irregular waves were generated using a linear superposition model using random initial phases and amplitudes as determined from a JONSWAP spectrum, given in Equation ( 1).The experimental condition was controlled by four parameters: water depth h, significant wave height H s , peak frequency f p and the peak enhancement factor γ.
where g denotes the gravitational acceleration; α denotes the adjustment factor to achieve the target H s ; and σ J denotes the spectral asymmetric parameter with respect to (w.r.t.) the spectral peak at f p , with σ J = 0.07 when f < f p and σ J = 0.09 when f > f p .The spectrum was discretised using 16,384 = 2 14 frequency components over the range [0, 6 f p ] to generate irregular sea states.During the experimental campaign, more than 30 cases were tested, which were chosen based on the nonlinearity (wave steepness k p H s , with k p being the wave number corresponding to the peak wave period T p ) and the dispersion (w.r.t. the relative water depth k p h) in the deeper and shallower regions.The relative importance of these two effects was measured by the Ursell number U r ≡ H s L2 p /h 3 .
For each test, the test duration lasted for 1200T p .Among all tests, the experimental results of three representative cases shown in Table 1 have already been studied [34].Only Case 2 was selected and simulated in this study for the following reasons: (1) the sea state of Case 1 is almost Gaussian and the number of measured large waves is insufficient to do reliable statistical analysis; (2) in Cases 2 and 3, the relative water depth k p h < 1.363 is verified for the whole bottom profile, which means the effect of modulation instability can be ignored; and (3) Cases 2 and 3 show similar behaviour in both the statistical and spectral analyses.Case 2 is also preferred because fewer wave breaking events happened during the experiment, and they are not accounted for in the two numerical models.Regarding the wave parameters, one difference between the work of Trulsen et al. [23] and the present work is that the nonlinearity of the waves in our case is significantly increased by the shoaling over the slope (note that k p H s ≈ 0.197 in our case is larger than 2k p a c ≈ 0.14 in the shallower side of Case 1 in [23]) with the dispersion effects, i.e., k p h and k p h being similar.Regarding the wave tank, our configuration was shorter in the deeper part but longer in the shallower part, also the slope is longer with the normalised length k p L s ≈ 18 (L s denotes the length of slope) than in [23] where k p L s ≈ 8. Finally, a larger number of wave gauges were set in the flume to follow the wave evolution in space.

General Modeling Approach
The wave problem is formulated in a 2D vertical Cartesian coordinate system (x, z), with x-axis located at the mean free surface level, and z-axis positive upwards.In the framework of potential flow theory, the fluid is assumed to be homogeneous and inviscid and the flow is assumed irrotational so that the wave breaking cannot be directly considered.With the velocity potential φ(x, z, t) introduced, the nonlinear potential flow problem is: (2) where the subscripts denote partial derivatives (i.e., φ x ≡ ∂φ ∂x ).The free surface boundary conditions (FSBCs) Equations ( 3) and (4) are expressed as functions of free surface variables η(x, t) and φ(x, t) ≡ φ(x, z = η, t), forming the so-called Zakharov equations: where the vertical velocity at the free surface is w(x, t) ≡ φ z (x, z = η(x, t), t), which is needed to integrate the Zakharov equations in time.Determining w from the free surface variables η and φ is known as a Dirichlet-to-Neumann (DtN) problem.It should be noticed that the wave flume is large so that the frictional dissipation is significant during the experiments.This is modeled here by adopting two additional viscous terms in Equations ( 6) and (7), following the work of Dias et al. [35].The value of ν needs to be calibrated: in our simulations, the value of ν = 0.0006 m 2 /s was chosen as it allows reproducing the mean decay rate of wave energy as waves propagate along the wave flume.

Boussinesq-Type Model
The Boussinesq-type model for the fully nonlinear and dispersive water waves with potential formulations [31] was used to simulate the experiments.In the model, the velocity potential φ is expanded around a certain depth ẑ(x) in the water column using a power series of the vertical coordinate z: where N B is related to the order of the model, and, in our simulations, N B = 2 ws chosen.φ(n) ≡ ∂ n φ/∂z n | z= ẑ for n = 0, 1, 2, 3, ..., ∞, with ẑ is usually fixed at the middle of the water column to optimise the dispersion properties.We define ŵ ≡ φ(1) as the vertical velocity at the chosen elevation.
The in-depth study of this model can be found in the literature [31,36,37].The main steps are briefly recalled here (see [31] for more information).Introducing Equation ( 8) into the Laplace equation (Equation (2)) and assuming ẑ is a function of slow variable δx (with δ 1), the velocity potential φ can be formulated as an expression of φ(0) and ŵ.In the N B = 2 model, the expression of φ is truncated at order O(δ 2 ) and the higher-order derivatives are up to 4N B + 1.The accuracy of the truncation can be significantly improved by adopting the enhancement technique.New expansion variables are obtained by applying the L-operator: where , with ∇ ≡ ∂/∂x being the horizontal derivative operator.The coefficients of a linear operator L 0 ( ẑ∇) are computed to let the high-order derivatives from 2N B + 2 to 4N B + 1 vanish.The shoaling enhancement operators L 1 ( ẑ∇) and L 2 ( ẑ∇) are used to optimise the linearised shoaling behaviour of the model.The expressions of these L-operators can be found in [31].Now, the velocity potential φ is expressed as a function of φ * and ŵ * via Equation (9).The expressions of the potential at the free surface and at the bottom can be derived using the chain rule.With the Dirichlet boundary condition at the free surface φ| z=η = φ and the bottom boundary condition in Equation (5) expressed in terms of φ * and ŵ * , a linear system is established.The potential φ can be obtained from the solutions of the linear system, φ * and ŵ * .Finally, the vertical velocity at the free surface w is computed to advance the model in time.
The numerical solution procedure is based on the finite difference method for the spatial derivative and the explicit fourth-order Runge-Kutta scheme for the time integration.In the horizontal direction, a stencil of seven points is used to apply up to fifth-order derivative operators.In the linearised N B = 2 model, the dispersion relationship is given explicitly by: where C N B denotes the phase velocity.The coefficients D 2n and E 2n depend on the expansion level ẑ(x)/h(x), and are given by [36] (Equation ( 25)).To evaluate the dispersion property of this model, Equation (10) has to be compared to the exact Airy phase velocity in flat bottom condition C 2 Airy /(gh) = tanh(kh)/(kh).

Whispers3D Model
The modeling approach of whispers3D is presented in previous works [32,33] and summarised hereafter.First, a change of the vertical coordinate from z ∈ [−h(x), η(x, t)] to s ∈ [−1, 1] is applied, mapping the varying domain to a rectangular one: The nonlinear potential system in Equations ( 2)-( 5) can be reformulated using in the new (x, s, t)-space, with φ(x, z, t) = ϕ(x, s(x, z, t), t).Following Tian and Sato (2008) [38], a spectral approach is used in the vertical to approximate the velocity potential.Using the set of orthogonal Chebyshev polynomials of the first kind, denoted T n (s), n = 0, 1, ..., N T , with s ∈ [−1, 1], as an expansion basis, the potential is approximated as: where the a n (x, t) coefficients of the T n terms are now the main unknowns of the problem.Then, the approximation in Equation ( 12) of ϕ is inserted into the governing equations composed of the Laplace equation, a Dirichlet free surface boundary condition (FSBC) with ϕ N T | s=1 ≈ φ and the bottom boundary condition expressed in the (x, s) coordinate system.The so-called Chebyshev-tau method, a variant of the Galerkin method, is used to project the Laplace equation onto the T n polynomials for n = 0, 1, ..., N T − 2 eliminating the s coordinate and giving a set of N T − 1 equations on the a n coefficients at each location x.These N T − 1 equations are supplemented with the Dirichlet FSBC and the bottom boundary condition, so that a system of N T + 1 linear equations with N T + 1 unknowns (a n , n = 0, ..., N T ) at each abscissa is formed.The solutions of this linear system are the a n coefficients, from which the free surface vertical velocity at time t is obtained as: In the whispers3D model, the first-and second-order horizontal derivatives are approximated using fourth-order finite difference formulas with stencils of five nodes.An explicit third-order Runge-Kutta scheme (SSP-RK3) is used for time marching.The dispersion relation of the linear model in constant water depth was derived analytically by Benoit et al. [39]: where the coefficients P n and Q n for different values of N T can be found in Table 1 of [39].

Simulations of the Experimental Case 2
The whispers3D model solves the fully nonlinear potential flow problem, with no assumptions made regarding the steepness of the surface waves or the bathymetry.In the Boussinesq-type model, with the linear operators applied, the model is "improved" in the dispersion relation accuracy.The rapidly varying bathymetry can be well simulated.The dispersion relations of two models are compared to C Airy in Figure 2 over a wide range of relative water depth kh.It is shown that the relative error varies significantly for different water depths.In intermediate and shallow water, i.e., kh ∈ (0, π], the relative errors of C N B and C N T are small (for instance, below 10 −5 for N T = 7 or N B = 2).However, the dispersion relation is less accurate for deep water waves in both models.One advantage of the whispers3D model is that the user can adjust the value of N T easily.In shallow water, a small N T can improve the efficiency of the model while keeping good accuracy.In deep water, by using large N T , the model has the possibility to achieve better accuracy for a large value of kh.The dispersion relations of both models are very accurate under the experimental condition (k p h ≈ 1.2 for Case 2), thus N T = 7 was chosen in the simulations with whispers3D.In the simulations with Boussinesq-type model, the expansion elevation was chosen at ẑ = −0.5h, the shoaling enhancement factors were chosen to be accurate for kh ∈ (0, 30]; the optimised coefficients are given in Equation ( 52) of [31].The chosen time step is ∆t ≈ 0.016 s and the horizontal grid size is ∆x ≈ 0.114 m, corresponding to the Courant-Friedrichs-Lewy number (CFL = C∆t/∆x, C = L/T) is CFL ≈ 0.417.Waves were generated in a relaxation zone of 3L p long at the left boundary and were absorbed in another of 0.5L p long at the right boundary.This set-up was chosen to achieve a better agreement of the low-frequency domain with the measured data.The measured data at x = 3 m were taken as the target spectrum and reconstructed through linear superposition of harmonic incident wave components.In the whispers3D simulation, N T = 7 was used for the approximation of velocity potential.The simulations were made with a time step of ∆t = 0.02 s and a grid size of ∆x = 0.05 m, resulting in CFL = 1.19.
In Figure 3, the spatial evolution of the largest wave group observed in the experiment is shown; the simulated results of Boussinesq-type model and whispers3D are also superposed for comparison.Note that the threshold of freak wave is considered locally since the sea state undergoes significant change due to the uneven bottom.For this reason, the free surface elevation in the figure is normalised by the local significant wave height.The origin of the two freak waves in Figure 3h is shown in the energetic wave packet generated around 33 min from the start, where the sea state has already been affected by the reflected long wave.Figure 3a-c shows that the wave form remains quite similar, from the formation of the wave packet until it propagates over the middle of the slope.The simulation results are in good agreement with the measurements in this region.As waves enter the shallower region and propagate within a certain distance, corresponding to the region between Probe 13 (Figure 3g) and Probe 17 (Figure 3f), two waves in the packet are amplified and can be identified as freak waves since their crest elevations exceed 1.25 times of the local H s .These two freak waves keep their wave profiles longer in the experiment than in the simulations.On the other hand, both numerical schemes seem to underestimate the group velocity as well as the amplitudes of large waves.

Data Processing: Spectral Analysis
Assuming that the free surface elevation is the sum of a large (infinite) number of statistically independent harmonic waves, each component having a random phase in [−π, π] and a constant positive amplitude, we operated with a Gaussian sea state and the statistical characteristics can be described through simple Fourier analysis.The spectral analysis was applied to both measured and simulated results.The simulated signals were resampled (by means of interpolation) to have the same sample points.A low-pass filter was applied to the measured signals to exclude the undesired high-frequency band which might be contaminated by the electronic noise of the wave probes.The time window where the analysis was carried out was translated for each wave gauge with the group velocity C g ( f p ), which ensured that all analysed signals recorded at different positions roughly started from the same wave event.The spectra were estimated via Welch's method.The overlapping factor was 50%, with which the signals were separated into a number of segments.First, the Hann window was applied to each segment of signal (approximately 80 s), and then it was Fourier-transformed through 2 13 -point fast Fourier transform (FFT), resulting in a high spectral resolution ∆ f = 0.0061 Hz.
In Figure 4, both the overall spatial evolution of the spectrum and the detailed spectra at four specific positions are shown.It was observed that, as waves propagate over the deeper region, the wave spectra are modulated mainly in the low-frequency range ( f < 0.2 Hz), and several low-frequency modes are formed before the bottom slope.The low-frequency part is a result of the reflection of unabsorbed long waves and the excitation of the natural modes of the wave flume.The wavemaker and the damping zone only "absorb" a part of the reflected wave energy.Thus, these low-frequency modes increase gradually during the test.Over the bottom slope, second-order effect gets increasingly significant especially around the end of the slope (see the yellow peak at about 2 f p in Figure 4a for x = 53.5 m and Figure 4b for the corresponding spectrum).After the slope segment, the increase of the second harmonic disappears rapidly.More and more wave energy transfers to the low-frequency part, and, as a result, the spectrum broadens.After some distance, due to the energy transfer the low-frequency peak even exceeds the "original" spectral peak and becomes the "new" highest peak of the spectrum (see, for instance, Figure 4b, Probe 24).Apart from the nonlinear wave-wave interactions, the evolution of the spectrum is also affected by the friction resulting from the boundaries, which play an important role in the dissipation of wave energy.The frictional dissipation works as a low-pass filter, i.e., the dissipation is stronger for the high frequency waves than for the low-frequency waves.The spectra simulated with the Boussinesq-type model are shown in Figure 5.The incident wave spectrum shape was well simulated, as waves propagate in the deeper region, the frequencies in the left side of the spectral peak remain in good agreement with the measurements, but differences appear in the higher frequency domain (see Figure 5b).Potentially, some wave-wave interactions occur between the wavemaker and the first probe generating high frequency components.These high-frequency components are considered as free modes in the simulations due to linear wave making, so that they are dissipated rapidly by friction.Around the end of the slope, the second harmonic clearly appears and is more energetic compared to the measured results.The broadening of the spectrum after the slope is well predicted.In total, the Boussinesq-type model provides a good prediction of spectral evolution, despite the slight overestimation of the nonlinear interaction.The spectra simulated with whispers3D are shown in Figure 6.It is shown that the results of whispers3D are in good agreement with the results of the Boussinesq-type model as well as the measurements.Indeed, this model also overestimates the nonlinear wave-wave interactions around the end of the slope so that second (and also higher) harmonic of the spectral peak receive more energy.

Data Processing: Bispectral Analysis
It is well-known that waves in nature can be strongly nonlinear when they enter in coastal areas.Especially in the case studied here, the nonlinear effects are significant, and a Fourier analysis is insufficient to interpret the data.A powerful tool to study the nonlinear process is the polyspectral analysis, in particular, the Fourier based bispectral analysis technique which has been used to study the triad wave-wave nonlinear interactions and quadratic phase coupling of nonlinear water waves [40,41].The bispectrum B( f 2 , f 1 ) is defined as the two dimensional Fourier transformation of the third-order auto-correlation function of the signal.Equivalently, it can also be defined as [42]: where E[•] is the ensemble average of the triple product of complex Fourier coefficients of two frequencies f 1 and f 2 and their sum f 1 + f 2 .The asterisk indicates complex conjugate.For time series with zero mean, the Fourier coefficients X n come from the following decomposition: This definition is economic regarding the computational effort thanks to the Fast Fourier transform (FFT) algorithm.The bispectrum has symmetry properties, i.e., B( f , which can further simplify the computation. In the following, the Fourier coefficients X n were computed in the same way as introduced in Section 4, except that 2 12 -point FFT was used to get a smoother bispectrum.In Figures 7-9, the imaginary parts of the bispectra are plotted, for both the experimental (upper row) and numerical results (second row for the Boussinesq-type model and lower row for the whispers3D model).The two component frequencies f 1 and f 2 were normalised by the incident peak frequency.
The positive and negative values are represented by the colors in the graphs indicating, respectively, the sum interactions and the difference interactions between the two frequency components.For instance, assume B(0.1, 1) < 0, it means that three wave components 0.1 Hz, 1 Hz and 1.1 Hz participate in the nonlinear interaction and the wave energy is transferred from 1.1 Hz to the 0.1 Hz and 1 Hz components.Similarly, if B(0.1, 1) > 0, the 1.1 Hz component gain energy from the 0.1 Hz and 1 Hz components.
Figure 7 shows that, in general, the wave-wave interaction between f p and 2 f p is weak in the deeper region.In Figure 7a1, two relatively strong interactions happen around B( f p , 0.1 f p ) (red) and B( f p , f p ) (blue), which indicates that the wave energy is transferred from both the low-frequency modes around 0.1 f p and the second harmonic 2 f p to the frequencies around the peak.As waves propagate in the deeper region, in Figure 7b1, it is shown that for the low-frequency domain both sum and difference interactions are present.This pattern indicates that the energy of peak frequency is transferred to both lower and higher frequencies and that the spectrum broadens.This pattern can hardly be seen in the corresponding spectrum due to weak interactions, but as shown below, the spectral broadening is more obvious in the shallower region.In the simulations, this energy transfer is overestimated, resulting in slightly larger long waves.When waves arrive at the end of the deeper region, Figure 7c1 shows that the spot around B( f p , 0.1 f p ) is blue, meaning that the energy transfers mainly from the peak frequency to the low-frequency components and that long waves are generated.In the deeper region, the simulated and measured bispectra are in overall good agreement.1)) of graphs presents measured bispectra, while bispectra obtained from numerical results are presented on the second (panels (2)) and third rows (panels (3)) for the Boussinesq-type and whispers3D models, respectively.All the panels share one colorbar, which is in scale (10 −8 m 3 ).
In Figure 8, the bispectra of waves propagating over the sloping bottom area are shown where strong nonlinear effects are expected.In Figure 8d1, it is observed that the nonlinear effect is present but still moderate.In Figure 8e1, the wave-wave interaction is strong, resulting in an increase of the second harmonic 2 f p and the low-frequency waves around 0.1 f p .In Figure 8f1, not only similar characteristics regarding the second harmonic and the long wave generation are seen, but also the wave-wave interactions between f p and 2 f p coupled to 3 f p is visible.The long waves are generated from both peak frequency and its second harmonic 2 f p (see the blue spots at B( f p , 0.1 f p ) and B(2 f p , 0.1 f p )).In the simulations, the agreement is still good but the intensity (or "magnitude") of nonlinear transfers is slightly overestimated.For instance, the sum interactions of low-frequency modes around B( f p , 0.1 f p ) at all three probes are stronger in the two numerical models than in the experiments.As is shown in Figure 8e2 and e3, the sum interactions at B( f p , 2 f p ) are also a bit stronger in the simulations.In Figure 9, three probes are chosen to illustrate the modulation of waves propagating over the shallower region.More wave components participate in the wave-wave interaction, and the contours become more complex.A negative band ranging from B(0.5 f p , 0.5 f p ) to B(0, f p ) is observed which indicates the strong energy transfer from the spectral peak to the bands with frequencies lower than the peak.This blue band is clear in Figure 9g1, and disappears in Figure 9i1 because, as the spectral peak keeps losing energy, this energy transfer becomes weaker and weaker.This explains partially the spatial evolution of the spectrum, namely the decrease of the "original" spectral peak, the increase of the low-frequency components (which eventually become new spectral peak) over the shallower region, and the increase of the long wave energy.Another possible reason for this spectral change is the different frictional dissipation rates for different frequencies.The friction dissipates more energy of high-frequency waves while low-frequency wave energy is less dissipated.Finally, a quasi-steady state is reached and the spectral shape keeps almost unchanged, as can be observed for the spectral shape after 100 m in Figure 4a.

Data Processing: Nonlinear Wave Parameters
The free surface elevation is regarded as a stationary stochastic process.In a Gaussian sea state, the characteristic of the statistics of the surface motion can be fully described by its mean η = E (η) and variance σ 2 = E (η − η) 2 .As waves propagate into coastal areas, the sea states become non-Gaussian and the nonlinearities affect the wave shape, resulting in sharper and higher wave crests.Higher-order moments can be used to characterise nonlinear effects.Usually, skewness (third-order normalised moment) and kurtosis (fourth-order normalised moment) are considered.They can be computed from time series η: The skewness is related to the bound mode harmonics and the asymmetry of the probability density function (PDF) of η, and the kurtosis is related to the sharpness of the PDF [43].For Gaussian sea states, λ 3 = 0 and λ 4 = 3 are expected.For non-Gaussian sea states, λ 3 > 0 represents asymmetric wave shape (sharper crests and f troughs) and the reverse for λ 3 < 0. λ 4 > 3 indicates an increased probability of occurrence of the largest wave heights [23], but it is not guaranteed that, when λ 4 < 3, the large waves (especially freak waves) will not manifest.
Statistical parameters can also be computed based on the bispectrum.Skewness can be defined as the normalised integral of the real part of the bispectrum [40]; the wave profile asymmetry w.r.t. the vertical axis is related to the normalised integral of the imaginary part of the bispectrum [44]: In Figure 10, the statistical parameters of the simulations as well as the measured results are shown, and the two different computation methods of the skewness are also compared.It can be seen that the two methods give very similar results; the small differences are probably due to the use of Hann function window (even though a correction factor has already been considered) when computing the Fourier coefficients.It can be seen that both the Boussinesq-type and whispers3D models have good agreement over the two flat bottom regions, but clearly overestimate the maximum skewness (triad wave-wave interaction) around the end of the slope.While the generated waves are almost symmetric w.r.t. the horizontal axis, they become skewed as they propagate over the slope.Around the end of the slope, the wave shape is highly asymmetric w.r.t. the horizontal axis, and the symmetry is not fully recovered in the shallower flat bottom.Regarding the asymmetry parameter, the simulated results agree well with the measured results.As can be seen, this parameter is almost zero except over a short region after the slope which implies that the waves generated are almost symmetric w.r.t. the vertical axis.It is also noted that the presence of the current bottom slope has limited effect on the asymmetry parameter.
In Figure 11, the evolution of the kurtosis is shown.It fluctuates around 3 except for the area around the end of the slope, where a local maximum value is observed.This corresponds to the location of Probe 13.The numerical models predict well the overall evolution of kurtosis along the bottom profile.However, they both overestimate the kurtosis around the end of the slope, and this trend is more marked with whispers3D.This is probably because the breaking events in the experiments, which are not considered in the simulations, limited the maximum wave heights.At the end of the slope, the prediction of whispers3D is more non-Gaussian compared to that of the Boussinesq-type model.According to previous studies (e.g., [23]), the local maximum of kurtosis indicates higher probability of the occurrence of extreme waves.This motivates further investigations of wave height statistics in similar, however, more complicated configurations.

Data Processing: Wave Height Distribution
The wave height distribution for the nonlinear shallow water waves has been of interest for oceanographers, as it represents a key input for the design of coastal structures.The understanding of the shallow water wave height distribution is, for the moment, limited not only due to the complexity of coastal hydrodynamics but also due to the complex nearshore environment.The existing distributions are either empirical/semi-empirical or analytical under strong assumptions.For this reason, more experimental data are needed to study the applicability and validity of these models.In the following, some of the most popular distributions are listed and compared with the measurements and the results of the simulations.[45] for narrow-banded linear deep water waves.The wave heights (H) follow the Rayleigh distribution, whose complementary cumulative distribution function (CCDF) is given by: The wave conditions of the current Case 2 obviously violate the aforementioned assumptions.The Rayleigh distribution is adopted here only for the purpose of comparison.Assume for the signals with zero mean, the relation H s ≈ 4σ is fulfilled throughout the domain, then the normalised wave height H/σ = 8 is the threshold for the freak waves.

Glukhovskiy Distribution and Its Modifications (GV91 and GK96)
As waves propagate into coastal areas, the shoaling effect, as well as the depth-induced wave breaking, are expected, which affect the shape of the distribution.Due to the insufficient understanding of wave breaking and high-order non-Gaussian effects, only empirical or semi-empirical models are available.A more general form of the Rayleigh distribution, namely the Weibull distribution, is used to account for the depth-induced wave breaking in shallow water, as initially proposed by Glukhovskiy [46]: where A, H * , and K are the Weibull parameters.which are defined differently in the literature.Van Vledder (1991) [47] (labelled as GV91) modified this model by considering the parameters as: with Γ denoting the gamma function, Γ(x) = ∞ 0 e −t t x−1 dt.Klopman (1996) [48] also proposed a modification (labelled as GK96) based on a fit to laboratory data: As input, the local water depth h and the significant wave height H s are required.Both modified Glukhovskiy's distributions are used and compared hereafter.

Composite Weibull Distribution (CWD)
Another well-known model for wave height distribution in shallow water is the semi-empirical composite Weibull distribution (labelled as CWD here) [49].In this model, the smaller wave heights are Rayleigh distributed, whereas the larger wave heights follow a Weibull distribution.The distribution requires not only h and H s but also the bottom slope m.This model also takes wave breaking into consideration: where H tr = (0.35 + 5.8m) h denotes the transition wave height between the two distributions, and H 1 and H 2 are scale parameters for wave heights and are given in [49] (see more details in Equation (6) or Table 2 in this reference).

Weibull-Generalised Pareto Distribution (WGP)
More recently, another two-part model, the Weibull-generalised Pareto wave height distribution (labelled as WGP), has been introduced [50].The parameters, including H s , local wave number k L (here the wave number k p corresponding to the peak frequency f p is used), and local water depth h are required in this semi-empirical model: where In this model, except for the required information regarding the waves and bathymetry, four additional parameters are needed to be calibrated: ρ for the transition wave height, α for the scale of Pareto distribution, β for the Miche wave breaking height limit, and λ for the shape of Weibull distribution.According to the results of Wu et al. [50], ρ = 1, α = 0.22, β ≈ 0.15, and λ = 1 are used.Unlike CWD, this two-part model is continuous and bounded (0 ≤ H ≤ H max ), and this model has a continuous PDF.

Gram-Charlier Type Model (GBD)
The relationship between the fourth-order moment of free surface elevation and the third-order nonlinearity is studied in [51].The effects resulting from third-order nonlinearity on the statistics of wave height distribution have been taken into consideration by using Edgeworth's form of a Gram-Charlier series [51][52][53][54][55]; this family of wave height distributions is sometimes referred to as modified Edgeworth Rayleigh distribution (MER).They do not seem to be of specific attraction in describing the statistics of coastal wave heights because the corresponding PDF of MER is not always positive (which is not realistic).In addition, these distributions are derived under the assumption of narrow-banded weakly nonlinear waves in the deep water condition, which is violated in nearshore areas.However, they are attractive in predicting abnormal waves because MER models showed good agreement for the tail of the wave height distribution when compared to some laboratory measurements [54].Thus, this model is included in the present paper anyways.
A recent model of this family proposed in [55] is a generalised Boccotti distribution (labelled as GBD) in the form of MER, taking third-order nonlinearity and finite spectral width into consideration.Recall the Boccotti's distribution [56]: where with S(ω) denoting the spectral density function, τ corresponding to the time lag of the global minimum of autocorrelation function R(τ) = E [η(t)η(t + τ)], and m i denoting the ith-order moment of S(ω), m i = ∞ 0 ω i S(ω)dω.Then, the formulation of GBD is: where with η denoting the Hilbert transform of the free surface elevation η, and [m, n] = 0, 2, or 4, but not 3 due to the absence of resonant three-wave interactions.It should be noticed that the CCDF starts with P GBD (0) = 1, which is not realistic, meaning that this model is not suitable to model the bulk of the distribution, but only its tail for large wave heights.

Comparative Evaluation of Wave Height Distributions
The wave height distribution at Probe 8 located at the beginning of the slope is first considered and displayed in Figure 12.It can be seen that the Rayleigh distribution predicts well the measured wave height distribution.The agreement between the simulated and the measured data at this position is good, as shown in Figure 12a.In Figure 12b, the other wave height distribution models are similar to the Rayleigh distribution, differences appear only in the prediction of larger waves with H > 6σ.
The measured results show fewer large waves than predicted by the Rayleigh distribution.This can be related to the work of Forristall (1978) [57] showing that the Rayleigh distribution overpredicts the heights of the largest waves in a deep water sea state.From the previous analysis of the nonlinear parameters skewness and kurtosis, the sea state here is quasi-Gaussian and the nonlinear wave-wave interactions are limited for the deeper region.Thus, it is anticipated that the wave heights follow a Rayleigh distribution.The data show better agreement with the modified Glukhovskiy's methods, GV91 and GK96.As waves propagate near the end of the slope, Figure 13 shows a clear increase of the probability of the large waves.The agreement between the measured and simulated results remains good, and the model predictions by CWD and GBD are in better agreement with other models.The two numerical models show similar results.In the Boussinesq-type model, note that one wave fulfils the criteria of freak waves (recalling that H > 2H s is equivalent to H > 8σ).However, the corresponding wave height is lower (below the above threshold) in the measured time series.At this position, skewness and kurtosis computed from the measured signal and from the simulated signals are in excellent agreement.The skewness increases approximately to 0.2, whereas the kurtosis is still around 3, indicating that second-order nonlinearity plays a more important role here.
In Figure 14, the wave height distribution after the slope (Probe 13), where the skewness and kurtosis reach their local maximum values, is illustrated.The amount of large waves is larger than the expectation of Rayleigh distribution, and more freak waves appear after the slope.At this location, it is also observed that some sporadic spilling breakers occurred, limiting the number and height of the large waves.It is clearly shown in Figure 14b that all the wave height distribution models except GBD underestimate the probability of wave height for H > H s .The GBD model shows nice agreement for the same wave height range, even though the case studied here lies out of its validity domain.Regarding the simulations, the Boussinesq-type model results agree well with the measurements.The wave heights are overestimated with whispers3D, as are the skewness and kurtosis.Then, in Figure 15b, it is shown that the tail of the wave height distribution lowers and the degeneration of the extreme waves starts.This change of wave height distribution is also part of the process of settling down towards a new equilibrium state.However, it was observed in the experiments that some of the generated freak waves will keep their wave form (without breaking) and propagate over a certain distance, even though these waves are quite steep, and propagate in relatively shallow water depth.The prediction of GBD is still good regarding large waves.Other models are lower than the Rayleigh distribution probably because the wave breaking effect is overestimated with these models.

Discussion and Conclusions
To achieve a better understanding of the generation mechanism of extreme waves in coastal areas, large scale experiments at THL and highly accurate numerical simulations with two advanced models were performed.A variable bottom profile was designed to study the effects of a constant bottom slope connecting two flat regions on the occurrence of extreme waves.The experimental conditions were chosen based on two non-dimensional parameters: the wave steepness and the relative water depth.Note that our experimental set-up is different from previous experiments [23,24] regarding the length of slope, the arrangement of the two flat regions, and, most importantly, the nonlinear effects were relatively strong in the experimental cases shown here.Besides, the large scale facility allowed investigating the whole life cycle of the extreme waves.Spectral, bispectral and statistical analyses of the collected data were conducted, through which a more comprehensive understanding of the shoaling process has been drawn.
In line with previous studies [23,[25][26][27], second-and third-order effects appear in the frequency spectrum, and local maximum values of the statistical parameters skewness and kurtosis were observed in the shallower region.It was experimentally observed that the second-order effects rapidly weaken after the end of the slope, as shown Figure 4a.The values of skewness and kurtosis also decrease with the weakening of second-order effect: the skewness decreases to a new level instead of zero due to the asymmetrical wave shape in the shallow water, while the kurtosis decreases back to 3.
The imaginary part of the bispectra was used to demonstrate the phase coupling of wave components and the energy transfer due to nonlinear wave-wave interactions.The nonlinear wave coupling is relatively strong around the end of the slope where mainly second and third harmonics are generated.The generation of third harmonic is small in magnitude, but it can be seen in the variance spectrum.Thus, we speculate that the increase of the occurrence probability of extreme waves is more related to the second-order effects.It was also observed in our experiments that a series of low-frequency waves appear.It is anticipated from the bispectra that long waves are mainly generated by the difference interactions as waves propagate in the flume.The frictional dissipation works as a low-pass filter on the wave spectrum, which has limited effects on these low-frequency modes.
In the deeper region, the sea state is almost Gaussian, and the nonlinear wave-wave interactions are moderate.The wave heights are well described by Rayleigh distribution, except that the latter slightly overpredicts the occurrence probability of the highest waves.The occurrence probability of large waves is significantly increased when waves propagate over the slope.This is due to the fact that the sea state is transforming from one equilibrium state (in the deeper region) to another (in the shallower region).During this process, some of the largest waves can be identified as freak waves and will propagate for a distance while keeping their wave form.We also emphasise that the process of settling down to the new state takes time and a certain distance (a distance-lag in response to the depth transition), thus the shallower region after the slope is also prone to freak waves due to this lag.This was clearly noticeable in our flume set-up, which contains a long flat shallower region.Since the nonlinearities obviously come into play, nonlinear wave parameters are useful in describing the sea state, especially the kurtosis is related to the occurrence probability of large waves.It was confirmed that strong modifications of the distributions of the largest wave are rather well correlated with larger values of the skewness and kurtosis (see Figures 10 and 11).The Rayleigh distribution is no longer suitable to describe such sea states, thus a number of wave height distribution models were introduced and compared.In the coastal areas with intermediate or shallow water depth, the wave height distribution models are either semi-empirical or analytical with strong assumptions.In this study, not only the modified Glukhovskiy distributions (GV91 and GK96) and two-part distributions such as CWD and WGP were considered, but also another deep water wave height distribution, GBD, was introduced.This distribution is known to have good predictions for extreme waves [54].As the non-equilibrium statistics appear and fade out over the slope, the probability for large waves increases then decreases in the same manner.It was seen that no distribution model can predict the measured wave height distribution equally well for all probes, but, in the prediction of the occurrence probability of large waves, the generalised Boccotti distribution works reasonably well, although it was originally designed for weakly nonlinear deep water waves.The listed shallow water wave height distribution models underestimate the wave height distribution mainly because they assume waves break during the shoaling process.
The numerical results of the Boussinesq-type model and the whispers3D model agree well with the experimental data.The results are in better agreement in the deeper region than in the shallower region, especially before the reflected long waves reach the probes.In particular, the agreement deteriorates due to the contamination of reflected long waves.A small number of wave breaking events also affect the agreement after the slope.These effects are insignificant at the beginning, but the disagreements accumulate as waves propagate in the flume, and eventually affect the simulated results in a significant way.The agreement is acceptable for the deeper region and the slope area, which is of most interest.Two models show very similar results as they both adopt fully nonlinear FSBC and are capable of accounting for variable bathymetry.We noticed that both models overestimate the nonlinear interactions around the end of the slope.The possible reasons for this are numerous; we anticipate that this is due to the non-equilibrium statistics provoked by the transition of water depth.The small number of wave breaking events may also play a role in constraining the maximum achievable skewness and kurtosis.As a further study, new unidirectional irregular wave experiments in the same facility, but with a new bottom profile containing two constant slopes will be conducted and studied in the future.Finally, it is recalled that wave obliquity and wave directionality were not considered in this study; the effects associated with multidirectional (three-dimensional) sea states are left for future work.
experiments were conducted in the 2-D wave flume of THL at Tainan, Taiwan.The flume is 200 m long and 2 m wide.The waves were initially generated in the deeper region with depths of h = 1.2 m and 1.3 m.The bathymetry was decomposed into three parts: a 30 m long flat part, a 20 m long part with a constant slope of 1/20 and again a 120 m long flat bottom part.Waves were generated at x = 0 m by a piston-type wave-maker, which was also able to absorb part of the wave energy reflected by the uneven bottom or far-field rip-rap mound.The last 30 m part of the flume worked as an effective wave absorber including a deep water area and a mound of stones placed at the end of the flume.Over the flume length, 30 capacitance gauges were used to measure the wave elevation with a sampling frequency of 100 Hz.The locations of the wave gauges are shown in Figure 1.

Figure 1 .
Figure 1.Sketch of the experimental setup, installed at THL.

Figure 2 .
Figure 2. Comparison of the dispersive effects between the Boussinesq-type model (N B = 2) and the whispers3D model (with different choices of N T ) in log-log scale.The two dash lines divide the domain into shallow water region (kh < π/10), intermediate water depth region, and deep water region (kh > π).

Figure 3 .
Figure 3.Comparison of the free surface elevation normalised by the local significant wave height at different positions (panel (a)-(g)) between the simulated results of the Boussinesq-type model (blue lines), whispers3D (red lines) and the measurements (black dot lines).Only the local relative water depth over the uneven region is marked.The threshold of freak wave basing on crest elevation is shown (dash lines).

Figure 4 .
Figure 4.The spectrum of the experimental data shown in different scales: (a) the spatial evolution of the spectrum in a log-scale colormap; and (b) the spectra of four specific positions.

Figure 5 .
Figure 5. Same plot as Figure 4 for the results of the Boussinesq-type simulation.

Figure 6 .
Figure 6.Same plot as Figure 4 for the results of whispers3D simulation.

Figure 7 .
Figure 7. Contour of the imaginary parts of bispectra of Probes 1, 5, and 8 in the deeper region (k p h = 1.19) (corresponding to panels (a)-(c)).The upper row (panels (1)) of graphs presents measured bispectra, while bispectra obtained from numerical results are presented on the second (panels (2)) and third rows (panels (3)) for the Boussinesq-type and whispers3D models, respectively.All the panels share one colorbar, which is in scale (10 −8 m 3 ).

Figure 8 .
Figure 8. Same plot as Figure 7 for Probes 10, 12, and 13 over the slope (corresponding to panels (d)-(f)), the corresponding relative water depths are k p h = 0.86, 0.54 and 0.50 respectively.

Figure 10 .
Figure10.Overall evolution of the nonlinear statistical parameters, skewness and asymmetry, computed from the time series (dash lines) and the bispectrum (solid lines).

Figure 11 .
Figure 11.Comparison of the overall evolution of kurtosis between the measured and the simulation results of the two numerical models.

7. 1 .
Brief Review of Existing Distribution Models 7.1.1.Rayleigh Distribution study of wave height distribution in the early 1950s starts with the work of Longuet-Higgins (1952)

Figure 12 .
Figure 12.The wave height distributions of measurements, simulations and different models at x 8 = 30 m: (a) the CCDF as a function of normalised wave height; and (b) the deviation from the Rayleigh is shown by the normalised CCDF.The horizontal solid line in (a) as well as the vertical solid line in (b), both at H/σ = 8, represent the usual threshold for freak waves.