Spatial Evolution of Skewness and Kurtosis of Unidirectional Extreme Waves Propagating over a Sloping Beach

: The understanding of the occurrence of extreme waves is crucial to simulate the growth of waves in coastal regions. Laboratory experiments were performed to study the spatial evolution of the statistics of group-focused waves that have a relatively broad-banded spectra propagating from intermediate water depth to shallow regions. Breaking waves with different spectral types, i.e., spectral bandwidths and wave nonlinearities, were generated in a wave flume using the dispersive focusing technique. The non-Gaussian behavior of the considered wave trains was demonstrated by the means of the skewness and kurtosis parameters estimated from a time series and was compared with the second-order theory. The skewness and kurtosis parameters were found to have an increasing trend during the focusing process. During both the downstream wave breaking and defocusing process, the wave train dispersed again and became less steep. As a result, both skewness and kurtosis almost returned to their initial values. This behavior is clearer for narrower wave train spectra. Additionally, the learning algorithm multilayer perceptron (MLP) was used to predict the spatial evolution of kurtosis. The predicted results are in satisfactory agreement with experimental findings. A machine learning MLP algorithm was used in order to predict the spatial evolution of kurtosis. We found that this MLP algorithm was able to identify patterns and to repro-duce the spatial evolution of kurtosis in a satisfactory manner. This study shows that the MLP algorithm is a promising tool for improving wave forecasting in field measurements and that it can accelerate kurtosis spatial evolution forecasting while retaining good pre-dictive accuracy. In the near future, we plan to perform new tests involving lower and higher wave steepness in order to extend simulations to wider wave train spectra and to improve the prediction of the spatial evolution of extreme wave statistics using the MLP machine learning algorithm.


Introduction
Extreme waves represent a serious problem in coastal engineering. An accurate description of the statistics of free surface elevation can improve our understanding of the propagation of large, steep waves in shallow regions and, consequently, enhance safety in coastal areas. Bitner [1] was a pioneer in the investigation of the statistics of irregular waves propagating toward a decreasing depth via field experiments. There are several mechanisms that can lead to the formation of extreme waves and, subsequently, to different shapes in the probability of the occurrence of wave heights. These mechanisms are extensively discussed in [2], while the topic of our paper is restricted to one of the wellknown mechanisms of extreme wave generation, which is dispersive focusing [3][4][5][6]. The identification of a wave train, as a freak wave train, is made according to a criterion, and the one given by [7], who proposed the amplification index, AI = Hmax/Hs > 2, is adopted in the present study. The parameter Hs represents the significant wave height and Hmax is the maximum wave height, both calculated from the time series.
Using full-scale data is the best approach to investigating the statistics related to extreme waves, but the availability of such data is limited, and these data cannot be used to study the direct effect of changing a specific parameter. Numerical models often complement and expand laboratory or field studies. Despite the recent advances in the computational power in CFD models (e.g., Open-source OpenFOAM), simpler models with more computationally efficient solvers, such as FUNWAVE 2.0 used in [8,9] (based on Boussinesq approximations), are better suited to gather data from various sea states propagating in various topographies. Thus, controlled physical tests will complement findings related to the context of extreme wave formation and might also form a good basis for testing numerical models. Wave characteristics in wave flumes have been extensively studied [10][11][12][13][14]. However, wave properties of extreme waves in shallow waters still require further verification, since the propagation of steep travelling wave trains include complex physical processes, such as energy transfer, nonlinear wave-wave interactions and wave breaking [15,16]. In a Gaussian sea state, the statistics of the free surface elevation can be described by its mean and variance. When waves propagate in shallow water regions and coastal areas, the sea state becomes non-Gaussian and nonlinearities are stronger. Consequently, higher-order moments can be used in order to characterize nonlinear effects. Third-and fourth-order moments, i.e., skewness λ3 and kurtosis λ4, are used to measure the deviation from the Gaussian process [17]. Quantitatively, λ3 characterizes the vertical asymmetry of the wave profile, and λ4 reflects the peakedness and the increase in the crest to trough amplitude.
The two form parameters, skewness and kurtosis, were investigated in several studies in the case of irregular waves [17][18][19][20]. Kashima et al. [20] studied the spatial evolution of these two parameters in the case of JONSWAP random waves propagating from deep water to shallow water regions. They showed that there is no link between the skewness and kurtosis in deep water. However, [21] found that kurtosis depends on the square of skewness (i.e., = 3 + 3 ). Bitner-Gregersen and Gramstad [22] proposed a spatial relationship (a second-order polynomial) between the kurtosis and skewness in the case of unidirectional JONSWAP deep-water random waves simulated using the Higher-Order Spectral Method (HOSM). In shallow regions, [20] demonstrated experimentally that there is a clear link between the two parameters regardless of the incident wave conditions. When the JONSWAP wave propagates in shallow water regions, both parameters increase and reach their maximum magnitude just prior to breaking. These results are consistent with field observations made by [1]. Recently, a similar experimental study was carried out in [23]. the authors studied the spatial evolution of the two form parameters in the case of long-crested JONSWAP random waves propagating over a shoal. They found that the surface elevation has a local minimum of skewness some distance into the down slope of the lee side of the shoal. However, a local maximum was found some distance inside the shallower side of the shoal. The authors also demonstrated that the locations of the maxima of skewness and kurtosis seemingly coincide. Moreover, experimental [24] and numerical studies [25,26] have investigated the spatial evolution of the skewness and kurtosis in the case of constant wave steepness (CWS) wave trains( [24]) and obtained qualitatively similar results to those found in [20,23].
Skewness and kurtosis are of great importance for navigation and coastal applications. The prediction of these parameters can be obtained numerically [27,28]. However, most of the existing wave models are used in the open seas (i.e., deep-water locations) and exhibit a low resolution in the nearshore zone. Considering the rapid advancement in the wave forecasting field, machine learning (ML) comes into play to offer a wealth of techniques to extract information from data that can then be translated into knowledge [29][30][31]. Among various machine learning algorithms, the Multilayer perceptron (MLP) algorithm [32,33] is selected in this study for kurtosis forecasting. MLP is a supervised learning algorithm that provides powerful information processing based on nonlinear regression by optimizing the squared error. The spatial evolution of the kurtosis parameter is a nonlinear problem; thus, linear regression is not suitable for this application, and it is for that reason that MLP was selected. Input datasets of the MLP included wave nonlinearity, abscissa along the flume and experimental kurtosis values. The framework is based on the Python library Scikit-Learn version 0.21.3.

Materials and Methods
Details of the experiments are presented in [15]. However, for the completeness of this study, a brief introduction of the generation of breaking wave trains and surface elevation measurements is provided below. The experiments were conducted in a two-dimensional wave flume at the Laboratory of M2C, Caen. The wave flume was 20 m long, 0.8 m wide, and filled with tap water to a depth of h0 = 0.3 m. A piston type wavemaker was located at one end of the flume to generate the wave trains, and an absorber beach was placed at the opposite end to help damp the incident waves. Temporal variation in surface elevations at desired locations along the wave flume was recorded by wave gauges, along with a high-speed camera. Two wave gauges were arranged to measure the free surface elevations at fifty pre-setting locations along the wave flume. The first wave gauge WG1 was located 4 m away from the wavemaker in order to measure the input wave parameters. The sampling rate for the wave gauges was 50 Hz. In our experiments, only one wave group was generated in each test, and the duration of the data acquisition was 35 s. Temporal surface elevation measurements were used to examine the maximum wave height Hm, the significant wave height Hs, the spectral bandwidth , the skewness λ3 and the kurtosis λ4 along the flat and the sloping beach, which had an angle of α = 1/25 ( Figure 1).
Preliminary tests were made in order to investigate the repeatability. Surface elevation measurements at the same location were conducted. Before breaking, differences were less than 2% (in order of error of the measurements). However, tests are less repeatable once the wave train breaks due to the entrapped compressed air present during the breaking process. Consequently, we are confident that our experiments are sufficiently repeatable. The wave trains were generated by imposing an input wave spectrum at the wavemaker. The energy distribution in the frequency domain is defined by the JONSWAP (γ = 3.3 and γ = 7) and the Pierson-Moskowitz spectra with peak frequency varying between fp = 0.66 Hz and fp = 0.75 Hz, which means that kph0 ranges between 0.8 and 0.93. This implies that the wave trains start propagating in intermediate water depths (kph0 > 0.5). The steepness parameter = ∑ , the same as the local wave steepness adopted in [15], and a modified version of the Benjamin-Feir Index (BFI) (Equation (1)), which measures the ratio of the wave steepness to the spectral bandwidth in a finite water depth [34], are used to characterize wave trains.

= √2
(1) Here, kp is the wavenumber related to the pic frequency fp and calculated using the linear finite water depth relation; ∑ a is the surface elevation at the focusing location according to linear wave theory; ks0 is the spectrally weighted wavenumber calculated at x = 4 m from the wavemaker; and is the spectral bandwidth, which is calculated using the definition given by [35] limited to wave groups with narrow-banded random seastates: Here, mi is the i th spectral moment, calculated as follows: where S(f) is the frequency spectrum. Longuet-Higgins [35] recommended imposing low and high cut-off frequencies in the case of broad-banded wave trains to calculate the i th spectral moment. In this study, fmin = 0.3Hz and fmax = 3Hz are, respectively, the lower and the upper cut-off frequencies. Figure 2 shows the spatial evolution of the spectral broadening of the studied wave trains, which is mainly related to the increasing trend of the spectral energy in high-frequency components, demonstrated in [15]. Generally, the initial spectral bandwidth is higher than 0.25 Hz; therefore, wave trains are considered as relatively broad-banded waves. Most of the generated wave trains in this study have more than one breaker, and so breaking locations (xb) are presented as intervals. Table 1 presents some key parameters of the generated wave trains, which are categorized via their spectrum type, their nonlinearity S0 and their BFI parameter. Twelve selected tests with different spectra and steepness, taken from [15], are studied. Segments of the measured free surface elevation of WTJ3 at six different locations along the wave flume are shown in Figure 3. Table 1. Wave train key parameters. γ is the peak enhancement factor. It is the ratio of the maximum spectral energy to the maximum of the corresponding Pierson-Moskowitz spectrum. The lower value γ = 3.3 corresponds to the standard JONSWAP formulation, and γ = 7 provides a narrower spectrum.  Skewness, λ3 (Equation (4)), is a statistical parameter that is contributed to primarily by the second-order nonlinear interactions between bound waves, and kurtosis, λ4 (Equation (5)), is a statistical parameter that includes a dynamic component due to the thirdorder interactions between wave components. The latter parameter indicates whether the probability density of a surface elevation's peak is higher or lower than that of a typical Gaussian distribution [25]. The kurtosis increase may indicate the formation of extreme waves [26]. In other words, kurtosis measures the weight of the peak relative to the rest of the distribution. In the case of a Gaussian sea state, λ3 = 0 and λ4 = 0 are normally expected. However, for a non-Gaussian sea state, λ3 > 0 is related to sharper crests and troughs, whereas λ3 < 0 is related to wider crests and troughs. Moreover, λ4 > 0 is related to an increased probability of the occurrence of extreme waves. In this paper, excess kurtosis (Equation (5)) is used instead of kurtosis [36].

Wave Train
In Equations (4) and (5), N is the number of samples used in the calculation. The number of samples should be carefully chosen because this may have an impact on the magnitude of the calculated skewness and kurtosis. No null values of η(t) are included in the calculation of skewness and kurtosis; i.e., only the part of the measurement whose value is greater than one-fiftieth is included. The mean value of sea surface is ~ 0.

Results and Discussions
In this section, the skewness and kurtosis coefficients are experimentally estimated at different wave stations along the flume. Figures 4-6 show the spatial evolution of the two shape parameters. The two vertical solid lines added in these figures represent the breaking zone limits, and the horizontal dashed line depicts the thresholds of a Gaussian sea state. Generally, both kurtosis and skewness deviate substantially from the Gaussian predicted values. These deviations can be attributed to the presence of bound waves, which do not fit in with the linear dispersion relationship [18].
When the wave train propagated on the flat bottom (4 m < x < 9.5 m), skewness λ3 ~ 0.25 remained approximately constant, indicating a slight crest-through asymmetry for all the Pierson-Moskowitz and JONSWAP wave trains. However, kurtosis gradually increased and reached ~3 for both the JONSWAP and the Pierson-Moskowitz wave trains. This result is qualitatively consistent with earlier studies [23].
When the wave train reached the slope (9.5 m < x < xb), kurtosis increased rapidly, and skewness increased slightly until reaching the surf zone, where most of the high-frequency waves finish breaking. The evolution of maxima of kurtosis is shown to depend on wave nonlinearity, S0. In other words, the wave trains with small nonlinearities become more nonlinear during the shoaling process than those with higher nonlinearities, which break in deeper regions before undergoing notable shoaling. The maxima of skewness and kurtosis were reached just prior to breaking (x ~ xb). Tian et al. [24] made the same conclusions in the case of constant steepness in wave trains propagating in deep water conditions. Just prior to breaking, the maximum kurtosis value reached was around λ4 = 7, which is greater than to that found in previous studies, such as in [19,20,37]. This might be explained by the sampling variability effect (See [38]). The duration of each experiment in our study is 35 s, whereas it is 20 min of continuous random JONSWAP waves in [20]. As mentioned in [24], when only one wave train is generated in each experiment, quantitative results found concerning skewness and kurtosis do not concern field measurements of continuous wave groups. In Tian et al. [24], the maximum skewness found was less than 1. The maximum skewness obtained in our study indicates a stronger crest-trough asymmetry (λ3 > 1). This can be explained by the shoaling phenomenon, since all the wave trains generated break on the slope. As mentioned in [20], the increase in λ3 and λ4 is mainly due to the secondorder nonlinear interactions associated with wave shoaling and shallow water effects.
After the breaking (x > xb), the wave train dispersed again and becomes less steep, and both skewness and kurtosis converged to zero. Similar conclusions were made by [1], who found that the maximum values of λ3 and λ4 can be obtained in the vicinity of the plunging breaker location. This may be explained by the fact that the breaking process is accompanied by an important energy dissipation, especially in high-frequency components [15]. The breaking process eliminated the strong crest-trough asymmetries, and, as a result, skewness approached zero. Near the coast (x > 13 m), the wave train nonlinearity increased again when the shoaling of the low-frequency waves became more significant. This is clearly illustrated in Figures 4, 5 and 6, where skewness is slightly greater than 0.  Figure 7 shows the initial and maximum values of , λ3 and λ4 as a function of the BFI parameter. Surface elevation measured at the first wave station (WG1 in Figure 1) was used to calculate the initial BFI value. However, local BFI values were calculated where the maximum spectral bandwidth, skewness and kurtosis were achieved. Generally, we found that the local BFI values were lower than the initial ones, and this is related to the spectrum broadening along the wave train propagation. The initial and maximum spectral bandwidth remained approximately constant with the increase in the Benjamin-Feir Index, which is in accordance with the results found in [24]. Additionally, no clear correlations between the BFI and the maximum values λ3 and λ4 could be identified.  Figure 8 shows the relationship between skewness and kurtosis in intermediate and shallow water depth before and after the toe of the slope. The solid black line (Equation (6)) depicts the second-order nonlinear theory first developed by [39], reflecting the water depth change proposed by [40]. The solid blue line(Equation (7)) represents an adjusted formula introduced by the data observed in the sloping bottom. We mention here that a narrow-banded wave approximation is adopted in the second-order theory, and this theoretical model shows that the kurtosis parameter increases when skewness increases.
The experimental data on the flat and on the sloping bottom have a qualitatively similar tendency. For the data on the flat bottom, i.e., kph0 > 0.8, the results are in qualitative agreement with the theoretical model of [40]. For the data on the sloping bottom, i.e., kph0 < 0.8, kurtosis is underestimated by the theoretical model. As expected, most of the skewness and kurtosis values of the experimental data are underpredicted along the flume if only the second-order theory is considered. Mori et al. [18] found a better agreement between experimental data and theory when the free wave effect was included [41]. It is important to note that the formulations in [18] were obtained using numerical simulations in the case of Gaussian deep-water waves, which explains the differences in field measurements, as pointed out by [22].  [40]. The solid blue line depicts an adjusted formula introduced by the data found in the sloping bottom. Figure 9 shows the relationship between the maximum wave height Hmax/Hs, skewness λ3 and kurtosis λ4, on the flat and on the sloping bottom. The solid black curve depicts the Hmax/Hs empirical formulae proposed by [41], and the dashed line represents the threshold of freak waves, i.e., Hmax/Hs = 2. The ratio Hmax/Hs is nearly flat between 1 and 2.5 according to the incident wave steepness, although kurtosis and skewness rapidly increased under the effects of the second-order nonlinear interactions due to the shoaling phenomenon. Consequently, the dependence between Hmax/Hs and the two form parameters is weak in shallow water locations.
Lastly, we tested the MLP algorithm in order to predict the spatial evolution of the kurtosis parameter. Figure 10 shows that the kurtosis predictions are satisfactory, and the accuracy score is around 0.88 for the Pierson-Moskowitz wave trains, 0.9 for JONSWAP (γ = 3.3) and 0.8 for JONSWAP (γ = 7). To gain a better understanding, more data (i.e., higher and lower nonlinearities) should be investigated. Aside from the intense data requirement, sensitivity tests on hyperparameters, such as the stepsize (α) and the exponential decay rates for momentum estimates (β1, β2), should be performed in order to achieve a more accurate MLP model.

Conclusions and Perspectives
The presence of extreme wave events can contribute to the significant deviations in the Gaussian sea state model and the Rayleigh statistics. The present study aimed to investigate the spatial evolution of maximum wave heights, skewness and kurtosis. Using the dispersive focusing technique, group focused waves derived from three different spectra (Pierson-Moskowitz, JONSWAP (γ = 3.3 and γ = 7)) and with the different breaking intensities were generated in the laboratory flume.
Generally, the skewness and the kurtosis for the JONSWAP and the Pierson-Moskowitz wave trains are qualitatively quite similar. Above the flat bottom (4 m < x < 9.5 m), we found that skewness remained approximately constant; however, kurtosis increased gradually during the wave train propagation over the range of conducted tests. Moreover, the skewness and kurtosis magnitudes along this portion of the flume depended on the wave steepness. This is illustrated in Figures 4-6, where the skewness values are slightly higher in the case of strong nonlinearities, S0. Nevertheless, there was a quick transition from intermediate to shallow water depth. Kurtosis depended on the evolution of skewness regardless of the incident wave steepness, (S0). It is shown that the second-order theory, suggested by [40] (Equation (6)), fits well with the data measured along the flat bottom. When the shoaling process started, the formula clearly underestimated kurtosis, and a modified version was proposed (Equation (7)). Generally, the scatter of the data was large, and this can be partially explained by important wave-wave interactions, which are not taken into account in the second-order theory. Furthermore, we found that the ratio Hmax/Hs was constant, although the two form parameters increased rapidly in shallow regions. Therefore, the link between Hmax/Hs and the two form parameters was very weak in shallow water regions in the case of our experiments.
Due to the randomness of the Pierson-Moskowitz and the JONSWAP spectra, wave statistics will depend on which portion of a record is used in the analysis. In other words, the length of a wave recorded in the investigated wave train is very important. An inherent disadvantage of this study is the problem of sampling variability. This kind of problem can be addressed by increasing the duration of temporal measurements in order to provide accurate form parameters and to quantify the effects of the sampling variability.
A machine learning MLP algorithm was used in order to predict the spatial evolution of kurtosis. We found that this MLP algorithm was able to identify patterns and to reproduce the spatial evolution of kurtosis in a satisfactory manner. This study shows that the MLP algorithm is a promising tool for improving wave forecasting in field measurements and that it can accelerate kurtosis spatial evolution forecasting while retaining good predictive accuracy. In the near future, we plan to perform new tests involving lower and higher wave steepness in order to extend simulations to wider wave train spectra and to improve the prediction of the spatial evolution of extreme wave statistics using the MLP machine learning algorithm.